Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 2267) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 2269) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.83 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-11-02 19:08:33 $ + * @version $Revision: 1.84 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-11-03 03:30:30 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -61,10 +61,9 @@ returned as a psVector sums. -XXX: change name XXX: Use a static vector. *****************************************************************************/ -void p_psBuildSums1D(double x, - psS32 polyOrder, - psVector* sums) +void psBuildSums1D(double x, + psS32 polyOrder, + psVector* sums) { psS32 i = 0; @@ -74,5 +73,4 @@ sums = psVectorAlloc(polyOrder, PS_TYPE_F64); } - if (polyOrder > sums->n) { sums = psVectorRealloc(sums, polyOrder); @@ -109,5 +107,4 @@ psS32 i; psS32 k; - psBool mustFreeX = false; float sig; float p; @@ -150,9 +147,5 @@ } - if (mustFreeX == true) { - psFree(X); - } psFree(u); - psTrace(".psLib.dataManip.CalculateSecondDerivs", 4, "---- CalculateSecondDerivs() end ----\n"); @@ -171,7 +164,7 @@ polynomials which are stored in psSpline1D. -XXX: check types/sizes XXX: This is F32 only *****************************************************************************/ +/* float p_psNRSpline1DEval(psSpline1D *spline, const psVector* restrict x, @@ -187,5 +180,5 @@ PS_VECTOR_CHECK_NULL(y, NAN); PS_VECTOR_CHECK_TYPE(y, PS_TYPE_F32, NAN); - + psS32 n; psS32 klo; @@ -197,7 +190,7 @@ float D; float Y; - + n = spline->n; - klo = p_psVectorBinDisectF32(spline->domains, (spline->n)+1, X); + klo = p_psVectorBinDisect32(spline->domains, (spline->n)+1, X); khi = klo + 1; H = (spline->domains)[khi] - (spline->domains)[klo]; @@ -206,13 +199,13 @@ 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 */ @@ -239,4 +232,5 @@ 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 domains specified in mySpline? @@ -414,5 +408,4 @@ /****************************************************************************** XXX: We assume unnormalized gaussians. -XXX: Currently, yErr is ignored. *****************************************************************************/ psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, @@ -441,6 +434,4 @@ x = ((psVector *) (coords->data[i]))->data.F32[0]; out->data.F32[i] = psGaussian(x, mean, stdev, false); - // psTrace(".psLib.dataManip.psMinimize", 6, - // "out->data.F32[%d]: (x, m, s) is (%f, %f, %f) is %f\n", i, x, mean, stdev, out->data.F32[i]); } @@ -449,5 +440,4 @@ float tmp = (x - mean) * psGaussian(x, mean, stdev, false); deriv->data.F32[i][0] = tmp / (stdev * stdev); - tmp = (x - mean) * (x - mean) * psGaussian(x, mean, stdev, 0); @@ -482,5 +472,5 @@ XXX: This must work for both F32 and F64. F32 is currently implemented. Note: since the LUD routines are only implemented in F64, then we -> will have to convert all F32 input vectors to F64 regardless. So, + will have to convert all F32 input vectors to F64 regardless. So, the F64 port might be. *****************************************************************************/ @@ -730,7 +720,7 @@ /****************************************************************************** -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. +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. @@ -738,5 +728,5 @@ XXX: Use private name? *****************************************************************************/ -psPolynomial1D *p_psVectorFitPolynomial1DCheby(psPolynomial1D* myPoly, +psPolynomial1D *VectorFitPolynomial1DCheby(psPolynomial1D* myPoly, const psVector* restrict x, const psVector* restrict y, @@ -785,5 +775,5 @@ psFree(fScalar); - psTrace(".psLib.dataManip.p_psVectorFitPolynomial1DCheby", 6, + 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]); @@ -809,12 +799,12 @@ /****************************************************************************** -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. +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. *****************************************************************************/ -psPolynomial1D* p_psVectorFitPolynomial1DOrd(psPolynomial1D* myPoly, +psPolynomial1D* VectorFitPolynomial1DOrd(psPolynomial1D* myPoly, const psVector* restrict x, const psVector* restrict y, @@ -833,7 +823,7 @@ psVector* xSums = NULL; - psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4, - "---- psVectorFitPolynomial1D() begin ----\n"); - // printf("psVectorFitPolynomial1D()\n"); + 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]); @@ -865,5 +855,5 @@ if (yErr == NULL) { for (i = 0; i < X->n; i++) { - p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); + psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); for (k = 0; k < polyOrder; k++) { @@ -879,5 +869,5 @@ } else { for (i = 0; i < X->n; i++) { - p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); + psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); for (k = 0; k < polyOrder; k++) { @@ -916,6 +906,6 @@ psFree(xSums); - psTrace(".psLib.dataManip.psVectorFitPolynomial1D", 4, - "---- psVectorFitPolynomial1D() begin ----\n"); + psTrace(".psLib.dataManip.VectorFitPolynomial1DOrd", 4, + "---- VectorFitPolynomial1DOrd() begin ----\n"); return (myPoly); } @@ -962,5 +952,5 @@ PS_VECTOR_GEN_X_INDEX_STATIC_F64(x64Static, y->n); if (myPoly->type == PS_POLYNOMIAL_CHEB) { - p_psNormalizeVector(x64Static); + p_psNormalizeVectorRange(x64Static, -1.0, 1.0); } x64 = x64Static; @@ -974,7 +964,7 @@ // Call the appropriate vector fitting routine. if (myPoly->type == PS_POLYNOMIAL_CHEB) { - p_psVectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64); + VectorFitPolynomial1DCheby(myPoly, x64, y64, yErr64); } else if (myPoly->type == PS_POLYNOMIAL_ORD) { - tmpPoly = p_psVectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64); + tmpPoly = VectorFitPolynomial1DOrd(myPoly, x64, y64, yErr64); } else { psError(__func__, "unknown polynomial type.\n");