Changeset 41896 for trunk/psLib/src/math/psPolynomial.c

Timestamp:

Nov 4, 2021, 6:10:51 PM (5 years ago)

Author:

eugene

Message:

merge changes from eam_branches/ipp-dev-20210817 (fix chebyshevs 1D, 2D, set DB field to NULL for inf, psFitsTableNew)

File:

• trunk/psLib/src/math/psPolynomial.c (modified) (12 diffs)

trunk/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;
+}
+
+# define CHEB_EVAL_0(OUT,IN) {OUT = 1.0;}
+# define CHEB_EVAL_1(OUT,IN) {                       OUT = IN; }
+# define CHEB_EVAL_2(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = 2.0*X2 - 1.0; }
+# define CHEB_EVAL_3(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN*(4.0*X2 - 3.0); }
+# define CHEB_EVAL_4(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(8.0*X2 - 8.0) + 1.0; }
+# define CHEB_EVAL_5(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(16.0*X2 - 20.0) + 5.0); }
+# define CHEB_EVAL_6(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(X2*(32.0*X2 - 48.0) + 18.0) - 1.0; }
+# define CHEB_EVAL_7(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(X2*(64.0*X2 - 112.0) + 56.0) - 7.0); }
+# define CHEB_EVAL_8(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = X2*(X2*(X2*(128.0*X2 - 256.0) + 160.0) - 32.0) + 1.0; }
+# define CHEB_EVAL_9(OUT,IN) {psF64 X2 = PS_SQR(IN); OUT = IN *(X2*(X2*(X2*(256.0*X2 - 576.0) + 432.0) - 129.0) + 9.0); }
+
+/** 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++) { CHEB_EVAL_0(out->data.F64[i], vec->data.F64[i]); } break;
+      case 1:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_1(out->data.F64[i], vec->data.F64[i]); } break;
+      case 2:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_2(out->data.F64[i], vec->data.F64[i]); } break;
+      case 3:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_3(out->data.F64[i], vec->data.F64[i]); } break;
+      case 4:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_4(out->data.F64[i], vec->data.F64[i]); } break;
+      case 5:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_5(out->data.F64[i], vec->data.F64[i]); } break;
+      case 6:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_6(out->data.F64[i], vec->data.F64[i]); } break;
+      case 7:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_7(out->data.F64[i], vec->data.F64[i]); } break;
+      case 8:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_8(out->data.F64[i], vec->data.F64[i]); } break;
+      case 9:
+        for (int i = 0; i < vec->n; i++) { CHEB_EVAL_9(out->data.F64[i], vec->data.F64[i]); } 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.
@@ -233 +343 @@
 }
 
-// XXX: You can do this without having to psAlloc() vector d.
-// XXX: How does the mask vector effect Crenshaw's formula?
-// NOTE: We assume that x is scaled between -1.0 and 1.0;
-// XXX: Create a faster version for low-order Chebyshevs.
-static psF64 chebPolynomial1DEval(
-    psF64 x,
-    const psPolynomial1D* poly)
-{
-    PS_ASSERT_DOUBLE_WITHIN_RANGE(x, -1.0, 1.0, NAN);
+static psF64 chebPolynomial1DEval(psF64 x, const psPolynomial1D* poly) {
+
     PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(poly->nX, 0, NAN);
+
+    psF64 xNorm = x*poly->scale[0] + poly->zero[0];
+
+    psF64 polySum = 0.0;
+
+    for (int ix = 0; ix <= poly->nX; ix++) {
+        if (poly->coeffMask[ix] & PS_POLY_MASK_SET) continue;
+        psF64 xCheb = NAN;
+        switch (ix) {
+          case 0: CHEB_EVAL_0 (xCheb, xNorm); break;
+          case 1: CHEB_EVAL_1 (xCheb, xNorm); break;
+          case 2: CHEB_EVAL_2 (xCheb, xNorm); break;
+          case 3: CHEB_EVAL_3 (xCheb, xNorm); break;
+          case 4: CHEB_EVAL_4 (xCheb, xNorm); break;
+          case 5: CHEB_EVAL_5 (xCheb, xNorm); break;
+          case 6: CHEB_EVAL_6 (xCheb, xNorm); break;
+          case 7: CHEB_EVAL_7 (xCheb, xNorm); break;
+          case 8: CHEB_EVAL_8 (xCheb, xNorm); break;
+          case 9: CHEB_EVAL_9 (xCheb, xNorm); break;
+          default:
+            break;
+        }
+        polySum += poly->coeff[ix] * xCheb;
+    }
+    return polySum;
+}
+
+/*** version 1 is a general case and could be used for Norder > 9.  ***/
+# ifdef CHEB_VERSION_1
+void oldcode_1(void) {
     psVector *d;
+    psF64 tmp = 0.0;
 
     unsigned int nTerms = 1 + poly->nX;
     unsigned int i;
-    psF64 tmp = 0.0;
 
     // Special case where the Chebyshev poly is constant.
@@ -268 +401 @@
     }
 
