Index: trunk/psLib/src/math/psMinimize.c =================================================================== --- trunk/psLib/src/math/psMinimize.c (revision 1907) +++ trunk/psLib/src/math/psMinimize.c (revision 1921) @@ -9,6 +9,6 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.49 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-09-27 23:41:42 $ + * @version $Revision: 1.50 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-09-28 23:27:37 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -656,43 +656,16 @@ -/***************************************************************************** -CreateChebyshevPolys(n): this routine takes as input the required order n, -and returns as output as a pointer to an array of n psPolynomial1D -structures, corresponding to the first n Chebyshev polynomials. - -XXX: The output should be static since the Chebyshev polynomials might be -used frequently and the data structure created here does not contain the -outer coefficients of the Chebyshev polynomials. +/****************************************************************************** +p_psVectorFitPolynomial1DCheb(): 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. + +XXX: must add type F32 (currently F64 only). + +XXX: Use private name? *****************************************************************************/ -static psPolynomial1D **CreateChebyshevPolys(int maxChebyPoly) -{ - psPolynomial1D **chebPolys = NULL; - int i = 0; - int j = 0; - - chebPolys = (psPolynomial1D **) psAlloc(maxChebyPoly * sizeof(psPolynomial1D *)); - for (i = 0; i < maxChebyPoly; i++) { - chebPolys[i] = psPolynomial1DAlloc(i + 1, PS_POLYNOMIAL_ORD); - } - - // Create the Chebyshev polynomials. - // Polynomial i has i-th order. - chebPolys[0]->coeff[0] = 1; - chebPolys[1]->coeff[1] = 1; - for (i = 2; i < maxChebyPoly; i++) { - for (j = 0; j < chebPolys[i - 1]->n; j++) { - chebPolys[i]->coeff[j + 1] = 2 * chebPolys[i - 1]->coeff[j]; - } - for (j = 0; j < chebPolys[i - 2]->n; j++) { - chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j]; - } - } - - return (chebPolys); -} - - - -psPolynomial1D* p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly, +psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly, const psVector* restrict x, const psVector* restrict y, @@ -702,42 +675,74 @@ int k; int n = x->n; - psPolynomial1D **chebPolys = CreateChebyshevPolys(myPoly->n); - ; - float fac; - float sum; - /* - fac = 2.0/((float) n); - for (j=0;jdata.F64[k] * - cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n)); - // sum+= y->data.F64[k] * - // cos(M_PI * ((float) j) * (-0.5 + ((float) k)) / ((float) n)) * - // cos(M_PI * (-0.5 + ((float) k)) / ((float) n)); - - } - myPoly->coeff[j] = fac * sum; - - } - return(myPoly); - */ - + psVector *f = psVectorAlloc(n, PS_TYPE_F64); + double fac; + double sum; + // XXX: Use static memory here. + psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); + psScalar *fScalar; + + + // XXX: These assignments appear too simple to warrant code and + // variable declarations. I retain them here to maintain coherence + // with the NR code. + double min = -1.0; + double max = 1.0; + double bma = 0.5 * (max-min); // 1 + double bpa = 0.5 * (max+min); // 0 + + // XXX: Eliminate this later by generating a F64 version of the + // LaGrange interpolation routines. + PS_VECTOR_F64_TO_F32(x, x32); + PS_VECTOR_F64_TO_F32(y, y32); + + // 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 (int i=0;idata.F32 = (float) X; + + // We interpolate against are tabluated x,y vectors to determine the + // function value at X. + fScalar = p_psVectorInterpolate((psVector *) x32, (psVector *) y32, + 3, tmpScalar); + f->data.F64[i] = (double) fScalar->data.F32; + + psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 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/((float) n); - for (j=0;jn;j++) { + for (j=0;jdata.F64[k] * - psPolynomial1DEval((float) x->data.F64[k], chebPolys[j])); + for (k=0;kdata.F64[k] * + cos(M_PI * ((float) j) * (0.5 + ((float) k)) / ((float) n)); } myPoly->coeff[j] = fac * sum; } + + // XXX: Must free memory. + psFree(x32); + psFree(y32); + psFree(tmpScalar); + psFree(fScalar); 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, as well as the error for each data point (yErr). +p_psVectorFitPolynomial1DOrd(): This routine will fit an ordinary +polynomial of degree myPoly to the data points (x, y) and return the +coefficients of that polynomial. XXX: yErr is currently ignored. @@ -745,23 +750,11 @@ XXX: must add type F32 (currently F64 only). -XXX: Must do: if x==NULL, use the index vector (???). - -XXX: Must do: if yErr==NULL, set all errors equal. - -XXX: type is currently ignored. Must implement this for Chebyshev -polynomials. +XXX: Use private name? *****************************************************************************/ -psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, - const psVector* restrict x, - const psVector* restrict y, - const psVector* restrict yErr) +psPolynomial1D* p_psVectorFitPolynomial1DOrd(psPolynomial1D* myPoly, + const psVector* restrict x, + const psVector* restrict y, + const psVector* restrict yErr) { - // XXX: Create a p_psVectorFitPolynomial1DOrd() function. Here, we - // should only be calling that or the Cheby function. - - if (myPoly->type == PS_POLYNOMIAL_CHEB) { - return(p_psVectorFitPolynomial1DCheby(myPoly, x, y, yErr)); - } - int polyOrder = myPoly->n; psImage* A = NULL; @@ -783,16 +776,8 @@ // } + // XXX: Some of these are redundant. PS_CHECK_NULL_1DPOLY(myPoly); PS_CHECK_NULL_VECTOR(y); PS_CHECK_EMPTY_VECTOR(y); - - // XXX: Verify that this is the correct action. - if (x == NULL) { - x = psVectorAlloc(y->n, PS_TYPE_F32); - for (i=0;in;i++) { - x->data.F32[i] = (float) i; - } - } - PS_CHECK_EMPTY_VECTOR(x); PS_CHECK_NULL_VECTOR(yErr); @@ -863,4 +848,113 @@ 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: yErr is currently ignored. + +XXX: must add type F32 (currently F64 only). + *****************************************************************************/ +psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly, + const psVector* restrict x, + const psVector* restrict y, + const psVector* restrict yErr) +{ + bool mustFreeMyYErr = false; + bool mustFreemyX = false; + int i; + psPolynomial1D *tmpPoly; + psVector *myX = NULL; + psVector *myYErr = NULL; + + PS_CHECK_NULL_1DPOLY(myPoly); + PS_CHECK_NULL_VECTOR(y); + PS_CHECK_EMPTY_VECTOR(y); + + // If yErr==NULL, set all errors equal. + if (yErr == NULL) { + myYErr = psVectorAlloc(y->n, y->type.type); + mustFreeMyYErr = true; + + if (y->type.type == PS_TYPE_F32) { + for (i=0;in;i++) { + myYErr->data.F32[i] = 1.0; + } + } else if (y->type.type == PS_TYPE_F64) { + for (i=0;in;i++) { + myYErr->data.F64[i] = 1.0; + } + } + } else { + myYErr = (psVector *) yErr; + } + + // If x==NULL, create an myX vector with x values set to (0:n), and if + // this is a CHebyshev polynomial, we must scale to (-1:1). + + // XXX: Verify that this is the correct action. + if (x == NULL) { + myX = psVectorAlloc(y->n, y->type.type); + mustFreemyX = true; + + if (y->type.type == PS_TYPE_F32) { + if (myPoly->type == PS_POLYNOMIAL_ORD) { + for (i=0;in;i++) { + myX->data.F32[i] = (float) i; + } + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + float min = 0.0; + float max = (float) (y->n - 1); + + for (i=0;in;i++) { + myX->data.F32[i] = (((float) i) - 0.5 * (min + max)) / + (0.5 * (max - min)); + } + } + } else if (y->type.type == PS_TYPE_F64) { + if (myPoly->type == PS_POLYNOMIAL_ORD) { + for (i=0;in;i++) { + myX->data.F64[i] = (float) i; + } + } else if (myPoly->type == PS_POLYNOMIAL_CHEB) { + double min = 0.0; + double max = (double) (y->n - 1); + + for (i=0;in;i++) { + myX->data.F64[i] = (((float) i) - 0.5 * (min + max)) / + (0.5 * (max - min)); + } + } + } + } else { + myX = (psVector *) x; + } + + PS_CHECK_VECTOR_SIZE_EQUAL(y, myX); + PS_CHECK_VECTOR_SIZE_EQUAL(y, myYErr); + + // Call the appropriate vector fitting routine. + if (myPoly->type == PS_POLYNOMIAL_CHEB) { + tmpPoly = p_psVectorFitPolynomial1DCheby(myPoly, myX, y, myYErr); + } else if (myPoly->type == PS_POLYNOMIAL_ORD) { + tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, myX, y, myYErr); + } else { + // XXX: psErrorMsg() + return(NULL); + } + + // Free any allocated memory. + if (mustFreeMyYErr == true) { + psFree(myYErr); + } + if (mustFreemyX == true) { + psFree(myX); + } + + return(myPoly); +} +