Changeset 41831 for branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c

Timestamp:

Oct 11, 2021, 11:33:33 AM (5 years ago)

Author:

eugene

Message:

add chebyshev polynomials to 2D fitting; add chebyshev support functions; use chebyshev polynomials for inverse transformation

File:

• branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c (modified) (13 diffs)

branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c

@@ -36 +36 @@
 #include "psLogMsg.h"
 #include "psPolynomial.h"
+#include "psAbort.h"
 #include "psAssert.h"
 
@@ -203 +204 @@
 
 
+/** This function calculates the appropriate scaling factors needed to normalize the
+ * input vector to the range -1 : +1.  These are stored on the polynomial in the given
+ * direction.
+ */
+bool psChebyshevSetScale (psPolynomial2D* myPoly, const psVector *vec, int dir) {
+
+    psAssert ((dir == 0) || (dir == 1), "invalid direction %d\n", dir);
+
+    // find the min and max of the vector
+    psF64 minValue = NAN;
+    psF64 maxValue = NAN;
+
+    for (int i = 0; i < vec->n; i++) {
+        if (isnan(vec->data.F64[i])) continue;
+        if (isnan(minValue)) { minValue = vec->data.F64[i]; }
+        if (isnan(maxValue)) { maxValue = vec->data.F64[i]; }
+        minValue = PS_MIN(minValue, vec->data.F64[i]);
+        maxValue = PS_MAX(maxValue, vec->data.F64[i]);
+    }
+    if (minValue == maxValue) {
+        psWarning ("insufficient data range to determine scale factors\n");
+        return false;
+    }
+
+    myPoly->scale[dir] = 2.0 / (maxValue - minValue);
+    myPoly->zero[dir]  = 1 - myPoly->scale[dir] * maxValue;
+    return true;
+}
+
+/** This function generates a normalized vector in the range -1 : +1 based on the input
+    vector using the scale factors stored in myPoly in the given direction.
+ */
+psVector *psChebyshevNormVector (const psPolynomial2D* myPoly, const psVector *vec, int dir) {
+
+    psVector *norm = psVectorAlloc (vec->n, PS_TYPE_F64);
+
+    if (vec->type.type == PS_TYPE_F64) {
+        for (int i = 0; i < vec->n; i++) {
+            norm->data.F64[i] = vec->data.F64[i]*myPoly->scale[dir] + myPoly->zero[dir];
+        }
+        return norm;
+    }
+    if (vec->type.type == PS_TYPE_F32) {
+        for (int i = 0; i < vec->n; i++) {
+            norm->data.F64[i] = vec->data.F32[i]*myPoly->scale[dir] + myPoly->zero[dir];
+        }
+        return norm;
+    }
+
+    psError(PS_ERR_UNKNOWN, true, "invalid type for chebyshev polynomial");
+    return NULL;
+}
+
+/** This function generates a vector containing the values of a Chebyshev polynomial of
+    the given order evaluated at the coordinates given by the input vector, i.e., this
+    function returns the vector T^n (x_i) where x_i is the input vector of values and n is
+    the polynomial order.
+ */
+psVector *psChebyshevPolyVector (const psVector *vec, int order) {
+
+    if (order > 9) {
+        psWarning ("Chebyshev orders higher than 9 are not yet coded\n");
+        return NULL;
+    }
+
+    psVector *out = psVectorAlloc (vec->n, PS_TYPE_F64);
+
+    // easy but non-general implementation
+    switch (order) {
+      case 0:
+        for (int i = 0; i < vec->n; i++) {                                                   out->data.F64[i] = 1.0; } break;
+      case 1:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i];                       out->data.F64[i] = x; } break;
+      case 2:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = 2.0*x2 - 1.0; } break;
+      case 3:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x*(4.0*x2 - 3.0); } break;
+      case 4:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x2*(8.0*x2 - 8.0) + 1.0; } break;
+      case 5:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x *(x2*(16.0*x2 - 20.0) + 5.0); } break;
+      case 6:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x2*(x2*(32.0*x2 - 48.0) + 18.0) - 1.0; } break;
+      case 7:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x *(x2*(x2*(64.0*x2 - 112.0) + 56.0) - 7.0); } break;
+      case 8:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x2*(x2*(x2*(128.0*x2 - 256.0) + 160.0) - 32.0) + 1.0; } break;
+      case 9:
+        for (int i = 0; i < vec->n; i++) { psF64 x = vec->data.F64[i]; psF64 x2 = PS_SQR(x); out->data.F64[i] = x *(x2*(x2*(x2*(256.0*x2 - 576.0) + 432.0) - 129.0) + 9.0); } break;
+      default:
+        psWarning ("Chebyshev orders higher than 9 are not yet coded\n");
+        psFree (out);
+        return NULL;
+    }
+
+    return out;
+}
+
 /*****************************************************************************
     Polynomial coefficients will be accessed in [w][x][y][z] fashion.
@@ -234 +333 @@
 
 // XXX: You can do this without having to psAlloc() vector d.
-// XXX: How does the mask vector effect Crenshaw's formula?
+// XXX: How does the mask vector affect Clenshaw's formula?
 // NOTE: We assume that x is scaled between -1.0 and 1.0;
 // XXX: Create a faster version for low-order Chebyshevs.
@@ -343 +442 @@
                                   const psPolynomial2D* poly)
 {
-    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
-    PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
+  // XXX transform x,y to chebyshev range
+  // PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, 0.0);
+  // PS_ASSERT_DOUBLE_WITHIN_RANGE(y, -1.0, 1.0, 0.0);
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
 
@@ -354 +454 @@
     unsigned int maxChebyPoly = 0;
 
+    psF64 xNorm = x*poly->scale[0] + poly->zero[0];
+    psF64 yNorm = y*poly->scale[1] + poly->zero[1];
+
     // Determine how many Chebyshev polynomials
     // are needed, then create them.
@@ -366 +469 @@
             if (!(poly->coeffMask[loop_x][loop_y] & PS_POLY_MASK_SET)) {
                 polySum += poly->coeff[loop_x][loop_y] *
-                           psPolynomial1DEval(chebPolys[loop_x], x) *
-                           psPolynomial1DEval(chebPolys[loop_y], y);
+                           psPolynomial1DEval(chebPolys[loop_x], xNorm) *
+                           psPolynomial1DEval(chebPolys[loop_y], yNorm);
             }
         }
@@ -603 +706 @@
     }
 
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    newPoly->scale[0] = NAN;
+    newPoly->zero[0]  = NAN;
+
     return(newPoly);
 }
@@ -638 +747 @@
             newPoly->coeffMask[x][y] = PS_POLY_MASK_NONE;
         }
+    }
+
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    for (int i = 0; i < 2; i++) {
+      newPoly->scale[i] = NAN;
+      newPoly->zero[i]  = NAN;
     }
 
@@ -756 +873 @@
             }
         }
+    }
+
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    for (int i = 0; i < 3; i++) {
+      newPoly->scale[i] = NAN;
+      newPoly->zero[i]  = NAN;
     }
 
@@ -820 +945 @@
     }
 
+    // scale & zero are used for Chebyshev polynomials to define the relationship between
+    // the independent variables and the normalized version with range -1 : +1.  These
+    // must be determined for a specific data set.
+    for (int i = 0; i < 4; i++) {
+      newPoly->scale[i] = NAN;
+      newPoly->zero[i]  = NAN;
+    }
+
     return(newPoly);
 }
@@ -841 +974 @@
 
 // this function must accept F32 and F64 input x vectors
+// EAM XXX these functions seem inefficiently implemented with many nested function calls.
+// they might benefit from unrolling.
 psVector *psPolynomial1DEvalVector(const psPolynomial1D *poly,
                                    const psVector *x)
@@ -888 +1023 @@
 }
 
+psVector *psPolynomial2DEvalChebVector(const psPolynomial2D *poly,
+                                       const psVector *x,
+                                       const psVector *y)
+{
+
+    if (!isfinite(poly->scale[0]) || !isfinite(poly->zero[0]) || !isfinite(poly->scale[1]) || !isfinite(poly->zero[1])) {
+        // re-calculate if not already determined?  
+        psError(PS_ERR_UNKNOWN, true, "normalization scales are not set for chebyshev polynomial");
+        return (NULL);
+    }
+
+    // Number of polynomial terms
+    int nXterm = 1 + poly->nX;      // Number of terms in x
+    int nYterm = 1 + poly->nY;      // Number of terms in y
+    if (nXterm > 9) {
+        psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n");
+        return NULL;
+    }
+    if (nYterm > 9) {
+        psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n");
+        return NULL;
+    }
+
+    // Generate normalized vectors for the range -1 : +1.  These functions cast to psF64
+    psVector *xNorm = psChebyshevNormVector (poly, x, 0);
+    psVector *yNorm = psChebyshevNormVector (poly, y, 1);
+    
+    // Generate the N cheb polynomials based on xNorm, yNorm
+    psArray *xPolySet = psArrayAlloc (nXterm);
+    for (int i = 0; i < nXterm; i++) {
+        xPolySet->data[i] = psChebyshevPolyVector (xNorm, i);
+    }
+    psArray *yPolySet = psArrayAlloc (nYterm);
+    for (int i = 0; i < nYterm; i++) {
+        yPolySet->data[i] = psChebyshevPolyVector (yNorm, i);
+    }
+
+    psVector *out = psVectorAlloc (x->n, PS_TYPE_F64);
+
+    psF64 *xData = xNorm->data.F64;
+    psF64 *yData = yNorm->data.F64;
+    psF64 *fData = out->data.F64;
+
+    // loop over all elements of the data vector
+    for (int i = 0; i < x->n; i++) {
+
+        if (!finite(xData[i])) {fData[i] = NAN; continue; }
+        if (!finite(yData[i])) {fData[i] = NAN; continue; }
+
+        psF64 sum = 0.0;
+        for (int jx = 0; jx < nXterm; jx++) {
+            psVector *jxCheb = xPolySet->data[jx];
+            for (int jy = 0; jy < nYterm; jy++) {
+                psVector *jyCheb = yPolySet->data[jy];
+                sum += poly->coeff[jx][jy] * jxCheb->data.F64[i] * jyCheb->data.F64[i];
+            }
+        }
+        fData[i] = sum;
+    }
+
+    psFree (xPolySet);
+    psFree (yPolySet);
+    psFree (xNorm);
+    psFree (yNorm);
+
+    return out;
+}
+
 // this function must support input data types of F32 and F64
 // all input vectors data types must match (all F32 or all F64)
@@ -901 +1104 @@
     PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL);
 
-    psVector *tmp;
-    unsigned int vecLen=x->n;
+    unsigned int vecLen = x->n;
+
+    // input vector types must match
+    if (y->type.type != x->type.type) {
+        psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors");
+        return (NULL);
+    }
 
     // Determine the length of the output vector to by the minimum of the x,y vectors
-    if (y->n < vecLen) {
-        vecLen = y->n;
+    // XXX shouldn't we require x & y to have the same length?  seems meaningless otherwise
+    if (y->n != vecLen) {
+        psError(PS_ERR_UNKNOWN, true, "length mismatch in data vectors");
+        return (NULL);
+    }
+
+    if (poly->type == PS_POLYNOMIAL_CHEB) {
+        psVector *out = psPolynomial2DEvalChebVector (poly, x, y);
+        return out;
     }
 
     switch (x->type.type) {
-    case PS_TYPE_F32:
-        if (y->type.type != x->type.type) {
-            psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors");
-            return (NULL);
-        }
-
-        // Create output vector to return
-        tmp = psVectorAlloc(vecLen, PS_TYPE_F32);
-
-        // Evaluate the polynomial at the specified points
-        for (unsigned int i=0; i<vecLen; i++) {
-            tmp->data.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]);
-        }
-        break;
-    case PS_TYPE_F64:
-        if (y->type.type != x->type.type) {
-            psError(PS_ERR_UNKNOWN, true, "type mismatch in data vectors");
-            return (NULL);
-        }
-
-        // Create output vector to return
-        tmp = psVectorAlloc(vecLen, PS_TYPE_F64);
-
-        // Evaluate the polynomial at the specified points
-        for (unsigned int i=0; i<vecLen; i++) {
-            tmp->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
-        }
-        break;
-    default:
+      case PS_TYPE_F32: {
+          // Create output vector to return
+          psVector *out = psVectorAlloc(vecLen, PS_TYPE_F32);
+
+          // Evaluate the polynomial at the specified points
+          for (unsigned int i = 0; i < vecLen; i++) {
+              out->data.F32[i] = psPolynomial2DEval(poly,x->data.F32[i],y->data.F32[i]);
+          }
+          return out;
+      }
+      case PS_TYPE_F64: {
+          // Create output vector to return
+          psVector *out = psVectorAlloc(vecLen, PS_TYPE_F64);
+
+          // Evaluate the polynomial at the specified points
+          for (unsigned int i = 0; i < vecLen; i++) {
+              out->data.F64[i] = psPolynomial2DEval(poly,x->data.F64[i],y->data.F64[i]);
+          }
+          return out;
+      }
+      default:
         psError(PS_ERR_UNKNOWN, false, "invalid input data type.\n");
         return (NULL);
     }
-    // Return output vector
-    return(tmp);
+    psAbort ("impossible");
+    return NULL;
 }