Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 4990) +++ /trunk/psLib/src/astro/psCoord.c (revision 4991) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.85 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-07 21:33:48 $ +* @version $Revision: 1.86 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -545,7 +545,4 @@ int nSamples) { - - // XXX: This does not yet use region and nSamples: need to modify -rdd - PS_ASSERT_PTR_NON_NULL(trans1, NULL); PS_ASSERT_PTR_NON_NULL(trans2, NULL); @@ -905,63 +902,61 @@ ) { - /* - PS_ASSERT_PTR_NON_NULL(transformation, NULL); - PS_ASSERT_POLY_NON_NULL(transformation->x, NULL); - PS_ASSERT_POLY_NON_NULL(transformation->y, NULL); - PS_ASSERT_PTR_NON_NULL(coord, NULL); - - if (out == NULL) { - out = psPlaneAlloc(); - } - - out->x = 0.0; - out->y = 0.0; - out->xErr = 0.0; - out->yErr = 0.0; - - psPolynomial2D *xPoly = transformation->x; - psPolynomial2D *yPoly = transformation->y; - - // - // Calculate the derivative with respect to x. - // - psF32 xSum = 1.0; - psF32 ySum = 1.0; - - // This loop starts at loop_x=1 since the derivative of the loop_x=0 terms are all 0.0 - for (psS32 loop_x = 1; loop_x < xPoly->nX; loop_x++) { - ySum = 1.0; - for (psS32 loop_y = 0; loop_y < xPoly->nY; loop_y++) { - // - // For each iteration of the loop, we multiple the (x, y) coefficient - // by (coord->x^(loop_x-1) * coord->y^loop_y) - // - - out->x+= xPoly->coeff[loop_x][loop_y] * xSum * ySum; - ySum*= coord->y; - } - xSum*= coord->x; - } - - // - // Calculate the derivative with respect to x. - // - xSum = 1.0; - - // This loop starts at loop_y=1 since the derivative of the loop_y=0 terms are all 0.0 - for (psS32 loop_x = 0; loop_x < yPoly->nX; loop_x++) { - ySum = 1.0; - for (psS32 loop_y = 1; loop_y < yPoly->nY; loop_y++) { - // - // For each iteration of the loop, we multiple the (x, y) coefficient - // by (coord->x^(loop_x-1) * coord->y^loop_y) - // - - out->y+= yPoly->coeff[loop_x][loop_y] * xSum * ySum; - ySum*= coord->y; - } - xSum*= coord->x; - } - */ + PS_ASSERT_PTR_NON_NULL(transformation, NULL); + PS_ASSERT_POLY_NON_NULL(transformation->x, NULL); + PS_ASSERT_POLY_NON_NULL(transformation->y, NULL); + PS_ASSERT_PTR_NON_NULL(coord, NULL); + + if (out == NULL) { + out = psPlaneAlloc(); + } + + out->x = 0.0; + out->y = 0.0; + out->xErr = 0.0; + out->yErr = 0.0; + + psPolynomial2D *xPoly = transformation->x; + psPolynomial2D *yPoly = transformation->y; + + // + // Calculate the derivative with respect to x. + // + psF32 xSum = 1.0; + psF32 ySum = 1.0; + + // This loop starts at loop_x=1 since the derivative of the loop_x=0 terms are all 0.0 + for (psS32 loop_x = 1; loop_x < xPoly->nX; loop_x++) { + ySum = 1.0; + for (psS32 loop_y = 0; loop_y < xPoly->nY; loop_y++) { + // + // For each iteration of the loop, we multiple the (x, y) coefficient + // by (coord->x^(loop_x-1) * coord->y^loop_y) + // + + out->x+= xPoly->coeff[loop_x][loop_y] * xSum * ySum; + ySum*= coord->y; + } + xSum*= coord->x; + } + + // + // Calculate the derivative with respect to x. + // + xSum = 1.0; + + // This loop starts at loop_y=1 since the derivative of the loop_y=0 terms are all 0.0 + for (psS32 loop_x = 0; loop_x < yPoly->nX; loop_x++) { + ySum = 1.0; + for (psS32 loop_y = 1; loop_y < yPoly->nY; loop_y++) { + // + // For each iteration of the loop, we multiple the (x, y) coefficient + // by (coord->x^(loop_x-1) * coord->y^loop_y) + // + + out->y+= yPoly->coeff[loop_x][loop_y] * xSum * ySum; + ySum*= coord->y; + } + xSum*= coord->x; + } return(out); } @@ -970,5 +965,5 @@ { if(sphere == NULL) { - // psError(); + psError( PS_ERR_UNKNOWN, true, "psSphere argument is NULL. Returning NULL.\n"); return NULL; } @@ -990,5 +985,5 @@ { if(cube == NULL) { - // psError(); + psError( PS_ERR_UNKNOWN, true, "psCube argument is NULL. Returning NULL.\n"); return NULL; } Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 4990) +++ /trunk/psLib/src/math/psConstants.h (revision 4991) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.75 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-07-20 01:21:13 $ + * @version $Revision: 1.76 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -357,15 +357,22 @@ #define PS_VECTOR_PRINT_F32(NAME) \ -for (int my_i=0;my_i<(NAME)->n;my_i++) { \ - printf("%s->data.F32[%d] is %f\n", #NAME, my_i, (NAME)->data.F32[my_i]); \ -} \ -printf("\n"); \ +if (NAME != NULL) { \ + for (int my_i=0;my_i<(NAME)->n;my_i++) { \ + printf("%s->data.F32[%d] is %f\n", #NAME, my_i, (NAME)->data.F32[my_i]); \ + } \ + printf("\n"); \ +} else {\ + printf("MACRO WARNING: vector %s is NULL.\n", #NAME); \ +}\ #define PS_VECTOR_PRINT_F64(NAME) \ -for (int my_i=0;my_i<(NAME)->n;my_i++) { \ - printf("%s->data.F64[%d] is %f\n", #NAME, my_i, (NAME)->data.F64[my_i]); \ -} \ -printf("\n"); \ - +if (NAME != NULL) { \ + for (int my_i=0;my_i<(NAME)->n;my_i++) { \ + printf("%s->data.F64[%d] is %f\n", #NAME, my_i, (NAME)->data.F64[my_i]); \ + } \ + printf("\n"); \ +} else {\ + printf("MACRO WARNING: vector %s is NULL.\n", #NAME); \ +}\ #define PS_VECTOR_CONVERT_F64_TO_F32_STATIC(OLD, NEW_PTR32, NEW_STATIC32) \ @@ -742,2 +749,4 @@ #define PS_SQR(A) \ ((A) * (A)) + +# define PS_SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;} Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 4990) +++ /trunk/psLib/src/math/psMinimize.c (revision 4991) @@ -8,7 +8,8 @@ * * @author GLG, MHPCC + * @author EAM, IfA * - * @version $Revision: 1.134 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-07 21:35:50 $ + * @version $Revision: 1.135 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -32,4 +33,5 @@ /*****************************************************************************/ /* DEFINE STATEMENTS */ +/* What are the following macros for? */ /*****************************************************************************/ #define PS_SEG psLib @@ -57,5 +59,11 @@ /*****************************************************************************/ +/****************************************************************************** + ****************************************************************************** + Levenberg-Marquadt routines. + ****************************************************************************** + *****************************************************************************/ // measure the distance to the minimum assuming a linear model +// XXX: Move this to the pub area. bool psMinimizeGaussNewtonDelta (psVector *delta, const psVector *params, @@ -112,597 +120,4 @@ } -/****************************************************************************** -buildSums1D(x, polyOrder, sums): this routine calculates the powers of input -parameter "x" between 0 and input parameter polyOrder. The result is returned -as a psVector sums. - -XXX: Use a static vector. - -XXX: should the argument be polyOrder or numTerms? - *****************************************************************************/ -static void buildSums1D(psF64 x, - psS32 polyOrder, - psVector* sums) -{ - psS32 i = 0; - psF64 xSum = 0.0; - psS32 numTerms = polyOrder + 1; - - if (sums == NULL) { - sums = psVectorAlloc(numTerms, PS_TYPE_F64); - } else if (numTerms > sums->n) { - sums = psVectorRealloc(sums, numTerms); - } - - xSum = 1.0; - for (i = 0; i < numTerms; i++) { - sums->data.F64[i] = xSum; - xSum *= x; - } -} - -/***************************************************************************** -CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a -tabulated function at n points, this routine calculates the second -derivatives of the interpolating cubic splines at those n points. - -The first and second derivatives at the endpoints, undefined in the SDR, are -here defined to be 0.0. They can be modified via ypo and yp1. - -This routine assumes that vectors x and y are of the appropriate types/sizes -(F32). - -XXX: This algorithm is derived from the Numerical Recipes. -XXX: use recycled vectors for internal data. -XXX: do an F64 version? - *****************************************************************************/ -static psF32 *calculateSecondDerivs(const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y) ///< Coordinates -{ - psTrace(".psLib.dataManip.calculateSecondDerivs", 4, - "---- calculateSecondDerivs() begin ----\n"); - - psS32 i; - psS32 k; - psF32 sig; - psF32 p; - psS32 n = y->n; - psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); - psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); - psF32 *X = (psF32 *) & (x->data.F32[0]); - psF32 *Y = (psF32 *) & (y->data.F32[0]); - psF32 qn; - - // XXX: The second derivatives at the endpoints, undefined in the SDR, - // are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV. - derivs2[0] = -0.5; - u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - PS_LEFT_SPLINE_DERIV); - - for (i=1;i<=(n-2);i++) { - sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]); - p = sig * derivs2[i-1] + 2.0; - derivs2[i] = (sig - 1.0) / p; - u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1])); - u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p; - - psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "X[%d] is %f\n", i, X[i]); - psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "Y[%d] is %f\n", i, Y[i]); - psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "u[%d] is %f\n", i, u[i]); - } - - qn = 0.5; - u[n-1] = (3.0/(X[n-1]-X[n-2])) * (PS_RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2])); - derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0); - - for (k=(n-2);k>=0;k--) { - derivs2[k] = derivs2[k] * derivs2[k+1] + u[k]; - - psTrace(".psLib.dataManip.calculateSecondDerivs", 6, - "derivs2[%d] is %f\n", k, derivs2[k]); - } - - psFree(u); - psTrace(".psLib.dataManip.calculateSecondDerivs", 4, - "---- calculateSecondDerivs() end ----\n"); - return(derivs2); -} - -/****************************************************************************** -p_psNRSpline1DEval(): This routine does NR-style evaluation of cubic splines. -It takes advantage of the 2nd derivatives of the cubic splines, which are -stored in the psSPline1D data structure, and computes the interpolated value -directly, without computing (or using) the interpolating cubic spline -polynomial. - -This routine is here mostly for a sanity check on the psLib function -evalSpline() which computes the interpolated value based on the cubic spline -polynomials which are stored in psSpline1D. - -XXX: This is F32 only - -XXX: spline->knots must be psF32 - -XXXX: Remove this for next code shipment. - *****************************************************************************/ -/* -psF32 p_psNRSpline1DEval(psSpline1D *spline, - const psVector* x, - const psVector* y, - psF32 X) -{ - PS_ASSERT_PTR_NON_NULL(spline, NAN); - PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN); - PS_ASSERT_VECTOR_NON_NULL(spline->domains, NAN); - PS_ASSERT_PTR_NON_NULL(spline->p_psDeriv2, NAN); - PS_ASSERT_VECTOR_NON_NULL(x, NAN); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_NON_NULL(y, NAN); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_TYPE(spline->knots, PS_TYPE_F32, NULL); - - psS32 n; - psS32 klo; - psS32 khi; - psF32 H; - psF32 A; - psF32 B; - psF32 C; - psF32 D; - psF32 Y; - - n = spline->n; - klo = p_psVectorBinDisect32(spline->knots->data.F32, (spline->n)+1, X); - if (klo < 0) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: psMinimize.c: p_psNRSpline1DEval(): p_psVectorBinDisect32 returned an error (%d).\n", klo); - return(NAN); - } - khi = klo + 1; - H = spline->knots->data.F32[khi] - spline->knots->data.F32[klo]; - A = (spline->knots->data.F32[khi] - X) / H; - B = (X - spline->knots->data.F32[klo]) / H; - C = ((A*A*A)-A) * (H*H/6.0); - D = ((B*B*B)-B) * (H*H/6.0); - - Y = (A * y->data.F32[klo]) + - (B * y->data.F32[khi]) + - (C * (spline->p_psDeriv2)[klo]) + - (D * (spline->p_psDeriv2)[khi]); - - return(Y); -} -*/ -/*****************************************************************************/ -/* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ - -/***************************************************************************** -psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y -vectors, this routine generates the linear or cublic splines which satisfy -those data points. - -The formula for calculating the spline polynomials is derived from Numerical -Recipes in C. The basic idea is that the polynomial is - (1) y = (A * y[0]) + - (2) (B * y[1]) + - (3) ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 + - (4) ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0 -Where: - H = x[1]-x[0] - A = (x[1]-x)/H - B = (x-x[0])/H -The bulk of the code in this routine is the expansion of the above equation -into a polynomial in terms of x, and then saving the coefficients of the -powers of x in the spline polynomials. This gets pretty complicated. - -XXX: usage of yErr is not specified in IfA documentation. - -XXX: Is the x argument redundant? What do we do if the x argument is -supplied, but does not equal the knots specified in mySpline? - -XXX: can psSpline be NULL? - -XXX: reimplement this assuming that mySpline is NULL? - -XXX: What happens if X is NULL, then an index vector is generated for X, but -that index vector lies outside the range vectors in mySpline? - -XXX: Assumes mySpline->knots is psF32. Must add psU32 and psF64. - *****************************************************************************/ -psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. - const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y, ///< Coordinates - const psVector* yErr) ///< Errors in coordinates, or NULL -{ - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - if (mySpline != NULL) { - PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL); - } - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() begin ----\n"); - psS32 numSplines = (y->n)-1; - psF32 tmp; - psF32 H; - psS32 i; - psF32 slope; - psVector *x32 = NULL; - psVector *y32 = NULL; - psVector *yErr32 = NULL; - static psVector *x32Static = NULL; - static psVector *y32Static = NULL; - static psVector *yErr32Static = NULL; - - PS_VECTOR_CONVERT_F64_TO_F32_STATIC(y, y32, y32Static); - - // If yErr==NULL, set all errors equal. - if (yErr == NULL) { - PS_VECTOR_GEN_YERR_STATIC_F32(yErr32Static, y->n); - yErr32 = yErr32Static; - } else { - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL); - PS_VECTOR_CONVERT_F64_TO_F32_STATIC(yErr, yErr32, yErr32Static); - } - - // If x==NULL, create an x32 vector with x values set to (0:n). - if (x == NULL) { - PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n); - x32 = x32Static; - } else { - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - PS_VECTOR_CONVERT_F64_TO_F32_STATIC(x, x32, x32Static); - } - PS_ASSERT_VECTORS_SIZE_EQUAL(x32, y32, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(yErr32, y32, NULL); - - /* - XXX: - This can not be implemented until SDR states what order spline should be - created. - Should we error if mySpline is not NULL? - Should we error if mySPline is not NULL? - */ - if (mySpline == NULL) { - mySpline = psSpline1DAllocGeneric(x32, 3); - } - PS_ASSERT_PTR_NON_NULL(mySpline, NULL); - PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL); - - if (y32->n != (1 + mySpline->n)) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "data size / spline size mismatch (%d %d)\n", - y32->n, mySpline->n); - return(NULL); - } - - // If these are linear splines, which means their polynomials will have - // two coefficients, then we do the simple calculation. - if (2 == (mySpline->spline[0])->n) { - for (i=0;in;i++) { - slope = (y32->data.F32[i+1] - y32->data.F32[i]) / - (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]); - (mySpline->spline[i])->coeff[0] = y32->data.F32[i] - - (slope * mySpline->knots->data.F32[i]); - - (mySpline->spline[i])->coeff[1] = slope; - psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- mySpline %d coeffs are (%f, %f)\n", i, - (mySpline->spline[i])->coeff[0], - (mySpline->spline[i])->coeff[1]); - } - psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- Exiting psVectorFitSpline1D()()\n"); - return((psSpline1D *) mySpline); - } - - // Check if these are cubic splines (n==4). If not, psError. - if (4 != (mySpline->spline[0])->n) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Don't know how to generate %d-order splines.", - (mySpline->spline[0])->n-1); - return(NULL); - } - - // If we get here, then we know these are cubic splines. We first - // generate the second derivatives at each data point. - mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32); - for (i=0;in;i++) - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); - - // We generate the coefficients of the spline polynomials. I can't - // concisely explain how this code works. See above function comments - // and Numerical Recipes in C. - for (i=0;idata.F32[i+1] - x32->data.F32[i]; - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "x data (%f - %f) (%f)\n", - x32->data.F32[i], - x32->data.F32[i+1], H); - // - // ******** Calculate 0-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H); - // From (2) - ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H); - // From (3) - tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); - tmp-= (x32->data.F32[i+1] / H); - tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - // From (4) - tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); - tmp+= (x32->data.F32[i] / H); - tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - - // - // ******** Calculate 1-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H; - // From (2) - ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H); - // From (3) - tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); - tmp+= (1.0 / H); - tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; - ((mySpline->spline[i])->coeff[1])+= tmp; - // From (4) - tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); - tmp-= (1.0 / H); - tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; - ((mySpline->spline[i])->coeff[1])+= tmp; - - // - // ******** Calculate 2-order term ******** - // - // From (3) - (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H); - // From (4) - ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H)); - - // - // ******** Calculate 3-order term ******** - // - // From (3) - (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); - // From (4) - ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); - - } - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() end ----\n"); - return(mySpline); -} - - - - - - - - - -/***************************************************************************** -psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y - -xF32 and yF32 are internal psVectors which are used to hold the psF32 versions -of the input data, if necessary. xPtr and yPtr are pointers to either xF32 or -the x argument. All computation is done on xPtr and yPtr. xF32 and yF32 will -simply be psFree() at the end. - -XXX: nKnots makes no sense. This number is always equal to the size of the x -an y vectors. - -XXX: How do we specify the spline order? For now, order=3. - *****************************************************************************/ -#define PS_XXX_SPLINE_ORDER 3 -psSpline1D *psVectorFitSpline1DNEW(const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y, ///< Coordinates - int nKnots) -{ - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n"); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - // PS_ASSERT_INT_EQUAL(y->n, nKnots); - - // - // The following code ensures that xPtr points to a psF32 version of the - // ordinate data. - // - psVector *xF32 = NULL; - psVector *xPtr = NULL; - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); - if (PS_TYPE_F64 == x->type.type) { - xF32 = psVectorAlloc(y->n, PS_TYPE_F32); - for (psS32 i = 0 ; i < x->n ; i++) { - xF32->data.F32[i] = (psF32) x->data.F64[i]; - } - xPtr = xF32; - } else if (PS_TYPE_F32 == x->type.type) { - xPtr = (psVector *) x; - } else { - printf("XXX: Gen Error message: x is wrong type.\n"); - return(NULL); - } - } else { - // If x==NULL, create an x32 vector with x values set to (0:n). - xF32 = psVectorAlloc(y->n, PS_TYPE_F32); - for (psS32 i = 0 ; i < x->n ; i++) { - xF32->data.F32[i] = (psF32) i; - } - xPtr = xF32; - } - - // - // If y is of type psF64, then create a new vector yF32 and convert the - // y elements. Regardless of y's type, we create a yPtr which will be - // used in the remainder of this function. - // - psVector *yF32 = NULL; - psVector *yPtr = NULL; - if (PS_TYPE_F64 == y->type.type) { - yF32 = psVectorAlloc(y->n, PS_TYPE_F32); - for (psS32 i = 0 ; i < y->n ; i++) { - yF32->data.F32[i] = (psF32) y->data.F64[i]; - } - yPtr = yF32; - } else { - yPtr = (psVector *) y; - } - - psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER); - - psS32 numSplines = nKnots - 1; - psF32 tmp; - psF32 H; - psS32 i; - psF32 slope; - // XXX: get rid of x32 and y32 (this is from old code) - psVector *x32 = xPtr; - psVector *y32 = yPtr; - - // If these are linear splines, which means their polynomials will have - // two coefficients, then we do the simple calculation. - if (1 == PS_XXX_SPLINE_ORDER) { - for (i=0;in;i++) { - slope = (y32->data.F32[i+1] - y32->data.F32[i]) / - (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]); - (mySpline->spline[i])->coeff[0] = y32->data.F32[i] - - (slope * mySpline->knots->data.F32[i]); - - (mySpline->spline[i])->coeff[1] = slope; - psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- mySpline %d coeffs are (%f, %f)\n", i, - (mySpline->spline[i])->coeff[0], - (mySpline->spline[i])->coeff[1]); - } - psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, - "---- Exiting psVectorFitSpline1D()()\n"); - return((psSpline1D *) mySpline); - } - - // - // Check if these are cubic splines (n==4). If not, psError. - // - if (3 != PS_XXX_SPLINE_ORDER) { - psError(PS_ERR_BAD_PARAMETER_SIZE, true, - "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER); - return(NULL); - } - - // - // If we get here, then we know these are cubic splines. We first - // generate the second derivatives at each data point. - // - mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32); - for (i=0;in;i++) - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); - - // - // We generate the coefficients of the spline polynomials. I can't - // concisely explain how this code works. See above function comments - // and Numerical Recipes in C. - // - for (i=0;idata.F32[i+1] - x32->data.F32[i]; - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "x data (%f - %f) (%f)\n", - x32->data.F32[i], - x32->data.F32[i+1], H); - // - // ******** Calculate 0-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H); - // From (2) - ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H); - // From (3) - tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); - tmp-= (x32->data.F32[i+1] / H); - tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - // From (4) - tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); - tmp+= (x32->data.F32[i] / H); - tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; - ((mySpline->spline[i])->coeff[0])+= tmp; - - // - // ******** Calculate 1-order term ******** - // - // From (1) - (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H; - // From (2) - ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H); - // From (3) - tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); - tmp+= (1.0 / H); - tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; - ((mySpline->spline[i])->coeff[1])+= tmp; - // From (4) - tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); - tmp-= (1.0 / H); - tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; - ((mySpline->spline[i])->coeff[1])+= tmp; - - // - // ******** Calculate 2-order term ******** - // - // From (3) - (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H); - // From (4) - ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H)); - - // - // ******** Calculate 3-order term ******** - // - // From (3) - (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); - // From (4) - ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); - psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, - "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); - - } - - psFree(xF32); - psFree(yF32); - - psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, - "---- psVectorFitSpline1D() end ----\n"); - - return(mySpline); -} - - - - -//XXX: What's this for? -psF64 p_psImageGetElementF64(psImage *a, int i, int j); /****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates @@ -717,4 +132,6 @@ will have to convert all F32 input vectors to F64 regardless. So, the F64 port might be. + +XXX: Change the whole thing to F64, if input data is F32, convert it. *****************************************************************************/ psBool psMinimizeLMChi2(psMinimization *min, @@ -816,5 +233,5 @@ if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) { // XXX: Where is this? - // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); + // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess"); } # endif /* PS_NO_TRACE */ @@ -876,5 +293,4 @@ // XXX EAM: this needs to respect the mask on params -// XXX EAM: check not NULL on alpha, beta, params // alpha, beta, params are already allocated psF64 p_psMinLM_SetABX (psImage *alpha, @@ -887,4 +303,10 @@ psMinimizeLMChi2Func func) { + PS_ASSERT_IMAGE_NON_NULL(alpha, NAN); + PS_ASSERT_VECTOR_NON_NULL(beta, NAN); + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_NON_NULL(x, NAN); + PS_ASSERT_VECTOR_NON_NULL(y, NAN); + PS_ASSERT_VECTOR_NON_NULL(dy, NAN); psF64 chisq; @@ -1011,7 +433,4 @@ } -// XXX: put this in constants.h -# define PS_SWAP(X,Y) {double tmp=(X); (X) = (Y); (Y) = tmp;} - // XXX EAM : temporary gauss-jordan solver based on gene's // version based on the Numerical Recipes version @@ -1044,6 +463,5 @@ for (j = 0; j < Nx; j++) { if (!isfinite(matrix[i][j])) { - // XXX EAM: this should use the psError stack - fprintf (stderr, "GAUSSJ: NaN\n"); + psError(PS_ERR_UNKNOWN, false, "Input matrix contains NaNs.\n"); goto fescape; } @@ -1058,6 +476,5 @@ } else { if (ipiv[k] > 1) { - // XXX EAM: this should use the psError stack - fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); + psError(PS_ERR_UNKNOWN, false, "Singular Matrix (1).\n"); goto fescape; } @@ -1076,6 +493,5 @@ indxc[i] = icol; if (matrix[icol][icol] == 0.0) { - // XXX EAM: this should use the psError stack - fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); + psError(PS_ERR_UNKNOWN, false, "Singular Matrix (2).\n"); goto fescape; } @@ -1118,266 +534,7 @@ } -/****************************************************************************** -vectorFitPolynomial1DCheb(): This routine will fit a Chebyshev polynomial of -degree myPoly to the data points (x, y) and return the coefficients of that -polynomial. - -XXX: yErr is currently ignored. -*****************************************************************************/ -static psPolynomial1D *vectorFitPolynomial1DCheby(psPolynomial1D* myPoly, - const psVector* x, - const psVector* y, - const psVector* yErr) -{ - psS32 j; - psS32 k; - psS32 n = x->n; - psF64 fac; - psF64 sum; - PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); - psScalar *fScalar; - psScalar tmpScalar; - tmpScalar.type.type = PS_TYPE_F64; - - // XXX: These assignments appear too simple to warrant code and - // variable declarations. I retain them here to maintain coherence - // with the NR code. - psF64 min = -1.0; - psF64 max = 1.0; - psF64 bma = 0.5 * (max-min); // 1 - psF64 bpa = 0.5 * (max+min); // 0 - - // In this loop, we first calculate the values of X for which the - // Chebyshev polynomials are zero (see NR, section 5.4). Then we - // calculate the value of the function we are fitting the Chebyshev - // polynomials to at those values of X. This is a bit tricky since - // we don't know that function. So, we instead do 3-order LaGrange - // interpolation at the point X for the psVectors x,y for which we - // are fitting this ChebyShev polynomial to. - - for (psS32 i=0;idata.F64[i] = fScalar->data.F64; - psFree(fScalar); - - psTrace(".psLib.dataManip.vectorFitPolynomial1DCheby", 6, - "(x, X, y, f(X)) is (%f, %f, %f, %f)\n", - x->data.F64[i], X, y->data.F64[i], f->data.F64[i]); - } - - // We have the values for f() at the zero points, we now calculate the - // coefficients of the Chebyshev polynomial: NR 5.8.7. - - fac = 2.0/((psF32) n); - // XXX: is this loop bound correct? - for (j=0;jn;j++) { - sum = 0.0; - for (k=0;kdata.F64[k] * - cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); - } - - myPoly->coeff[j] = fac * sum; - } - - return(myPoly); -} - -/****************************************************************************** -VectorFitPolynomial1DOrd(): This routine will fit an ordinary polynomial of -degree myPoly to the data points (x, y) and return the coefficients of that -polynomial. - -XXX: Use private name? -XXX: Use recycled vectors. - *****************************************************************************/ -static psPolynomial1D* vectorFitPolynomial1DOrd(psPolynomial1D* myPoly, - const psVector* x, - const psVector* y, - const psVector* yErr) -{ - psS32 polyOrder = myPoly->n; - psImage* A = NULL; - psImage* ALUD = NULL; - psVector* B = NULL; - psVector* outPerm = NULL; - psVector* X = NULL; // NOTE: do we need this? - psVector* coeffs = NULL; - psS32 i = 0; - psS32 j = 0; - psS32 k = 0; - psVector* xSums = NULL; - - psTrace(".psLib.dataManip.vectorFitPolynomial1DOrd", 4, - "---- vectorFitPolynomial1DOrd() begin ----\n"); - // printf("VectorFitPolynomial1D()\n"); - // for (i=0;in;i++) { - // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]); - // } - - A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); - ALUD = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); - - B = psVectorAlloc(polyOrder, PS_TYPE_F64); - coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64); - X = psVectorAlloc(x->n, PS_TYPE_F64); - xSums = psVectorAlloc(1 + 2 * polyOrder, PS_TYPE_F64); - - // Initialize data structures. - for (i = 0; i < polyOrder; i++) { - B->data.F64[i] = 0.0; - coeffs->data.F64[i] = 0.0; - for (j = 0; j < polyOrder; j++) { - A->data.F64[i][j] = 0.0; - ALUD->data.F64[i][j] = 0.0; - } - } - for (i = 0; i < X->n; i++) { - X->data.F64[i] = x->data.F64[i]; - } - - // Build the B and A data structs. - if (yErr == NULL) { - for (i = 0; i < X->n; i++) { - buildSums1D(X->data.F64[i], 2 * polyOrder, xSums); - - for (k = 0; k < polyOrder; k++) { - B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k]; - } - - for (k = 0; k < polyOrder; k++) { - for (j = 0; j < polyOrder; j++) { - A->data.F64[k][j] += xSums->data.F64[k + j]; - } - } - } - } else { - for (i = 0; i < X->n; i++) { - buildSums1D(X->data.F64[i], 2 * polyOrder, xSums); - - for (k = 0; k < polyOrder; k++) { - B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k] / - yErr->data.F64[i]; - } - - for (k = 0; k < polyOrder; k++) { - for (j = 0; j < polyOrder; j++) { - A->data.F64[k][j] += xSums->data.F64[k + j] / - yErr->data.F64[i]; - } - } - } - } - - // XXX: How do we know if these routines were successful? - ALUD = psMatrixLUD(ALUD, &outPerm, A); - coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - - for (k = 0; k < polyOrder; k++) { - myPoly->coeff[k] = coeffs->data.F64[k]; - // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]); - } - - psFree(A); - psFree(ALUD); - psFree(B); - psFree(coeffs); - psFree(X); - psFree(outPerm); - psFree(xSums); - - psTrace(".psLib.dataManip.vectorFitPolynomial1DOrd", 4, - "---- vectorFitPolynomial1DOrd() begin ----\n"); - return (myPoly); -} - -/****************************************************************************** -psVectorFitPolynomial1D(): This routine must fit a polynomial of degree -myPoly to the data points (x, y) and return the coefficients of that -polynomial. - -XXX: type F32 is done via vector conversion only. - -XXX: Get rid of this. Use new argument list. - *****************************************************************************/ -psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* poly, - const psVector* x, - const psVector* y, - const psVector* yErr) -{ - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(y, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); - - psS32 i; - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *yErr64 = NULL; - static psVector *x64Static = NULL; - static psVector *y64Static = NULL; - static psVector *yErr64Static = NULL; - - PS_VECTOR_CONVERT_F32_TO_F64_STATIC(y, y64, y64Static); - // If yErr==NULL, set all errors equal. - if (yErr == NULL) { - PS_VECTOR_GEN_YERR_STATIC_F64(yErr64Static, y->n); - yErr64 = yErr64Static; - } else { - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL); - PS_VECTOR_CONVERT_F32_TO_F64_STATIC(yErr, yErr64, yErr64Static); - } - - // If x==NULL, create an x64 vector with x values set to (0:n). - if (x == NULL) { - PS_VECTOR_GEN_X_INDEX_STATIC_F64(x64Static, y->n); - if (poly->type == PS_POLYNOMIAL_CHEB) { - p_psNormalizeVectorRangeF64(x64Static, -1.0, 1.0); - } - x64 = x64Static; - } else { - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - PS_VECTOR_CONVERT_F32_TO_F64_STATIC(x, x64, x64Static); - if (poly->type == PS_POLYNOMIAL_CHEB) { - p_psNormalizeVectorRangeF64(x64, -1.0, 1.0); - } - } - PS_ASSERT_VECTORS_SIZE_EQUAL(x64, y64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(yErr64, y64, NULL); - - // Call the appropriate vector fitting routine. - psPolynomial1D *rc = NULL; - if (poly->type == PS_POLYNOMIAL_CHEB) { - rc = vectorFitPolynomial1DCheby(poly, x64, y64, yErr64); - } else if (poly->type == PS_POLYNOMIAL_ORD) { - rc = vectorFitPolynomial1DOrd(poly, x64, y64, yErr64); - } else { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - "unknown polynomial type.\n"); - return(NULL); - } - if (rc == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - - return(poly); -} - static void minimizationFree(psMinimization *min) { - // There are non dynamic allocated items + // There are no dynamically allocated items } @@ -1437,4 +594,10 @@ /****************************************************************************** + ****************************************************************************** +Powell routines. +****************************************************************************** +*****************************************************************************/ + +/****************************************************************************** p_psDetermineBracket(): This routine takes as input an arbitrary function, and the parameter to vary, and the line along which it must vary. This @@ -1446,21 +609,17 @@ Algorithm: -XXX completely ad hoc: -start with the user-supplied starting parameter and +XXX completely ad hoc: start with the user-supplied starting parameter and call that b. Calculate a/c as a fractional amount smaller/larger than b. Repeat this process until a local minimum is found. -XXX: -new algorithm: -start at x=0, expand in one direction until the function +XXX: new algorithm: start at x=0, expand in one direction until the function decreases. Then you have two points in the bracket. Keep going until it increases, or x is too large. If thst does not work, expand in the other direction. -XXX: -This is F32 only. +XXX: This is F32 only. Must add F64 support (actually, make the defaults F64, +and convert F32 vectors to F64). -XXX: -output bracket vector should be an input as well. +XXX: output bracket vector should be an input as well. *****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, @@ -1716,5 +875,5 @@ XXX: This routine is not very efficient in terms of total evaluations of the function. -XXX: This is F32 only +XXX: This is F32 only (make it F64). XXX: Since this is an internal function, many of the parameter checks are redundant. @@ -1871,6 +1030,5 @@ XXX: The SDR is silent about data types. F32 is implemented here. - -XXX: Check for F32 types? +Reimplement in F64, convert F32 vectors to F64. *****************************************************************************/ #define PS_MINIMIZE_POWELL_LINEMIN_MAX_ITERATIONS 20 @@ -2083,5 +1241,5 @@ This functions uses global variables to receive the function pointer, the data values, and the data errors. -XXX: This is F32 only +XXX: This is F32 only. Must implement F64. *****************************************************************************/ static psF32 myPowellChi2Func(const psVector *params, @@ -2149,13 +1307,25 @@ /****************************************************************************** ****************************************************************************** - ****************************************************************************** - ****************************************************************************** -EAM Code: - ****************************************************************************** - ****************************************************************************** + Analytical 1-D fitting routines. ****************************************************************************** *****************************************************************************/ - -static psVector *EAMBuildSums1D( +// XXX: Make this a general type conversion macro, or function +#define PS_VECTOR_GEN_F64_FROM_F32(VECF64, VECF32) \ +VECF64 = psVectorAlloc(VECF32->n, PS_TYPE_F64); \ +for (psS32 i = 0 ; i < VECF32->n ; i++) { \ + VECF64->data.F64[i] = (psF64) VECF32->data.F32[i]; \ +} \ + +#define PS_VECTOR_GEN_CHEBY_INDEX(VECF64, SIZE) \ +VECF64 = psVectorAlloc(SIZE, PS_TYPE_F64); \ +for (psS32 i = 0 ; i < SIZE ; i++) { \ + VECF64->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \ +}\ +/****************************************************************************** +BuildSums1D(sums, x, polyOrder, sums): this routine calculates the powers of +input parameter "x" between 0 and input parameter nTerms*2. The result is +returned as a psVector sums. +*****************************************************************************/ +static psVector *BuildSums1D( psVector* sums, psF64 x, @@ -2166,6 +1336,6 @@ // - // XXX: Why do we multiply by 2 here? Better to do it outside and have the - // definition of this function remain sensible. + // XXX: Why do we multiply by 2 here? It's better to do it outside and + // have the definition of this function remain sensible. // nSum = 2*nTerm; @@ -2184,13 +1354,218 @@ } -// XXX EAM : test version of 1d fitting -psPolynomial1D* VectorFitPolynomial1DOrd_EAM( +/****************************************************************************** +BuildSums2D(sums, x, y, nXterm, nYterm): this routine calculates the powers of +input parameter "x" and "y" between 0 and input parameter nXterms*2 and +nYterm*2. The result is returned as a psImage sums. + *****************************************************************************/ +static psImage *BuildSums2D( + psImage *sums, + psF64 x, + psF64 y, + psS32 nXterm, + psS32 nYterm) +{ + psS32 nXsum = 0; + psS32 nYsum = 0; + psF64 xSum = 1.0; + psF64 ySum = 1.0; + + nXsum = 2*nXterm; + nYsum = 2*nYterm; + if (sums == NULL) { + sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64); + } + if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) { + psFree (sums); + sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64); + } + + ySum = 1.0; + for (int j = 0; j < nYsum; j++) { + xSum = ySum; + for (int i = 0; i < nXsum; i++) { + sums->data.F64[i][j] = xSum; + xSum *= x; + } + ySum *= y; + } + return (sums); +} + + +/****************************************************************************** +Polynomial2DEvalVectorD(myPoly, x, y): This routine takes as input two +psVectors x and y, and evaluates myPoly for each pair of (x, y), and stores it +in the output vector. This routine works on single-precision polynomials with +double precision data. + *****************************************************************************/ +psVector *Polynomial2DEvalVectorD( + const psPolynomial2D *myPoly, + const psVector *x, + const psVector *y) +{ + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + + // + // Set the length of the output vector to the minimum of the x,y vectors. + // + psS32 vecLen=x->n; + if (y->n < vecLen) { + vecLen = y->n; + } + + // + // Create output vector to return + // + psVector *tmp = psVectorAlloc(vecLen, PS_TYPE_F64); + + // + // Evaluate the polynomial at the specified points + // + for (psS32 i=0; idata.F64[i] = psPolynomial2DEval(myPoly, x->data.F64[i], y->data.F64[i]); + } + + return(tmp); +} + +/****************************************************************************** + ****************************************************************************** + 1-D Vector Fitting Code. + ****************************************************************************** + *****************************************************************************/ + +/****************************************************************************** +vectorFitPolynomial1DCheb(): This routine will fit a Chebyshev polynomial of +degree myPoly to the data points (x, y) and return the coefficients of that +polynomial. + +XXX: mask, maskValue, yErr are currently ignored. + +XXX: Change arg order to that of the psLib function. +*****************************************************************************/ +static psPolynomial1D *vectorFitPolynomial1DCheby( psPolynomial1D* myPoly, - psVector *mask, - const psVector *x, - const psVector *y, - const psVector *yErr) -{ - // I think this is 1 dimension down + const psVector *mask, + psMaskType maskValue, + const psVector* y, + const psVector* yErr, + const psVector* x) +{ + // XXX: these ASSERTS are redundant. + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + if (yErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, NULL); + PS_ASSERT_VECTOR_TYPE(yErr, PS_TYPE_F64, NULL); + } + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + } + + psS32 j; + psS32 k; + psS32 n = y->n; + psF64 fac; + psF64 sum; + PS_VECTOR_GEN_STATIC_RECYCLED(f, n, PS_TYPE_F64); + psScalar *fScalar; + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F64; + + // XXX: These assignments appear too simple to warrant code and + // variable declarations. I retain them here to maintain coherence + // with the NR code. + psF64 min = -1.0; + psF64 max = 1.0; + psF64 bma = 0.5 * (max-min); // 1 + psF64 bpa = 0.5 * (max+min); // 0 + + // In this loop, we first calculate the values of X for which the + // Chebyshev polynomials are zero (see NR, section 5.4). Then we + // calculate the value of the function we are fitting the Chebyshev + // polynomials to at those values of X. This is a bit tricky since + // we don't know that function. So, we instead do 3-order LaGrange + // interpolation at the point X for the psVectors x,y for which we + // are fitting this ChebyShev polynomial to. + + for (psS32 i=0;idata.F64[i] = fScalar->data.F64; + psFree(fScalar); + + psTrace(".psLib.dataManip.vectorFitPolynomial1DCheby", 6, + "(x, X, y, f(X)) is (%f, %f, %f, %f)\n", + x->data.F64[i], X, y->data.F64[i], f->data.F64[i]); + } + + // We have the values for f() at the zero points, we now calculate the + // coefficients of the Chebyshev polynomial: NR 5.8.7. + + fac = 2.0/((psF32) n); + // XXX: is this loop bound correct? + for (j=0;jn;j++) { + sum = 0.0; + for (k=0;kdata.F64[k] * + cos(M_PI * ((psF32) j) * (0.5 + ((psF32) k)) / ((psF32) n)); + } + + myPoly->coeff[j] = fac * sum; + } + + return(myPoly); +} + +/****************************************************************************** +VectorFitPolynomial1DOrd(myPoly, *mask, maskValue, *y, *yErr, *x): This is a +private routine which will fit a 1-D polynomial to a set of (x, f) pairs. The +x and fErr vectors may be NULL. All non-NULL vectors must be of type +PS_TYPE_F64. + *****************************************************************************/ +psPolynomial1D* VectorFitPolynomial1DOrd( + psPolynomial1D* myPoly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x) +{ + // XXX: these ASSERTS are redundant. + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->n, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); + } + psImage* A = NULL; psVector* B = NULL; @@ -2204,10 +1579,20 @@ if (psTraceGetLevel (".psLib.dataManip.VectorFitPolynomial1DOrd") >= 5) { - FILE *f = fdopen((psTraceGetDestination ()), "a+"); - fprintf (f, "VectorFitPolynomial1D()\n"); - for (int i = 0; i < x->n; i++) { - fprintf (f, "(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]); - } - fclose(f); + FILE *fp = fdopen((psTraceGetDestination ()), "a+"); + fprintf(fp, "VectorFitPolynomial1D()\n"); + for (int i = 0; i < f->n; i++) { + fprintf(fp, "(x, f, fErr) is ("); + if (x != NULL) { + fprintf(fp, "%f, %f, ", x->data.F64[i], f->data.F64[i]); + } else { + fprintf(fp, "%f, %f, ", (psF64) i, f->data.F64[i]); + } + if (fErr != NULL) { + fprintf(fp, "%f)\n", fErr->data.F64[i]); + } else { + fprintf(fp, "NULL)\n"); + } + } + fclose(fp); } @@ -2227,17 +1612,21 @@ // xSums look like: 1, x, x^2, ... x^(2n+1) // Build the B and A data structs. - for (int k = 0; k < x->n; k++) { + for (int k = 0; k < f->n; k++) { if ((mask != NULL) && mask->data.U8[k]) continue; - xSums = EAMBuildSums1D(xSums, x->data.F64[k], nTerm); - - if (yErr == NULL) { + if (x != NULL) { + xSums = BuildSums1D(xSums, x->data.F64[k], nTerm); + } else { + xSums = BuildSums1D(xSums, (psF64) k, nTerm); + } + + if (fErr == NULL) { wt = 1.0; } else { - // this should filter yErr == 0 values - wt = 1.0 / PS_SQR(yErr->data.F64[k]); + // this should filter fErr == 0 values + wt = 1.0 / PS_SQR(fErr->data.F64[k]); } for (int i = 0; i < nTerm; i++) { - B->data.F64[i] += y->data.F64[k] * xSums->data.F64[i] * wt; + B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt; } @@ -2261,7 +1650,6 @@ myPoly->coeff[k] = B->data.F64[k]; } - } else + } else { // LUD version of the fit - { psImage *ALUD = NULL; psVector* outPerm = NULL; @@ -2274,4 +1662,7 @@ myPoly->coeff[k] = coeffs->data.F64[k]; } + psFree(ALUD); + psFree(coeffs); + psFree(outPerm); } @@ -2281,84 +1672,188 @@ psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4, - "---- VectorFitPolynomial1DOrd() begin ----\n"); + "---- VectorFitPolynomial1DOrd() End ----\n"); return (myPoly); } -// XXX EAM : this version uses the F64 vectors -psVector *Polynomial2DEvalVectorD( - const psPolynomial2D *myPoly, + +/****************************************************************************** +psVectorFitPolynomial1D(): This routine fits a polynomial of arbitrary degree +(specified in poly) to the data points (x, y) and return that polynomial. +Types F32 and F64 are supported, however, type F32 is done via vector +conversion only. + *****************************************************************************/ +psPolynomial1D *psVectorFitPolynomial1D( + psPolynomial1D *poly, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x) +{ + // Internal pointers for possibly NULL or mis-typed vectors. + psVector *x64 = NULL; + psVector *f64 = NULL; + psVector *fErr64 = NULL; + + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_INT_NONNEGATIVE(poly->n, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_NON_EMPTY(f, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); + if (f->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(f64, f); + } else { + f64 = (psVector *) f; + } + + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); + if (x->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(x64, x); + } else { + x64 = (psVector *) x; + } + } + + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + if (fErr->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr); + } else { + fErr64 = (psVector *) fErr; + } + } + + if (poly->type == PS_POLYNOMIAL_ORD) { + poly = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); + return(NULL); + } + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + if (mask != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + } + if (x == NULL) { + // If x==NULL, create an x64 vector with x values set to (-1:1). + PS_VECTOR_GEN_CHEBY_INDEX(x64, f64->n); + } + // XXX: Change arg order. + // XXX: Must modify this routine so that x64 or fErr64 can be NULL. + poly = vectorFitPolynomial1DCheby(poly, NULL, 0, f64, fErr64, x64); + if (x == NULL) { + psFree(x64); + } + } else { + psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); + } + + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if ((x != NULL) && (x->type.type != PS_TYPE_F64)) { + psFree(x64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + + return(NULL); +} + + +psPolynomial1D *psVectorClipFitPolynomial1D( + psPolynomial1D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x) +{ + // Internal pointers for possibly NULL vectors. + psVector *x32 = NULL; + + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_PTR_NON_NULL(stats, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } + + psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + // Free psVectors that were created for NULL arguments. + if (x == NULL) { + psFree(x32); + } + return(NULL); +} + + +/****************************************************************************** + ****************************************************************************** + 2-D Vector Fitting Code. + ****************************************************************************** + *****************************************************************************/ + +/****************************************************************************** +VectorFitPolynomial2DOrd(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. + +// XXX: What about the maskValue? + *****************************************************************************/ +psPolynomial2D* VectorFitPolynomial2DOrd( + psPolynomial2D* myPoly, + const psVector* mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, const psVector *x, const psVector *y) { + // These ASSERTS are redundant. PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); + + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - - psVector *tmp; - psS32 vecLen=x->n; - - // Determine the length of the output vector to by the minimum of the x,y vectors - if (y->n < vecLen) { - vecLen = y->n; - } - - // Create output vector to return - tmp = psVectorAlloc(vecLen, PS_TYPE_F64); - - // Evaluate the polynomial at the specified points - for (psS32 i=0; idata.F64[i] = psPolynomial2DEval(myPoly, x->data.F64[i], y->data.F64[i]); - } - - // Return output vector - return(tmp); -} - -// XXX EAM : EAMBuildSums2D in analogy with EAMBuildSums1D -static psImage *EAMBuildSums2D( - psImage *sums, - psF64 x, - psF64 y, - psS32 nXterm, - psS32 nYterm) -{ - psS32 nXsum = 0; - psS32 nYsum = 0; - psF64 xSum = 1.0; - psF64 ySum = 1.0; - - nXsum = 2*nXterm; - nYsum = 2*nYterm; - if (sums == NULL) { - sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64); - } - if ((nXsum != sums->numCols) || (nYsum != sums->numRows)) { - psFree (sums); - sums = psImageAlloc(nXsum, nYsum, PS_TYPE_F64); - } - - ySum = 1.0; - for (int j = 0; j < nYsum; j++) { - xSum = ySum; - for (int i = 0; i < nXsum; i++) { - sums->data.F64[i][j] = xSum; - xSum *= x; - } - ySum *= y; - } - return (sums); -} - -// XXX EAM : test version of 2d fitting -psPolynomial2D* VectorFitPolynomial2DOrd_EAM( - psPolynomial2D* myPoly, - psVector* mask, - const psVector* x, - const psVector* y, - const psVector* z, - const psVector* zErr) -{ + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + // I think this is 1 dimension down psImage* A = NULL; @@ -2389,13 +1884,13 @@ if ((mask != NULL) && mask->data.U8[k]) continue; - Sums = EAMBuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm); - - if (zErr == NULL) { + Sums = BuildSums2D(Sums, x->data.F64[k], y->data.F64[k], nXterm, nYterm); + + if (fErr == NULL) { wt = 1.0; } else { - // XXX: this should probably by zErr^2 !! - // this should filter zErr == 0 values - // XXX: Why isn't this zErr^2? - wt = 1.0 / zErr->data.F64[k]; + // XXX: this should probably by fErr^2 !! + // this should filter fErr == 0 values + // XXX: Why isn't this fErr^2? + wt = 1.0 / fErr->data.F64[k]; } @@ -2404,5 +1899,5 @@ for (int n = 0; n < nXterm; n++) { for (int m = 0; m < nYterm; m++) { - B->data.F64[n+m*nXterm] += z->data.F64[k] * Sums->data.F64[n][m] * wt; + B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt; } } @@ -2447,6 +1942,6 @@ const psVector* x, const psVector* y, - const psVector* z, - const psVector* dz) + const psVector* f, + const psVector* df) { psVector *X; @@ -2455,33 +1950,33 @@ psVector *dZ; - psVector *zFit = NULL; - psVector *zResid = NULL; + psVector *fFit = NULL; + psVector *fResid = NULL; psStats *stats = NULL; X = psVectorCopy (NULL, x, PS_TYPE_F64); Y = psVectorCopy (NULL, y, PS_TYPE_F64); - Z = psVectorCopy (NULL, z, PS_TYPE_F64); - dZ = psVectorCopy (NULL, dz, PS_TYPE_F64); + Z = psVectorCopy (NULL, f, PS_TYPE_F64); + dZ = psVectorCopy (NULL, df, PS_TYPE_F64); for (int N = 0; N < 3; N++) { // XXX EAM : this would be better defined with an element mask - poly = VectorFitPolynomial2DOrd_EAM(poly, NULL, X, Y, Z, dZ); - zFit = Polynomial2DEvalVectorD(poly, x, y); - zResid = (psVector *) psBinaryOp(NULL, (void *) z, "-", (void *) zFit); + poly = VectorFitPolynomial2DOrd(poly, NULL, 0, X, Y, Z, dZ); + fFit = Polynomial2DEvalVectorD(poly, x, y); + fResid = (psVector *) psBinaryOp(NULL, (void *) f, "-", (void *) fFit); stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); - stats = psVectorStats (stats, zResid, NULL, NULL, 0); + stats = psVectorStats (stats, fResid, NULL, NULL, 0); psTrace (".psphot.RobustFit", 4, "residual stats for robust fit: %g +/- %g (%d pts)\n", stats->clippedMean, stats->clippedStdev, stats->clippedNvalues); // re-create X, Y, Z, dZ if pts are valid int n = 0; - for (int i = 0; i < zResid->n; i++) { - if (fabs(zResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) { + for (int i = 0; i < fResid->n; i++) { + if (fabs(fResid->data.F64[i] - stats->clippedMean) > 3*stats->clippedStdev) { continue; } X->data.F64[n] = x->data.F64[i]; Y->data.F64[n] = y->data.F64[i]; - Z->data.F64[n] = z->data.F64[i]; - dZ->data.F64[n] = dz->data.F64[i]; + Z->data.F64[n] = f->data.F64[i]; + dZ->data.F64[n] = df->data.F64[i]; n++; } @@ -2507,26 +2002,26 @@ psPolynomial2D* RobustFit2D(psPolynomial2D* poly, - psVector* mask, + const psVector* mask, const psVector* x, const psVector* y, - const psVector* z, - const psVector* dz) + const psVector* f, + const psVector* df) { PS_ASSERT_VECTOR_NON_NULL(mask, NULL); PS_ASSERT_VECTOR_NON_NULL(x, NULL); PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_NON_NULL(dz, NULL); - - psVector *zFit = NULL; - psVector *zResid = psVectorAlloc (x->n, PS_TYPE_F64); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_NON_NULL(df, NULL); + + psVector *fFit = NULL; + psVector *fResid = psVectorAlloc (x->n, PS_TYPE_F64); psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); for (int N = 0; N < 3; N++) { - poly = VectorFitPolynomial2DOrd_EAM (poly, mask, x, y, z, dz); - zFit = Polynomial2DEvalVectorD (poly, x, y); - zResid = (psVector *) psBinaryOp(zResid, (void *) z, "-", (void *) zFit); - - stats = psVectorStats (stats, zResid, NULL, mask, 1); + poly = VectorFitPolynomial2DOrd (poly, mask, 0, x, y, f, df); + fFit = Polynomial2DEvalVectorD (poly, x, y); + fResid = (psVector *) psBinaryOp(fResid, (void *) f, "-", (void *) fFit); + + stats = psVectorStats (stats, fResid, NULL, mask, 1); psTrace (".psphot.RobustFit", 4, "residual stats for robust fit: %g +/- %g\n", stats->sampleMean, stats->sampleStdev); @@ -2535,92 +2030,28 @@ // we are masking out any point which is out of range // recovery is not allowed with this scheme - for (int i = 0; i < zResid->n; i++) { + for (int i = 0; i < fResid->n; i++) { if (mask->data.U8[i]) continue; - if (fabs(zResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) { + if (fabs(fResid->data.F64[i] - stats->sampleMean) > 3*stats->sampleStdev) { mask->data.U8[i] = 1; continue; } } - psFree (zFit); - } - psFree (zResid); + psFree (fFit); + } + psFree (fResid); psFree (stats); return (poly); } + /****************************************************************************** - ****************************************************************************** - ****************************************************************************** - ****************************************************************************** -NEW Code: - ****************************************************************************** - ****************************************************************************** - ****************************************************************************** +psVectorFitPolynomial2D(): This routine fits a 2D polynomial of arbitrary +degree (specified in poly) to the data points (x, y)-(f) and returns that +polynomial. Types F32 and F64 are supported, however, type F32 is done via +vector conversion only. *****************************************************************************/ - -#define PS_VECTOR_GEN_INDEX_F32(VEC, SIZE) \ -VEC = psVectorAlloc(SIZE, PS_TYPE_F32); \ -for (psS32 i = 0 ; i < SIZE ; i++) { \ - VEC->data.F32[i] = (psF32) i; \ -} \ - -psPolynomial1D *psVectorFitPolynomial1D_NEW( - psPolynomial1D *poly, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x) -{ - // Internal pointers for possibly NULL vectors. - psVector *x32 = NULL; - psVector *fErr32 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); - } else { - PS_VECTOR_GEN_INDEX_F32(x32, f->n); - } - - if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } else { - fErr = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr, 1.0); - } - - psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); - if (poly->type == PS_POLYNOMIAL_ORD) { - // XXX: Call EAM code - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - // XXX: Call my code - } else { - printf("XXX: ERROR: incorrect polynomial type.\n"); - } - - // Free psVectors that were created for NULL arguments. - if (x == NULL) { - psFree(x32); - } - if (fErr == NULL) { - psFree(fErr32); - } - - return(NULL); -} - - -psPolynomial2D *psVectorFitPolynomial2D_NEW( - psPolynomial1D *poly, +psPolynomial2D *psVectorFitPolynomial2D( + psPolynomial2D *poly, const psVector *mask, psMaskType maskValue, @@ -2630,9 +2061,146 @@ const psVector *y) { - // Internal pointers for possibly NULL vectors. - psVector *fErr32 = NULL; + // Internal pointers for possibly NULL or mis-typed vectors. + psVector *x64 = NULL; + psVector *y64 = NULL; + psVector *f64 = NULL; + psVector *fErr64 = NULL; PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + + // + // f + // + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); + if (f->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(f64, f); + } else { + f64 = (psVector *) f; + } + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + + // + // x + // + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + if (x->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(x64, x); + } else { + x64 = (psVector *) x; + } + + // + // y + // + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + if (y->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(y64, y); + } else { + y64 = (psVector *) y; + } + + // + // fErr + // + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + if (fErr->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr); + } else { + fErr64 = (psVector *) fErr; + } + } + + if (poly->type == PS_POLYNOMIAL_ORD) { + poly = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64); + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + return(NULL); + } + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + if (mask != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + psFree(poly); + poly = NULL; + } else { + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); + } + + + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + + return(poly); +} + +psPolynomial2D *psVectorClipFitPolynomial2D( + psPolynomial2D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y) +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_PTR_NON_NULL(stats, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); @@ -2648,27 +2216,87 @@ PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + } + + if (mask == NULL) { + // XXX: Change argument order. + poly = RobustFit2D_nomask(poly, x, y, f, fErr); } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); - } - - psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); - if (poly->type == PS_POLYNOMIAL_ORD) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D polynomial vector fitting has not been implemented. Returning NULL.\n"); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - } else { - printf("XXX: ERROR: incorrect polynomial type.\n"); - } - if (fErr == NULL) { - psFree(fErr32); - } - return(NULL); -} - -psPolynomial3D *psVectorFitPolynomial3D_NEW( - psPolynomial1D *poly, + // XXX: Use maskValue. + // XXX: Change argument order. + poly = RobustFit2D(poly, mask, x, y, f, fErr); + } + + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); + return(NULL); + } + return(poly); +} + +/****************************************************************************** + ****************************************************************************** + 3-D Vector Fitting Code. + ****************************************************************************** + *****************************************************************************/ + +/****************************************************************************** +VectorFitPolynomial3DOrd(myPoly, *mask, maskValue, *f, *fErr, *x, *y, *z): +This is a private routine which will fit a 3-D polynomial to a set of (x, +y, z)-(f) pairs. All non-NULL vectors must be of type PS_TYPE_F64. + +XXX: This routine needs to be written. Currently, this is simply a shell. We +can assume that all vectors have been converted to F64, that (f, x, y, z) are +non-null and F64. fErr may be NULL, but will be F64 is not. + *****************************************************************************/ +psPolynomial3D* VectorFitPolynomial3DOrd( + psPolynomial3D* myPoly, + const psVector* mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z) +{ + // These ASSERTS are redundant. + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); + + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + + psError(PS_ERR_UNKNOWN, true, "3-D Polynomial Fitting is not Implemented.\n"); + return (NULL); +} + +/****************************************************************************** +psVectorFitPolynomial3D(): This routine fits a 3D polynomial of arbitrary +degree (specified in poly) to the data points (x, y, z)-(f) and returns that +polynomial. Types F32 and F64 are supported, however, type F32 is done via +vector conversion only. + *****************************************************************************/ +psPolynomial3D *psVectorFitPolynomial3D( + psPolynomial3D *poly, const psVector *mask, psMaskType maskValue, @@ -2679,9 +2307,171 @@ const psVector *z) { - // Internal pointers for possibly NULL vectors. - psVector *fErr32 = NULL; + // Internal pointers for possibly NULL or mis-typed vectors. + psVector *x64 = NULL; + psVector *y64 = NULL; + psVector *z64 = NULL; + psVector *f64 = NULL; + psVector *fErr64 = NULL; PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + + // + // f + // + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); + if (f->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(f64, f); + } else { + f64 = (psVector *) f; + } + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + + // + // x + // + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + if (x->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(x64, x); + } else { + x64 = (psVector *) x; + } + + // + // y + // + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + if (y->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(y64, y); + } else { + y64 = (psVector *) y; + } + + // + // z + // + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + if (z->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(y64, z); + } else { + z64 = (psVector *) z; + } + + // + // fErr + // + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + if (fErr->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr); + } else { + fErr64 = (psVector *) fErr; + } + } + + if (poly->type == PS_POLYNOMIAL_ORD) { + poly = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64); + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if (z->type.type != PS_TYPE_F64) { + psFree(z64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + return(NULL); + } + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + if (mask != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + psFree(poly); + poly = NULL; + } else { + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if (z->type.type != PS_TYPE_F64) { + psFree(z64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); + } + + + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if (z->type.type != PS_TYPE_F64) { + psFree(z64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + + return(poly); +} + +psPolynomial3D *psVectorClipFitPolynomial3D( + psPolynomial3D *poly, + psStats *stats, + const psVector *mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z) +{ + PS_ASSERT_POLY_NON_NULL(poly, NULL); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_PTR_NON_NULL(stats, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); @@ -2696,31 +2486,83 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); if (fErr != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); } psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); - if (poly->type == PS_POLYNOMIAL_ORD) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D polynomial vector fitting has not been implemented. Returning NULL.\n"); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - } else { - printf("XXX: ERROR: incorrect polynomial type.\n"); - } - if (fErr == NULL) { - psFree(fErr32); - } return(NULL); } -psPolynomial4D *psVectorFitPolynomial4D_NEW( - psPolynomial1D *poly, + +/****************************************************************************** + ****************************************************************************** + 4-D Vector Fitting Code. + ****************************************************************************** + *****************************************************************************/ +/****************************************************************************** +VectorFitPolynomial4DOrd(myPoly, *mask, maskValue, *f, *fErr, *x, *y, *z, *t): +This is a private routine which will fit a 4-D polynomial to a set of (x, +y, z, t)-(f) pairs. All non-NULL vectors must be of type PS_TYPE_F64. + +XXX: This routine needs to be written. Currently, this is simply a shell. We +can assume that all vectors have been converted to F64, that (f, x, y, z) are +non-null and F64. fErr may be NULL, but will be F64 is not. + *****************************************************************************/ +psPolynomial4D* VectorFitPolynomial4DOrd( + psPolynomial4D* myPoly, + const psVector* mask, + psMaskType maskValue, + const psVector *f, + const psVector *fErr, + const psVector *x, + const psVector *y, + const psVector *z, + const psVector *t) +{ + // These ASSERTS are redundant. + PS_ASSERT_POLY_NON_NULL(myPoly, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); + } + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTOR_NON_NULL(z, NULL); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTOR_NON_NULL(t, NULL); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + } + + psError(PS_ERR_UNKNOWN, true, "4-D Polynomial Fitting is not Implemented.\n"); + return (NULL); +} + +/****************************************************************************** +psVectorFitPolynomial4D(): This routine fits a 4D polynomial of arbitrary +degree (specified in poly) to the data points (x, y, z, t)-(f) and returns +that polynomial. Types F32 and F64 are supported, however, type F32 is done +via vector conversion only. + *****************************************************************************/ +psPolynomial4D *psVectorFitPolynomial4D( + psPolynomial4D *poly, const psVector *mask, psMaskType maskValue, @@ -2732,183 +2574,184 @@ const psVector *t) { - // Internal pointers for possibly NULL vectors. - psVector *fErr32 = NULL; + // Internal pointers for possibly NULL or mis-typed vectors. + psVector *x64 = NULL; + psVector *y64 = NULL; + psVector *z64 = NULL; + psVector *t64 = NULL; + psVector *f64 = NULL; + psVector *fErr64 = NULL; PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + + // + // f + // PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); + if (f->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(f64, f); + } else { + f64 = (psVector *) f; + } if (mask != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); } + + // + // x + // PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); + if (x->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(x64, x); + } else { + x64 = (psVector *) x; + } + + // + // y + // PS_ASSERT_VECTOR_NON_NULL(y, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + if (y->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(y64, y); + } else { + y64 = (psVector *) y; + } + + // + // z + // PS_ASSERT_VECTOR_NON_NULL(z, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + if (z->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(y64, z); + } else { + z64 = (psVector *) z; + } + + // + // t + // PS_ASSERT_VECTOR_NON_NULL(t, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F32, NULL); - if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + if (t->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(y64, t); } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D polynomial vector fitting has not been implemented. Returning NULL.\n"); - } else if (poly->type == PS_POLYNOMIAL_CHEB) { - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - } else { - printf("XXX: ERROR: incorrect polynomial type.\n"); - } - if (fErr == NULL) { - psFree(fErr32); - } - return(NULL); -} - - -psPolynomial1D *psVectorClipFitPolynomial1D_NEW( - psPolynomial1D *poly, - psStats *stats, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x) -{ - // Internal pointers for possibly NULL vectors. - psVector *x32 = NULL; - psVector *fErr32 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); - } + t64 = (psVector *) t; + } + + // + // fErr + // if (fErr != NULL) { PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); + if (fErr->type.type != PS_TYPE_F64) { + PS_VECTOR_GEN_F64_FROM_F32(fErr64, fErr); + } else { + fErr64 = (psVector *) fErr; + } + } + + if (poly->type == PS_POLYNOMIAL_ORD) { + poly = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64); + if (poly == NULL) { + psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if (z->type.type != PS_TYPE_F64) { + psFree(z64); + } + + if (t->type.type != PS_TYPE_F64) { + psFree(t64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + return(NULL); + } + } else if (poly->type == PS_POLYNOMIAL_CHEB) { + if (mask != NULL) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + } + psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + psFree(poly); + poly = NULL; } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); - } - - psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); + // Free psVectors that were created for NULL arguments. + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if (z->type.type != PS_TYPE_F64) { + psFree(z64); + } + + if (t->type.type != PS_TYPE_F64) { + psFree(t64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); + } + + // Free psVectors that were created for NULL arguments. - if (x == NULL) { - psFree(x32); - } - if (fErr == NULL) { - psFree(fErr32); - } - return(NULL); -} - -psPolynomial2D *psVectorClipFitPolynomial2D_NEW( - psPolynomial1D *poly, - psStats *stats, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x, - const psVector *y) -{ - // Internal pointers for possibly NULL vectors. - psVector *fErr32 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); - } - - psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); - if (fErr == NULL) { - psFree(fErr32); - } - return(NULL); -} - -psPolynomial3D *psVectorClipFitPolynomial3D_NEW( - psPolynomial1D *poly, - psStats *stats, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x, - const psVector *y, - const psVector *z) -{ - // Internal pointers for possibly NULL vectors. - psVector *fErr32 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); - if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); - } - - psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); - if (fErr == NULL) { - psFree(fErr32); - } - return(NULL); -} - -psPolynomial4D *psVectorClipFitPolynomial4D_NEW( - psPolynomial1D *poly, + if (f->type.type != PS_TYPE_F64) { + psFree(f64); + } + + if (x->type.type != PS_TYPE_F64) { + psFree(x64); + } + + if (y->type.type != PS_TYPE_F64) { + psFree(y64); + } + + if (z->type.type != PS_TYPE_F64) { + psFree(z64); + } + + if (t->type.type != PS_TYPE_F64) { + psFree(t64); + } + + if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { + psFree(fErr64); + } + + return(poly); +} + + +psPolynomial4D *psVectorClipFitPolynomial4D( + psPolynomial4D *poly, psStats *stats, const psVector *mask, @@ -2921,9 +2764,7 @@ const psVector *t) { - // Internal pointers for possibly NULL vectors. - psVector *fErr32 = NULL; - PS_ASSERT_POLY_NON_NULL(poly, NULL); PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); + PS_ASSERT_PTR_NON_NULL(stats, NULL); PS_ASSERT_VECTOR_NON_NULL(f, NULL); PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); @@ -2938,7 +2779,7 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(f, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, f, NULL); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F32, NULL); PS_ASSERT_VECTOR_NON_NULL(t, NULL); PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); @@ -2947,13 +2788,8 @@ PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, mask, NULL); PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F32, NULL); - } else { - fErr32 = psVectorAlloc(f->n, PS_TYPE_F32); - PS_VECTOR_SET_F32(fErr32, 1.0); } psLogMsg(__func__, PS_LOG_WARN, "WARNING: This function has not been implemented. Returning NULL.\n"); - if (fErr == NULL) { - psFree(fErr32); - } + return(NULL); } Index: /trunk/psLib/src/math/psMinimize.h =================================================================== --- /trunk/psLib/src/math/psMinimize.h (revision 4990) +++ /trunk/psLib/src/math/psMinimize.h (revision 4991) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.56 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-07 23:46:04 $ + * @version $Revision: 1.57 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -85,16 +85,6 @@ * @return psPolynomial1D* polynomial fit */ -psPolynomial1D* psVectorFitPolynomial1D( - psPolynomial1D* poly, ///< Polynomial to fit - const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y, ///< Coordinates - const psVector* yErr ///< Errors in coordinates, or NULL -); - - - - - -psPolynomial1D *psVectorFitPolynomial1D_NEW( + +psPolynomial1D *psVectorFitPolynomial1D( psPolynomial1D *poly, const psVector *mask, @@ -105,6 +95,6 @@ ); -psPolynomial2D *psVectorFitPolynomial2D_NEW( - psPolynomial1D *poly, +psPolynomial2D *psVectorFitPolynomial2D( + psPolynomial2D *poly, const psVector *mask, psMaskType maskValue, @@ -115,6 +105,6 @@ ); -psPolynomial3D *psVectorFitPolynomial3D_NEW( - psPolynomial1D *poly, +psPolynomial3D *psVectorFitPolynomial3D( + psPolynomial3D *poly, const psVector *mask, psMaskType maskValue, @@ -126,6 +116,6 @@ ); -psPolynomial4D *psVectorFitPolynomial4D_NEW( - psPolynomial1D *poly, +psPolynomial4D *psVectorFitPolynomial4D( + psPolynomial4D *poly, const psVector *mask, psMaskType maskValue, @@ -139,5 +129,5 @@ -psPolynomial1D *psVectorClipFitPolynomial1D_NEW( +psPolynomial1D *psVectorClipFitPolynomial1D( psPolynomial1D *poly, psStats *stats, @@ -149,6 +139,6 @@ ); -psPolynomial2D *psVectorClipFitPolynomial2D_NEW( - psPolynomial1D *poly, +psPolynomial2D *psVectorClipFitPolynomial2D( + psPolynomial2D *poly, psStats *stats, const psVector *mask, @@ -160,6 +150,6 @@ ); -psPolynomial3D *psVectorClipFitPolynomial3D_NEW( - psPolynomial1D *poly, +psPolynomial3D *psVectorClipFitPolynomial3D( + psPolynomial3D *poly, psStats *stats, const psVector *mask, @@ -172,6 +162,6 @@ ); -psPolynomial4D *psVectorClipFitPolynomial4D_NEW( - psPolynomial1D *poly, +psPolynomial4D *psVectorClipFitPolynomial4D( + psPolynomial4D *poly, psStats *stats, const psVector *mask, @@ -183,26 +173,4 @@ const psVector *z, const psVector *t -); - - - - -/** Derive a one-dimensional spline fit. - * - * Given a psSpline1D data structure and a set of x,y vectors, this routine - * generates the linear splines which satisfy those data points. - * - * @return psSpline1D*: the calculated one-dimensional splines - */ -psSpline1D *psVectorFitSpline1D( - psSpline1D *mySpline, ///< The spline which will be generated. - const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y, ///< Coordinates - const psVector* yErr ///< Errors in coordinates, or NULL -); -psSpline1D *psVectorFitSpline1DNEW( - const psVector* x, ///< Ordinates (or NULL to just use the indices) - const psVector* y, ///< Coordinates - int nKnots ); @@ -233,4 +201,14 @@ const psVector *yErr, ///< Errors in the measurement coordinates psMinimizeLMChi2Func func ///< Specified function +); + +bool psMinimizeGaussNewtonDelta ( + psVector *delta, + const psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func ); Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 4990) +++ /trunk/psLib/src/math/psPolynomial.c (revision 4991) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.119 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-08 00:14:31 $ +* @version $Revision: 1.120 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -564,5 +564,5 @@ of numbers with the specified mean and sigma. -XXX: It's possible to have a different seed everutime. However, for now, +XXX: It's possible to have a different seed everytime. However, for now, for testability, we use a common seed. *****************************************************************************/ Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 4990) +++ /trunk/psLib/src/math/psSpline.c (revision 4991) @@ -7,6 +7,6 @@ * splines. * -* @version $Revision: 1.122 $ $Name: not supported by cvs2svn $ -* @date $Date: 2005-09-08 00:14:32 $ +* @version $Revision: 1.123 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-11 22:18:40 $ * * @@ -207,4 +207,502 @@ } + +/***************************************************************************** +CalculateSecondDerivs(): Given a set of x/y vectors corresponding to a +tabulated function at n points, this routine calculates the second +derivatives of the interpolating cubic splines at those n points. + +The first and second derivatives at the endpoints, undefined in the SDR, are +here defined to be 0.0. They can be modified via ypo and yp1. + +This routine assumes that vectors x and y are of the appropriate types/sizes +(F32). + +XXX: This algorithm is derived from the Numerical Recipes. +XXX: use recycled vectors for internal data. +XXX: do an F64 version? + *****************************************************************************/ +static psF32 *calculateSecondDerivs(const psVector* x, ///< Ordinates (or NULL to just use the indices) + const psVector* y) ///< Coordinates +{ + psTrace(".psLib.dataManip.calculateSecondDerivs", 4, + "---- calculateSecondDerivs() begin ----\n"); + + psS32 i; + psS32 k; + psF32 sig; + psF32 p; + psS32 n = y->n; + psF32 *u = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *derivs2 = (psF32 *) psAlloc(n * sizeof(psF32)); + psF32 *X = (psF32 *) & (x->data.F32[0]); + psF32 *Y = (psF32 *) & (y->data.F32[0]); + psF32 qn; + + // XXX: The second derivatives at the endpoints, undefined in the SDR, + // are set in psConstants.h: PS_LEFT_SPLINE_DERIV, PS_RIGHT_SPLINE_DERIV. + derivs2[0] = -0.5; + u[0]= (3.0/(X[1]-X[0])) * ((Y[1]-Y[0])/(X[1]-X[0]) - PS_LEFT_SPLINE_DERIV); + + for (i=1;i<=(n-2);i++) { + sig = (X[i] - X[i-1]) / (X[i+1] - X[i-1]); + p = sig * derivs2[i-1] + 2.0; + derivs2[i] = (sig - 1.0) / p; + u[i] = ((Y[i+1] - Y[i])/(X[i+1]-X[i])) - ((Y[i]-Y[i-1])/(X[i]-X[i-1])); + u[i] = ((6.0 * u[i] / (X[i+1] - X[i-1])) - (sig * u[i-1])) / p; + + psTrace(".psLib.dataManip.calculateSecondDerivs", 6, + "X[%d] is %f\n", i, X[i]); + psTrace(".psLib.dataManip.calculateSecondDerivs", 6, + "Y[%d] is %f\n", i, Y[i]); + psTrace(".psLib.dataManip.calculateSecondDerivs", 6, + "u[%d] is %f\n", i, u[i]); + } + + qn = 0.5; + u[n-1] = (3.0/(X[n-1]-X[n-2])) * (PS_RIGHT_SPLINE_DERIV - (Y[n-1]-Y[n-2])/(X[n-1]-X[n-2])); + derivs2[n-1] = (u[n-1] - (qn * u[n-2])) / ((qn * derivs2[n-2]) + 1.0); + + for (k=(n-2);k>=0;k--) { + derivs2[k] = derivs2[k] * derivs2[k+1] + u[k]; + + psTrace(".psLib.dataManip.calculateSecondDerivs", 6, + "derivs2[%d] is %f\n", k, derivs2[k]); + } + + psFree(u); + psTrace(".psLib.dataManip.calculateSecondDerivs", 4, + "---- calculateSecondDerivs() end ----\n"); + return(derivs2); +} + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +/***************************************************************************** +psVectorFitSpline1D(): given a psSpline1D data structure and a set of x/y +vectors, this routine generates the linear or cublic splines which satisfy +those data points. + +The formula for calculating the spline polynomials is derived from Numerical +Recipes in C. The basic idea is that the polynomial is + (1) y = (A * y[0]) + + (2) (B * y[1]) + + (3) ((((A*A*A)-A) * mySpline->p_psDeriv2[0]) * H^2)/6.0 + + (4) ((((B*B*B)-B) * mySpline->p_psDeriv2[1]) * H^2)/6.0 +Where: + H = x[1]-x[0] + A = (x[1]-x)/H + B = (x-x[0])/H +The bulk of the code in this routine is the expansion of the above equation +into a polynomial in terms of x, and then saving the coefficients of the +powers of x in the spline polynomials. This gets pretty complicated. + +XXX: usage of yErr is not specified in IfA documentation. + +XXX: Is the x argument redundant? What do we do if the x argument is +supplied, but does not equal the knots specified in mySpline? + +XXX: can psSpline be NULL? + +XXX: reimplement this assuming that mySpline is NULL? + +XXX: What happens if X is NULL, then an index vector is generated for X, but +that index vector lies outside the range vectors in mySpline? + +XXX: Assumes mySpline->knots is psF32. Must add psU32 and psF64. + *****************************************************************************/ +psSpline1D *psVectorFitSpline1D(psSpline1D *mySpline, ///< The spline which will be generated. + const psVector* x, ///< Ordinates (or NULL to just use the indices) + const psVector* y, ///< Coordinates + const psVector* yErr) ///< Errors in coordinates, or NULL +{ + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); + if (mySpline != NULL) { + PS_ASSERT_VECTOR_TYPE(mySpline->knots, PS_TYPE_F32, NULL); + } + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() begin ----\n"); + psS32 numSplines = (y->n)-1; + psF32 tmp; + psF32 H; + psS32 i; + psF32 slope; + psVector *x32 = NULL; + psVector *y32 = NULL; + psVector *yErr32 = NULL; + static psVector *x32Static = NULL; + static psVector *y32Static = NULL; + static psVector *yErr32Static = NULL; + + PS_VECTOR_CONVERT_F64_TO_F32_STATIC(y, y32, y32Static); + + // If yErr==NULL, set all errors equal. + if (yErr == NULL) { + PS_VECTOR_GEN_YERR_STATIC_F32(yErr32Static, y->n); + yErr32 = yErr32Static; + } else { + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(yErr, NULL); + PS_VECTOR_CONVERT_F64_TO_F32_STATIC(yErr, yErr32, yErr32Static); + } + + // If x==NULL, create an x32 vector with x values set to (0:n). + if (x == NULL) { + PS_VECTOR_GEN_X_INDEX_STATIC_F32(x32Static, y->n); + x32 = x32Static; + } else { + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); + PS_VECTOR_CONVERT_F64_TO_F32_STATIC(x, x32, x32Static); + } + PS_ASSERT_VECTORS_SIZE_EQUAL(x32, y32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(yErr32, y32, NULL); + + /* + XXX: + This can not be implemented until SDR states what order spline should be + created. + Should we error if mySpline is not NULL? + Should we error if mySPline is not NULL? + */ + if (mySpline == NULL) { + mySpline = psSpline1DAllocGeneric(x32, 3); + } + PS_ASSERT_PTR_NON_NULL(mySpline, NULL); + PS_ASSERT_INT_NONNEGATIVE(mySpline->n, NULL); + + if (y32->n != (1 + mySpline->n)) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "data size / spline size mismatch (%d %d)\n", + y32->n, mySpline->n); + return(NULL); + } + + // If these are linear splines, which means their polynomials will have + // two coefficients, then we do the simple calculation. + if (2 == (mySpline->spline[0])->n) { + for (i=0;in;i++) { + slope = (y32->data.F32[i+1] - y32->data.F32[i]) / + (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]); + (mySpline->spline[i])->coeff[0] = y32->data.F32[i] - + (slope * mySpline->knots->data.F32[i]); + + (mySpline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, + "---- mySpline %d coeffs are (%f, %f)\n", i, + (mySpline->spline[i])->coeff[0], + (mySpline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, + "---- Exiting psVectorFitSpline1D()()\n"); + return((psSpline1D *) mySpline); + } + + // Check if these are cubic splines (n==4). If not, psError. + if (4 != (mySpline->spline[0])->n) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Don't know how to generate %d-order splines.", + (mySpline->spline[0])->n-1); + return(NULL); + } + + // If we get here, then we know these are cubic splines. We first + // generate the second derivatives at each data point. + mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32); + for (i=0;in;i++) + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); + + // We generate the coefficients of the spline polynomials. I can't + // concisely explain how this code works. See above function comments + // and Numerical Recipes in C. + for (i=0;idata.F32[i+1] - x32->data.F32[i]; + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "x data (%f - %f) (%f)\n", + x32->data.F32[i], + x32->data.F32[i+1], H); + // + // ******** Calculate 0-order term ******** + // + // From (1) + (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H); + // From (2) + ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H); + // From (3) + tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); + tmp-= (x32->data.F32[i+1] / H); + tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + // From (4) + tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); + tmp+= (x32->data.F32[i] / H); + tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + + // + // ******** Calculate 1-order term ******** + // + // From (1) + (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H; + // From (2) + ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H); + // From (3) + tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); + tmp+= (1.0 / H); + tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; + ((mySpline->spline[i])->coeff[1])+= tmp; + // From (4) + tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); + tmp-= (1.0 / H); + tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; + ((mySpline->spline[i])->coeff[1])+= tmp; + + // + // ******** Calculate 2-order term ******** + // + // From (3) + (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H); + // From (4) + ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H)); + + // + // ******** Calculate 3-order term ******** + // + // From (3) + (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); + // From (4) + ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); + + } + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() end ----\n"); + return(mySpline); +} + + + + + + + + + +/***************************************************************************** +psVectorFitSpline1DNEW(): given a psSpline1D data structure and a set of x/y + +xF32 and yF32 are internal psVectors which are used to hold the psF32 versions +of the input data, if necessary. xPtr and yPtr are pointers to either xF32 or +the x argument. All computation is done on xPtr and yPtr. xF32 and yF32 will +simply be psFree() at the end. + +XXX: nKnots makes no sense. This number is always equal to the size of the x +an y vectors. + +XXX: How do we specify the spline order? For now, order=3. + *****************************************************************************/ +#define PS_XXX_SPLINE_ORDER 3 +psSpline1D *psVectorFitSpline1DNEW(const psVector* x, ///< Ordinates (or NULL to just use the indices) + const psVector* y, ///< Coordinates + int nKnots) +{ + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, "---- psVectorFitSpline1DNEW() begin ----\n"); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(y, NULL); + // PS_ASSERT_INT_EQUAL(y->n, nKnots); + + // + // The following code ensures that xPtr points to a psF32 version of the + // ordinate data. + // + psVector *xF32 = NULL; + psVector *xPtr = NULL; + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL); + if (PS_TYPE_F64 == x->type.type) { + xF32 = psVectorAlloc(y->n, PS_TYPE_F32); + for (psS32 i = 0 ; i < x->n ; i++) { + xF32->data.F32[i] = (psF32) x->data.F64[i]; + } + xPtr = xF32; + } else if (PS_TYPE_F32 == x->type.type) { + xPtr = (psVector *) x; + } else { + psError(PS_ERR_UNKNOWN, true, "psVector x is wrong type.\n"); + return(NULL); + } + } else { + // If x==NULL, create an x32 vector with x values set to (0:n). + xF32 = psVectorAlloc(y->n, PS_TYPE_F32); + for (psS32 i = 0 ; i < x->n ; i++) { + xF32->data.F32[i] = (psF32) i; + } + xPtr = xF32; + } + + // + // If y is of type psF64, then create a new vector yF32 and convert the + // y elements. Regardless of y's type, we create a yPtr which will be + // used in the remainder of this function. + // + psVector *yF32 = NULL; + psVector *yPtr = NULL; + if (PS_TYPE_F64 == y->type.type) { + yF32 = psVectorAlloc(y->n, PS_TYPE_F32); + for (psS32 i = 0 ; i < y->n ; i++) { + yF32->data.F32[i] = (psF32) y->data.F64[i]; + } + yPtr = yF32; + } else { + yPtr = (psVector *) y; + } + + psSpline1D *mySpline = psSpline1DAllocGeneric(xPtr, PS_XXX_SPLINE_ORDER); + + psS32 numSplines = nKnots - 1; + psF32 tmp; + psF32 H; + psS32 i; + psF32 slope; + // XXX: get rid of x32 and y32 (this is from old code) + psVector *x32 = xPtr; + psVector *y32 = yPtr; + + // If these are linear splines, which means their polynomials will have + // two coefficients, then we do the simple calculation. + if (1 == PS_XXX_SPLINE_ORDER) { + for (i=0;in;i++) { + slope = (y32->data.F32[i+1] - y32->data.F32[i]) / + (mySpline->knots->data.F32[i+1] - mySpline->knots->data.F32[i]); + (mySpline->spline[i])->coeff[0] = y32->data.F32[i] - + (slope * mySpline->knots->data.F32[i]); + + (mySpline->spline[i])->coeff[1] = slope; + psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, + "---- mySpline %d coeffs are (%f, %f)\n", i, + (mySpline->spline[i])->coeff[0], + (mySpline->spline[i])->coeff[1]); + } + psTrace(".psLib.dataManip.psMinimize.psVectorFitSpline1D", 4, + "---- Exiting psVectorFitSpline1D()()\n"); + return((psSpline1D *) mySpline); + } + + // + // Check if these are cubic splines (n==4). If not, psError. + // + if (3 != PS_XXX_SPLINE_ORDER) { + psError(PS_ERR_BAD_PARAMETER_SIZE, true, + "Don't know how to generate %d-order splines.", PS_XXX_SPLINE_ORDER); + return(NULL); + } + + // + // If we get here, then we know these are cubic splines. We first + // generate the second derivatives at each data point. + // + mySpline->p_psDeriv2 = calculateSecondDerivs(x32, y32); + for (i=0;in;i++) + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "Second deriv[%d] is %f\n", i, mySpline->p_psDeriv2[i]); + + // + // We generate the coefficients of the spline polynomials. I can't + // concisely explain how this code works. See above function comments + // and Numerical Recipes in C. + // + for (i=0;idata.F32[i+1] - x32->data.F32[i]; + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "x data (%f - %f) (%f)\n", + x32->data.F32[i], + x32->data.F32[i+1], H); + // + // ******** Calculate 0-order term ******** + // + // From (1) + (mySpline->spline[i])->coeff[0] = (y32->data.F32[i] * x32->data.F32[i+1]/H); + // From (2) + ((mySpline->spline[i])->coeff[0])-= ((y32->data.F32[i+1] * x32->data.F32[i])/H); + // From (3) + tmp = (x32->data.F32[i+1] * x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); + tmp-= (x32->data.F32[i+1] / H); + tmp*= (mySpline->p_psDeriv2)[i] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + // From (4) + tmp = -(x32->data.F32[i] * x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); + tmp+= (x32->data.F32[i] / H); + tmp*= (mySpline->p_psDeriv2)[i+1] * H * H / 6.0; + ((mySpline->spline[i])->coeff[0])+= tmp; + + // + // ******** Calculate 1-order term ******** + // + // From (1) + (mySpline->spline[i])->coeff[1] = -(y32->data.F32[i]) / H; + // From (2) + ((mySpline->spline[i])->coeff[1])+= (y32->data.F32[i+1] / H); + // From (3) + tmp = -3.0 * (x32->data.F32[i+1] * x32->data.F32[i+1]) / (H * H * H); + tmp+= (1.0 / H); + tmp*= ((mySpline->p_psDeriv2)[i]) * H * H / 6.0; + ((mySpline->spline[i])->coeff[1])+= tmp; + // From (4) + tmp = 3.0 * (x32->data.F32[i] * x32->data.F32[i]) / (H * H * H); + tmp-= (1.0 / H); + tmp*= ((mySpline->p_psDeriv2)[i+1]) * H * H / 6.0; + ((mySpline->spline[i])->coeff[1])+= tmp; + + // + // ******** Calculate 2-order term ******** + // + // From (3) + (mySpline->spline[i])->coeff[2] = ((mySpline->p_psDeriv2)[i]) * 3.0 * x32->data.F32[i+1] / (6.0 * H); + // From (4) + ((mySpline->spline[i])->coeff[2])-= (((mySpline->p_psDeriv2)[i+1]) * 3.0 * x32->data.F32[i] / (6.0 * H)); + + // + // ******** Calculate 3-order term ******** + // + // From (3) + (mySpline->spline[i])->coeff[3] = -((mySpline->p_psDeriv2)[i]) / (6.0 * H); + // From (4) + ((mySpline->spline[i])->coeff[3])+= ((mySpline->p_psDeriv2)[i+1]) / (6.0 * H); + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[0] is %f\n", i, (mySpline->spline[i])->coeff[0]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[1] is %f\n", i, (mySpline->spline[i])->coeff[1]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[2] is %f\n", i, (mySpline->spline[i])->coeff[2]); + psTrace(".psLib.dataManip.psVectorFitSpline1D", 6, + "(mySpline->spline[%d])->coeff[3] is %f\n", i, (mySpline->spline[i])->coeff[3]); + + } + + psFree(xF32); + psFree(yF32); + + psTrace(".psLib.dataManip.psVectorFitSpline1D", 4, + "---- psVectorFitSpline1D() end ----\n"); + + return(mySpline); +} + + + + + /*****************************************************************************/ /* FUNCTION IMPLEMENTATION - PUBLIC */ @@ -283,5 +781,5 @@ } } else { - printf("XXX: Generate an error here.\n"); + psError(PS_ERR_UNKNOWN, true, "psVector in has an incorrect type.\n"); return(NULL); } Index: /trunk/psLib/src/math/psSpline.h =================================================================== --- /trunk/psLib/src/math/psSpline.h (revision 4990) +++ /trunk/psLib/src/math/psSpline.h (revision 4991) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.54 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-08 00:02:48 $ + * @version $Revision: 1.55 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -118,4 +118,24 @@ ); +/** Derive a one-dimensional spline fit. + * + * Given a psSpline1D data structure and a set of x,y vectors, this routine + * generates the linear splines which satisfy those data points. + * + * @return psSpline1D*: the calculated one-dimensional splines + */ +psSpline1D *psVectorFitSpline1D( + psSpline1D *mySpline, ///< The spline which will be generated. + const psVector* x, ///< Ordinates (or NULL to just use the indices) + const psVector* y, ///< Coordinates + const psVector* yErr ///< Errors in coordinates, or NULL +); + +psSpline1D *psVectorFitSpline1DNEW( + const psVector* x, ///< Ordinates (or NULL to just use the indices) + const psVector* y, ///< Coordinates + int nKnots +); + /** \} */ // End of MathGroup Functions Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 4990) +++ /trunk/psLib/src/math/psStats.c (revision 4991) @@ -17,6 +17,6 @@ * * - * @version $Revision: 1.144 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-07 21:40:09 $ + * @version $Revision: 1.145 $ $Name: not supported by cvs2svn $ + * @date $Date: 2005-09-11 22:18:40 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -746,4 +746,7 @@ XXX: Write a general routine which smoothes a psVector. This routine should call that. Is that possible? + +XXX: This is, or will be, prettier than the previous +psVectorSmoothHistGaussian(). However, it is not being used, yet. *****************************************************************************/ psVector *p_psVectorSmoothHistGaussianNEW(psHistogram *histogram, @@ -1549,5 +1552,6 @@ // Determine the coefficients of the polynomial. - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); + // myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); + myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x); if (myPoly == NULL) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1825,5 +1829,5 @@ psPolynomial1D *tmpPoly = psPolynomial1DAlloc(3, PS_POLYNOMIAL_ORD); // XXX: What about the NULL x argument? - tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, y, NULL); + tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, 0, y, NULL, NULL); if (tmpPoly == NULL) { psLogMsg(__func__, PS_LOG_WARN, @@ -1936,5 +1940,5 @@ } } else { - printf("XXX: Generate an error here.\n"); + psError(PS_ERR_UNKNOWN, true, "Unallowable vector type.\n"); return(NULL); } @@ -1950,6 +1954,8 @@ psVector *Fit1DGaussian(psVector *x, psVector*y) { - printf("XXX: Generate an error here.\n"); - printf("XXX: Error: This function was previously part of psStats.c, was removed, was purged from CVS, and now needs to be retrieved from tape.\n"); + psError(PS_ERR_UNKNOWN, true, "This code has not been implemented yet.\n"); + // XXX: This function was previously part of psStats.c, was removed, was + // purged from CVS, and now needs to be retrieved from tape. + return(NULL); } @@ -2168,8 +2174,6 @@ stats->robustLQ = binLo25F32; stats->robustUQ = binHi25F32; - // XXX: No idea how to calculate stats->stdev // PS_BIN_MIDPOINT(robustHistogram, modeBinNum); - // XXX: I think sumNfit == sumN50 here. stats->robustNfit = -1;