Index: branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.c (revision 41831) @@ -40,4 +40,10 @@ #include "psMinimizePolyFit.h" +# define TEST_SAVE_INVERSE_TRANSFORM 0 + +# if (TEST_SAVE_INVERSE_TRANSFORM) +# include "psBinaryOp.h" +# endif + # define ELIXIR_CODE 1 @@ -86,5 +92,6 @@ } - psPlaneTransform *out = psPlaneTransformAlloc(1, 1); + // since the output polynomial is 1st order, a Chebyshev is not really useful + psPlaneTransform *out = psPlaneTransformAlloc(1, 1, PS_POLYNOMIAL_ORD); psF64 r12 = 0.0; @@ -191,5 +198,5 @@ } -psPlaneTransform* psPlaneTransformAlloc(int order1, int order2) +psPlaneTransform* psPlaneTransformAlloc(int order1, int order2, psPolynomialType type) { PS_ASSERT_INT_NONNEGATIVE(order1, NULL); @@ -198,6 +205,6 @@ psPlaneTransform *pt = psAlloc(sizeof(psPlaneTransform)); - pt->x = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order1, order2); - pt->y = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, order1, order2); + pt->x = psPolynomial2DAlloc(type, order1, order2); + pt->y = psPolynomial2DAlloc(type, order1, order2); psMemSetDeallocator(pt, (psFreeFunc) planeTransformFree); @@ -797,4 +804,10 @@ } + // both polynomials in both input transforms must match type -- and for now be ORD + psAssert (trans1->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB"); + psAssert (trans1->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB"); + psAssert (trans2->x->type == PS_POLYNOMIAL_ORD, "fix for CHEB"); + psAssert (trans2->y->type == PS_POLYNOMIAL_ORD, "fix for CHEB"); + // // Determine the size of the new psPlaneTransform. @@ -814,5 +827,5 @@ psPlaneTransform *myPT = NULL; if (out == NULL) { - myPT = psPlaneTransformAlloc(orderX, orderY); + myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD); } else { if ((out->x->nX == orderX) && @@ -834,5 +847,5 @@ } else { psFree(out); - myPT = psPlaneTransformAlloc(orderX, orderY); + myPT = psPlaneTransformAlloc(orderX, orderY, PS_POLYNOMIAL_ORD); } } @@ -1037,5 +1050,5 @@ psPlaneTransform *myPT = NULL; if (out == NULL) { - myPT = psPlaneTransformAlloc(order, order); + myPT = psPlaneTransformAlloc(order, order, PS_POLYNOMIAL_CHEB); } else { // the user has supplied a model with a specific order : fit that order @@ -1071,4 +1084,59 @@ result &= psVectorFitPolynomial2D(myPT->x, NULL, 0, xOut, NULL, xIn, yIn); result &= psVectorFitPolynomial2D(myPT->y, NULL, 0, yOut, NULL, xIn, yIn); + +# if (TEST_SAVE_INVERSE_TRANSFORM) + + psVector *xFit = psPolynomial2DEvalVector (myPT->x, xIn, yIn); + psVector *xRes = (psVector *) psBinaryOp (NULL, xOut, "-", xFit); + + psVector *yFit = psPolynomial2DEvalVector (myPT->y, xIn, yIn); + psVector *yRes = (psVector *) psBinaryOp (NULL, yOut, "-", yFit); + + static int nOut = 0; + char filename[1024]; + snprintf (filename, 1024, "test.fit.%03d.ply", nOut); + FILE *fp = fopen (filename, "w"); + + fprintf (fp, " ---- xFit ---- \n"); + + for (int ix = 0; ix < myPT->x->nX + 1; ix++) { + for (int iy = 0; iy < myPT->x->nY + 1; iy++) { + fprintf (fp, "%18.12e ", myPT->x->coeff[ix][iy]); + } + fprintf (fp, "\n"); + } + fprintf (fp, " ---- yFit ---- \n"); + + for (int ix = 0; ix < myPT->y->nX + 1; ix++) { + for (int iy = 0; iy < myPT->y->nY + 1; iy++) { + fprintf (fp, "%18.12e ", myPT->y->coeff[ix][iy]); + } + fprintf (fp, "\n"); + } + fclose (fp); + + snprintf (filename, 1024, "test.fit.%03d.dat", nOut); nOut ++; + FILE *f1 = fopen (filename, "w"); + for (int i = 0; i < xFit->n; i++) { + fprintf (f1, "%d : %f %f : %f %f : %f %f : %f %f\n", i, + xIn->data.F64[i], yIn->data.F64[i], + xOut->data.F64[i], yOut->data.F64[i], + xFit->data.F64[i], yFit->data.F64[i], + xRes->data.F64[i], yRes->data.F64[i]); + } + fclose (f1); + + psStats *myStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV); + psVectorStats (myStats, xRes, NULL, NULL, 0); float dX = myStats->sampleStdev; + psVectorStats (myStats, yRes, NULL, NULL, 0); float dY = myStats->sampleStdev; + fprintf (stderr, "xRes Sigma: %f -- yRes Sigma %f\n", dX, dY); + + psFree (myStats); + psFree (xFit); + psFree (yFit); + psFree (xRes); + psFree (yRes); + +# endif psFree(inCoord); Index: branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.h =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.h (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/psLib/src/astro/psCoord.h (revision 41831) @@ -178,6 +178,7 @@ psPlaneTransform* psPlaneTransformAlloc( - int order1, ///< The order of the x term in the transform. - int order2 ///< The order of the y term in the transform. + int order1, ///< The order of the x term in the transform. + int order2, ///< The order of the y term in the transform. + psPolynomialType type ///< The polynomial type (ORD or CHEB) for this transform ) PS_ATTR_MALLOC; Index: branches/eam_branches/ipp-dev-20210817/psLib/src/db/psDB.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/src/db/psDB.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/psLib/src/db/psDB.c (revision 41831) @@ -2875,37 +2875,36 @@ switch (pType) { - case PS_DATA_S8: + case PS_DATA_S8: isNaN = PS_IS_NAN(psS8, data, PS_MAX_S8); break; - case PS_DATA_S16: + case PS_DATA_S16: isNaN = PS_IS_NAN(psS16, data, PS_MAX_S16); break; - case PS_DATA_S32: + case PS_DATA_S32: isNaN = PS_IS_NAN(psS32, data, PS_MAX_S32); break; - case PS_DATA_S64: + case PS_DATA_S64: isNaN = PS_IS_NAN(psS64, data, PS_MAX_S64); break; - case PS_DATA_U8: + case PS_DATA_U8: isNaN = PS_IS_NAN(psU8, data, PS_MAX_U8); break; - case PS_DATA_U16: + case PS_DATA_U16: isNaN = PS_IS_NAN(psU16, data, PS_MAX_U16); break; - case PS_DATA_U32: + case PS_DATA_U32: isNaN = PS_IS_NAN(psU32, data, PS_MAX_U32); break; - case PS_DATA_U64: + case PS_DATA_U64: isNaN = PS_IS_NAN(psU64, data, PS_MAX_U64); break; - case PS_DATA_F32: - isNaN = isnan(*((psF32 *) data)); - break; - case PS_DATA_F64: - isNaN = isnan(*((psF64 *) data)); - break; - case PS_DATA_BOOL: - // XXX: what is NaN for a bool? - isNaN = PS_IS_NAN(psU8, data, PS_MAX_U8); + case PS_DATA_F32: + isNaN = !isfinite(*((psF32 *) data)); // trap nan, +inf, -inf + break; + case PS_DATA_F64: + isNaN = !isfinite(*((psF64 *) data)); // trap nan, +inf, -inf + break; + case PS_DATA_BOOL: + isNaN = PS_IS_NAN(psU8, data, PS_MAX_U8); // probably meaningless break; } Index: branches/eam_branches/ipp-dev-20210817/psLib/src/math/psMinimizePolyFit.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/src/math/psMinimizePolyFit.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/psLib/src/math/psMinimizePolyFit.c (revision 41831) @@ -35,4 +35,5 @@ #include "psMinimizePolyFit.h" +#include "psAbort.h" #include "psAssert.h" #include "psMinimizeLMM.h" // For Gauss-Jordan routines @@ -1225,4 +1226,153 @@ /****************************************************************************** +VectorFitPolynomial2DCheb(myPoly, *mask, maskValue, *f, *fErr, *x, *y): This is +a private routine which will fit a 2-D polynomial to a set of (x, y)-(f) +pairs. All non-NULL vectors must be of type PS_TYPE_F64. + + *****************************************************************************/ +static bool VectorFitPolynomial2DCheb( + psPolynomial2D* myPoly, + const psVector *f, + const psVector *x, + const psVector *y) +{ + psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + + // Number of polynomial terms + int nXterm = 1 + myPoly->nX; // Number of terms in x + int nYterm = 1 + myPoly->nY; // Number of terms in y + int nTerm = nXterm * nYterm; // Total number of terms + if (nXterm > 9) { + psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n"); + return false; + } + if (nYterm > 9) { + psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit: orders higher than 9 are not yet coded\n"); + return false; + } + + // determine scale factors + if (!psChebyshevSetScale (myPoly, x, 0)) { psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit.\n"); return false; } + if (!psChebyshevSetScale (myPoly, y, 1)) { psError(PS_ERR_UNKNOWN, false, "failed 2D chebyshev fit.\n"); return false; } + + // generate normalized vectors + psVector *xNorm = psChebyshevNormVector (myPoly, x, 0); + psVector *yNorm = psChebyshevNormVector (myPoly, 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); + } + + psF64 *fData = f->data.F64; // Dereference f + + psImage *A = psImageAlloc(nTerm, nTerm, PS_TYPE_F64); // Least-squares matrix + psVector *B = psVectorAlloc(nTerm, PS_TYPE_F64); // Least-squares vector + + // Initialize data structures (should not be able to fail) + psAssert (psImageInit(A, 0.0), "Could initialize data structures A"); + psAssert (psVectorInit(B, 0.0), "Could initialize data structures B"); + + // Dereference stuff, to make the loop go faster + psF64 **matrix = A->data.F64; // Dereference the least-squares matrix + psF64 *vector = B->data.F64; // Dereference the least-squares vector + + // loop over all elements of the data vector + for (int k = 0; k < x->n; k++) { + + if (!finite(fData[k])) continue; + + // XXX can we only calculate the upper diagonal? + int nelem = 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]; + psF64 chebValue = jxCheb->data.F64[k] * jyCheb->data.F64[k]; + + vector[nelem] += fData[k] * chebValue; + + int melem = 0; + for (int kx = 0; kx < nXterm; kx++) { + psVector *kxCheb = xPolySet->data[kx]; + for (int ky = 0; ky < nYterm; ky++) { + psVector *kyCheb = yPolySet->data[ky]; + matrix[nelem][melem] += chebValue * kxCheb->data.F64[k]*kyCheb->data.F64[k]; + melem++; + } + } + nelem++; + } + } + } + + if (psTraceGetLevel("psLib.math") >= 4) { + printf("Least-squares vector:\n"); + for (int i = 0; i < nTerm; i++) { + printf("%f ", B->data.F64[i]); + } + printf("\n"); + printf("Least-squares matrix:\n"); + for (int i = 0; i < nTerm; i++) { + for (int j = 0; j < nTerm; j++) { + printf("%f ", A->data.F64[i][j]); + } + printf("\n"); + } + } + + bool status = false; + if (USE_GAUSS_JORDAN) { + status = psMatrixGJSolve(A, B); + } else { + status = psMatrixLUSolve(A, B); + } + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.\n"); + goto escape; + } + + // unroll the result: + int nelem = 0; + for (int jx = 0; jx < nXterm; jx++) { + for (int jy = 0; jy < nYterm; jy++) { + myPoly->coeff[jx][jy] = B->data.F64[nelem]; + myPoly->coeffErr[jx][jy] = sqrt(A->data.F64[nelem][nelem]); + nelem ++; + } + } + psFree(A); + psFree(B); + + psFree (xNorm); + psFree (yNorm); + psFree (xPolySet); + psFree (yPolySet); + + return true; + +escape: + psFree (A); + psFree (B); + return false; +} + +/****************************************************************************** psVectorFitPolynomial2D(): This routine fits a 2D polynomial of arbitrary degree (specified in poly) to the data points (x, y)-(f) and returns that @@ -1240,5 +1390,5 @@ { PS_ASSERT_POLY_NON_NULL(poly, false); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + // PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); PS_ASSERT_VECTOR_NON_NULL(f, false); @@ -1277,10 +1427,15 @@ break; case PS_POLYNOMIAL_CHEB: - if (mask != NULL) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); - } - psError(PS_ERR_UNKNOWN, true, "2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - result = false; - break; + if (mask != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + } + if (fErr != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring error values for Chebyshev polynomials.\n"); + } + result = VectorFitPolynomial2DCheb(poly, f64, x64, y64); + if (!result) { + psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); + } + break; default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); Index: branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.c (revision 41831) @@ -36,4 +36,5 @@ #include "psLogMsg.h" #include "psPolynomial.h" +#include "psAbort.h" #include "psAssert.h" @@ -203,4 +204,102 @@ +/** 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,5 +333,5 @@ // 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,6 +442,7 @@ 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,4 +454,7 @@ 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,6 +469,6 @@ 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,4 +706,10 @@ } + // 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,4 +747,12 @@ 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,4 +873,12 @@ } } + } + + // 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,4 +945,12 @@ } + // 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,4 +974,6 @@ // 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,4 +1023,72 @@ } +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,47 +1104,51 @@ 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; idata.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; idata.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; } Index: branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.h =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.h (revision 41830) +++ branches/eam_branches/ipp-dev-20210817/psLib/src/math/psPolynomial.h (revision 41831) @@ -68,4 +68,6 @@ psF64 *coeffErr; ///< Error in coefficients psMaskType *coeffMask; ///< Coefficient mask + double scale[1]; ///< Chebyshev scale factor + double zero[1]; ///< Chebyshev zero point } psPolynomial1D; @@ -80,4 +82,6 @@ psF64 **coeffErr; ///< Error in coefficients psMaskType **coeffMask; ///< Coefficients mask + double scale[2]; ///< Chebyshev scale factor + double zero[2]; ///< Chebyshev zero point } psPolynomial2D; @@ -93,4 +97,6 @@ psF64 ***coeffErr; ///< Error in coefficients psMaskType ***coeffMask; ///< Coefficients mask + double scale[3]; ///< Chebyshev scale factor + double zero[3]; ///< Chebyshev zero point } psPolynomial3D; @@ -107,4 +113,6 @@ psF64 ****coeffErr; ///< Error in coefficients psMaskType ****coeffMask; ///< Coefficients mask + double scale[4]; ///< Chebyshev scale factor + double zero[4]; ///< Chebyshev zero point } psPolynomial4D; @@ -305,4 +313,9 @@ p_chebyPolys; + +// chebyshev support functions: +bool psChebyshevSetScale (psPolynomial2D* myPoly, const psVector *vec, int dir); +psVector *psChebyshevNormVector (const psPolynomial2D* myPoly, const psVector *vec, int dir); +psVector *psChebyshevPolyVector (const psVector *vec, int order); /***************************************************************************** Index: branches/eam_branches/ipp-dev-20210817/psLib/test/math/tap_psPolyFit2DCheb.c =================================================================== --- branches/eam_branches/ipp-dev-20210817/psLib/test/math/tap_psPolyFit2DCheb.c (revision 41831) +++ branches/eam_branches/ipp-dev-20210817/psLib/test/math/tap_psPolyFit2DCheb.c (revision 41831) @@ -0,0 +1,84 @@ +#include +#include +#include +//#include "tap.h" +//#include "pstap.h" + +int main(int argc, char **argv) +{ + psLogSetFormat("HLNM"); + psLogSetLevel(PS_LOG_INFO); + // plan_tests(88); + + if (argc != 4) { + fprintf (stderr, "USAGE: tap_psPolyFit2DCheb (input) (Nx) (Ny)\n"); + exit (2); + } + + FILE *f = fopen (argv[1], "r"); + int Nx = atoi(argv[2]); + int Ny = atoi(argv[3]); + + // data file should contain x y p + // determine the best fit 2D polynomial from x y to p + psVector *x = psVectorAllocEmpty(1000, PS_TYPE_F64); + psVector *y = psVectorAllocEmpty(1000, PS_TYPE_F64); + psVector *p = psVectorAllocEmpty(1000, PS_TYPE_F64); + + double xv, yv, pv; + + while (fscanf (f, "%lf %lf %lf", &xv, &yv, &pv) == 3) { + psVectorAppend (x, xv); + psVectorAppend (y, yv); + psVectorAppend (p, pv); + } + + // Fit an ordinary polynomial + { + psPolynomial2D *polyOrd = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, Nx, Ny); + int result = psVectorFitPolynomial2D(polyOrd, NULL, 0, p, NULL, x, y); + + fprintf (stdout, " ---- Ord ---- \n"); + + for (int ix = 0; ix < polyOrd->nX + 1; ix++) { + for (int iy = 0; iy < polyOrd->nY + 1; iy++) { + fprintf (stdout, "%18.12e ", polyOrd->coeff[ix][iy]); + } + fprintf (stdout, "\n"); + } + fprintf (stdout, "\n"); + + psVector *P = psPolynomial2DEvalVector (polyOrd, x, y); + + FILE *g1 = fopen ("ord.dat", "w"); + for (int i = 0; i < x->n; i++) { + fprintf (g1, "%lf %lf %lf %lf\n", x->data.F64[i], y->data.F64[i], p->data.F64[i], P->data.F64[i]); + } + fclose (g1); + } + + // Fit a chebychev polynomial + { + psPolynomial2D *polyCheb = psPolynomial2DAlloc(PS_POLYNOMIAL_CHEB, Nx, Ny); + + int result = psVectorFitPolynomial2D(polyCheb, NULL, 0, p, NULL, x, y); + + fprintf (stdout, " ---- Cheb ---- \n"); + + for (int ix = 0; ix < polyCheb->nX + 1; ix++) { + for (int iy = 0; iy < polyCheb->nY + 1; iy++) { + fprintf (stdout, "%18.12e ", polyCheb->coeff[ix][iy]); + } + fprintf (stdout, "\n"); + } + + psVector *P = psPolynomial2DEvalVector (polyCheb, x, y); + + FILE *g1 = fopen ("cheb.dat", "w"); + for (int i = 0; i < x->n; i++) { + fprintf (g1, "%lf %lf %lf %lf\n", x->data.F64[i], y->data.F64[i], p->data.F64[i], P->data.F64[i]); + } + fclose (g1); + } +} +/// gcc -o tap_psPolyFit2DCheb tap_psPolyFit2DCheb.c -lm -lpslib -I /home/real/eugene/src/psconfig/ipp-trunk-20210518.lin64/include/pslib -L /home/real/eugene/src/psconfig/ipp-trunk-20210518.lin64/lib