-    if (1) {
-        // General case where the Chebyshev poly has 2 or more terms.
-        d = psVectorAlloc(nTerms, PS_TYPE_F64);
-        if (!(poly->coeffMask[nTerms-1] & PS_POLY_MASK_SET)) {
-            d->data.F64[nTerms-1] = poly->coeff[nTerms-1];
-        } else {
-            d->data.F64[nTerms-1] = 0.0;
-        }
-
-        d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]);
-        if (!(poly->coeffMask[nTerms-2] & PS_POLY_MASK_SET)) {
-            d->data.F64[nTerms-2] += poly->coeff[nTerms-2];
-        }
-
-        for (i=nTerms-3;i>=1;i--) {
-            d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]);
-            if (!(poly->coeffMask[i] & PS_POLY_MASK_SET)) {
-                d->data.F64[i] += poly->coeff[i];
-            }
-        }
-
-        tmp = (x * d->data.F64[1]) - (d->data.F64[2]);
-        if (!(poly->coeffMask[0] & PS_POLY_MASK_SET)) {
-            tmp += (0.5 * poly->coeff[0]);
-        }
-        psFree(d);
+    // General case where the Chebyshev poly has 2 or more terms.
+    d = psVectorAlloc(nTerms, PS_TYPE_F64);
+    if (!(poly->coeffMask[nTerms-1] & PS_POLY_MASK_SET)) {
+        d->data.F64[nTerms-1] = poly->coeff[nTerms-1];
     } else {
-        // XXX: This is old code that does not use Clenshaw's formula.  Get rid of it.
-        psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX);
-
-        tmp = 0.0;
-        for (psS32 i=0;i<(1 + poly->nX);i++) {
-            tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x));
-        }
-        tmp-= (poly->coeff[0]/2.0);
-
-        for (psS32 i=0;i<(1 + poly->nX);i++) {
-            psFree(chebPolys[i]);
-        }
-        psFree(chebPolys);
-    }
+        d->data.F64[nTerms-1] = 0.0;
+    }
+
+    d->data.F64[nTerms-2] = (2.0 * x * d->data.F64[nTerms-1]);
+    if (!(poly->coeffMask[nTerms-2] & PS_POLY_MASK_SET)) {
+        d->data.F64[nTerms-2] += poly->coeff[nTerms-2];
+    }
+
+    for (i=nTerms-3;i>=1;i--) {
+        d->data.F64[i] = (2.0 * x * d->data.F64[i+1]) - (d->data.F64[i+2]);
+        if (!(poly->coeffMask[i] & PS_POLY_MASK_SET)) {
+            d->data.F64[i] += poly->coeff[i];
+        }
+    }
+
+    tmp = (x * d->data.F64[1]) - (d->data.F64[2]);
+    if (!(poly->coeffMask[0] & PS_POLY_MASK_SET)) {
+        tmp += (0.5 * poly->coeff[0]);
+    }
+    psFree(d);
+}
+# endif
+
+/*** version 0 should be removed when version 2 is ready ***/
+# ifdef CHEB_VERSION_0
+void oldcode_0(void) {
+    // XXX: This is old code that does not use Clenshaw's formula.  Get rid of it.
+    psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX);
+
+    tmp = 0.0;
+    for (psS32 i=0;i<(1 + poly->nX);i++) {
+        tmp+= (poly->coeff[i] * psPolynomial1DEval(chebPolys[i], x));
+    }
+    tmp-= (poly->coeff[0]/2.0);
+
+    for (psS32 i=0;i<(1 + poly->nX);i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
 
     return(tmp);
 }
+# endif
 
 static psF64 ordPolynomial2DEval(psF64 x,
@@ -343 +480 @@
                                   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);
     PS_ASSERT_POLY_NON_NULL(poly, NAN);
 
-    unsigned int loop_x = 0;
-    unsigned int loop_y = 0;
-    unsigned int i = 0;
+    psF64 xNorm = x*poly->scale[0] + poly->zero[0];
+    psF64 yNorm = y*poly->scale[1] + poly->zero[1];
+
     psF64 polySum = 0.0;
-    psPolynomial1D* *chebPolys = NULL;
-    unsigned int maxChebyPoly = 0;
-
-    // Determine how many Chebyshev polynomials
-    // are needed, then create them.
-    maxChebyPoly = poly->nX;
-    if (poly->nY > maxChebyPoly) {
-        maxChebyPoly = poly->nY;
-    }
-    chebPolys = p_psCreateChebyshevPolys(maxChebyPoly + 1);
-
-    for (loop_x = 0; loop_x < (1 + poly->nX); loop_x++) {
-        for (loop_y = 0; loop_y < (1 + poly->nY); loop_y++) {
-            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);
-            }
-        }
-    }
-    for (i=0;i<maxChebyPoly+1;i++) {
-        psFree(chebPolys[i]);
-    }
-    psFree(chebPolys);
+
+    // XXX this could be quicker if we saved the N xvalues are re-used the resuls
+    for (int ix = 0; ix <= poly->nX; ix++) {
+        psF64 xCheb = NAN;
+        switch (ix) {
+          case 0: CHEB_EVAL_0 (xCheb, xNorm); break;
+          case 1: CHEB_EVAL_1 (xCheb, xNorm); break;
+          case 2: CHEB_EVAL_2 (xCheb, xNorm); break;
+          case 3: CHEB_EVAL_3 (xCheb, xNorm); break;
+          case 4: CHEB_EVAL_4 (xCheb, xNorm); break;
+          case 5: CHEB_EVAL_5 (xCheb, xNorm); break;
+          case 6: CHEB_EVAL_6 (xCheb, xNorm); break;
+          case 7: CHEB_EVAL_7 (xCheb, xNorm); break;
+          case 8: CHEB_EVAL_8 (xCheb, xNorm); break;
+          case 9: CHEB_EVAL_9 (xCheb, xNorm); break;
+          default:
+            break;
+        }
+        for (int iy = 0; iy <= poly->nY; iy++) {
+            if (poly->coeffMask[ix][iy] & PS_POLY_MASK_SET) continue;
+            psF64 yCheb = NAN;
+            switch (iy) {
+              case 0: CHEB_EVAL_0 (yCheb, yNorm); break;
+              case 1: CHEB_EVAL_1 (yCheb, yNorm); break;
+              case 2: CHEB_EVAL_2 (yCheb, yNorm); break;
+              case 3: CHEB_EVAL_3 (yCheb, yNorm); break;
+              case 4: CHEB_EVAL_4 (yCheb, yNorm); break;
+              case 5: CHEB_EVAL_5 (yCheb, yNorm); break;
+              case 6: CHEB_EVAL_6 (yCheb, yNorm); break;
+              case 7: CHEB_EVAL_7 (yCheb, yNorm); break;
+              case 8: CHEB_EVAL_8 (yCheb, yNorm); break;
+              case 9: CHEB_EVAL_9 (yCheb, yNorm); break;
+              default:
+                break;
+            }
+            polySum += poly->coeff[ix][iy] * xCheb * yCheb;
+        }
+    }
     return(polySum);
 }
@@ -603 +752 @@
     }
 
+    // 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 +793 @@
             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 +919 @@
             }
         }
+    }
+
+    // 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 +991 @@
     }
 
+    // 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);
 }
 
+/* note these functions accept unscaled values and apply the scaling saved on poly */
 psF64 psPolynomial1DEval(const psPolynomial1D* poly,
                          psF64 x)
@@ -830 +1010 @@
     if (poly->type == PS_POLYNOMIAL_ORD) {
         return(ordPolynomial1DEval(x, poly));
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+    }
+    if (poly->type == PS_POLYNOMIAL_CHEB) {
         return(chebPolynomial1DEval(x, poly));
-    } else {
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
-                _("Unknown polynomial type 0x%x found.  Evaluation failed."),
-                poly->type);
-    }
+    } 
+    psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+            _("Unknown polynomial type 0x%x found.  Evaluation failed."),
+            poly->type);
+
     return(NAN);
 }
 
 // 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)
@@ -878 +1061 @@
     if (poly->type == PS_POLYNOMIAL_ORD) {
         return(ordPolynomial2DEval(x, y, poly));
-    } else if (poly->type == PS_POLYNOMIAL_CHEB) {
+    }
+    if (poly->type == PS_POLYNOMIAL_CHEB) {
         return(chebPolynomial2DEval(x, y, poly));
-    } else {
-        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
-                _("Unknown polynomial type 0x%x found.  Evaluation failed."),
-                poly->type);
-    }
+    } 
+    psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+            _("Unknown polynomial type 0x%x found.  Evaluation failed."),
+            poly->type);
     return(NAN);
+}
+
+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];
+                if (poly->coeffMask[jx][jy] & PS_POLY_MASK_SET) continue;
+                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;
 }
 
@@ -901 +1153 @@
     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;
 }