Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 2227) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 2228) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-28 00:22:53 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -208,4 +208,10 @@ p_psMemSetPersistent(NAME->data.V, true); \ +#define PS_VECTOR_DECLARE_ALLOC_STATIC(NAME, SIZE, TYPE) \ +static psVector *NAME = NULL; \ +if (NAME == NULL) { \ + NAME = psVectorAlloc(SIZE, TYPE); \ + p_psMemSetPersistent(NAME, true); \ +} \ Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 2227) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 2228) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.72 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:43 $ + * @version $Revision: 1.73 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -408,4 +408,54 @@ } + +/****************************************************************************** +XXX: We assume unnormalized gaussians. +XXX: Currently, yErr is ignored. + *****************************************************************************/ +psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, + const psVector *params, + const psArray *coords) +{ + psTrace(".psLib.dataManip.psMinimize", 4, + "---- psMinimizeLMChi2Gauss1D() begin ----\n"); + + PS_PTR_CHECK_NULL(coords, NULL); + PS_PTR_CHECK_NULL(params, NULL); + float x; + int i; + float mean = params->data.F32[0]; + float stdev = params->data.F32[1]; + psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); + + psTrace(".psLib.dataManip.psMinimize", 6, + "(mean, stdev) is (%f, %f)\n", mean, stdev); + + if (deriv == NULL) { + deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); + } + + for (i=0;in;i++) { + 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]); + } + + for (i=0;in;i++) { + x = ((psVector *) (coords->data[i]))->data.F32[0]; + 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); + deriv->data.F32[i][1] = tmp / (stdev * stdev * stdev); + } + + psTrace(".psLib.dataManip.psMinimize", 4, + "---- psMinimizeLMChi2Gauss1D() end ----\n"); + return(out); +} + + /****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates @@ -429,5 +479,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. *****************************************************************************/ @@ -457,5 +507,5 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2() begin ----\n"); psS32 numData = y->n; @@ -490,8 +540,9 @@ float newChi2 = 0.0; float lamda = 0.00005; - - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + lamda = 0.05; + + psTrace(".psLib.dataManip.psMinimize", 6, "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "min->tol is %f\n", min->tol); @@ -504,7 +555,7 @@ while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "------------------------------------------------------\n"); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); @@ -515,12 +566,16 @@ currValueVec = func(deriv, params, x); for (n=0;ndata.F32[n] * currValueVec->data.F32[n]); + currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) * + (currValueVec->data.F32[n] - y->data.F32[n]); + psTrace(".psLib.dataManip.psMinimize", 6, + "data[%d], chi2 calculation+= (%f * %f)\n", n, + currValueVec->data.F32[n], y->data.F32[n]); } for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "Current chi-squared is (%f)\n", currChi2); @@ -548,5 +603,5 @@ // XXX: multiple by -1 here? (beta->data.F64[p])*= -1.0; - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); } @@ -580,9 +635,9 @@ } for (j=0;jdata.F64[j][k]); - } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n"); + psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]); + } + psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n"); } @@ -597,5 +652,5 @@ // for (i=0;idata.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]); if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { @@ -607,8 +662,8 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "Calling func() with new parameters:\n"); for (i=0;idata.F32[%d] is %f.\n", i, newParams->data.F32[i]); } @@ -621,9 +676,11 @@ newValueVec = func(deriv, newParams, x); for (n=0;ndata.F32[n] * newValueVec->data.F32[n]); + newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) * + (newValueVec->data.F32[n] - y->data.F32[n]); + } psFree(newValueVec); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2); @@ -640,8 +697,10 @@ } lamda*= 0.1; + psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n"); } else { lamda*= 10.0; - } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n"); + } + psTrace(".psLib.dataManip.psMinimize", 4, "lamda is %f\n", lamda); min->iter++; @@ -662,5 +721,5 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2() end (false) ----\n"); return(false); @@ -942,4 +1001,14 @@ } +#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ +for (int i=0;in;i++) { \ + if (BASEMASK->data.U8[i] == 0) { \ + OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ + } else { \ + OUT->data.F32[i] = BASE->data.F32[i]; \ + } \ +} \ + + /****************************************************************************** p_psDetermineBracket(): This routine takes as input an arbitrary function, @@ -947,18 +1016,26 @@ function produces as output a bracket [a, b, c] such that - f(param + b * line) is less than f(param + a * line) and - f(param + c * line). - -Algorithm: XXX completely ad hoc: start with the user-supplied starting +f(param + b * line) is less than f(param + a * line) and +f(param + c * line). + +Algorithm: +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. + +XXX: +output bracket vector should be an input as well. +*****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, psVector *line, @@ -966,4 +1043,155 @@ const psArray *coords, psMinimizePowellFunc func) +{ + float a = 0.0; + float b = 0.0; + float c = 0.0; + float fa = 0.0; + float fb = 0.0; + float fc = 0.0; + psS32 iter = 100; + float aDir = 0.0; + float cDir = 0.0; + float new_aDir = 0.0; + float new_cDir = 0.0; + psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); + float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; + psVector *tmp = NULL; + psS32 i = 0; + psS32 null = 0; + + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() begin ----\n"); + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + if (line->data.F32[i] >= FLT_EPSILON) { + null = 1; + } + } + } + + if (null == 0) { + psTrace(".psLib.dataManip.p_psDetermineBracket", 2, + "p_psDetermineBracket() called with zero line vector.\n"); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end (NULL) ----\n"); + psFree(bracket); + return(NULL); + } + + tmp = psVectorAlloc(params->n, PS_TYPE_F32); + + b = 0; + a = -stepSize; + c = stepSize; + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a); + fa = func(tmp, coords); + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b); + fb = func(tmp, coords); + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c); + fc = func(tmp, coords); + + if (fa < fb) { + aDir = -1; + } else { + aDir = 1; + } + + if (fc < fb) { + cDir = -1; + } else { + cDir = 1; + } + + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); + + while (iter > 0) { + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "psDetermineBracket(): iteration %d\n", iter); + if ((fb < fa) && (fb < fc)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + stepSize*= (1.0 + stepSize); + a =- stepSize; + c =+ stepSize; + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]); + } else { + tmp->data.F32[i] = params->data.F32[i]; + } + } + fa = func(tmp, coords); + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]); + } else { + tmp->data.F32[i] = params->data.F32[i]; + } + } + fc = func(tmp, coords); + + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc); + + if (fa < fb) { + new_aDir = -1; + } else { + new_aDir = 1; + } + + if (fc < fb) { + new_cDir = -1; + } else { + new_cDir = 1; + } + if ((new_aDir == 1) && (aDir == -1)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + + if ((new_cDir == 1) && (cDir == -1)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + aDir = new_aDir; + cDir = new_cDir; + iter--; + } + psFree(tmp); + psFree(bracket); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end (NULL) ----\n"); + return(NULL); +} + +psVector *p_psDetermineBracket2(psVector *params, + psVector *line, + const psVector *paramMask, + const psArray *coords, + psMinimizePowellFunc func) { float a = 0.0; @@ -1365,4 +1593,5 @@ float fqp = 0.0; float diff = 0.0; + int iterationNumber = 0; psTrace(".psLib.dataManip.psMinimizePowell", 4, @@ -1407,9 +1636,8 @@ } - min->iter = 0; - while (min->iter < min->maxIter) { - min->iter++; + while (iterationNumber < min->maxIter) { + iterationNumber++; psTrace(".psLib.dataManip.psMinimizePowell", 6, - "psMinimizePowell() iteration %d\n", min->iter); + "psMinimizePowell() iteration %d\n", iterationNumber); // 3: For each dimension in params, move Q only in the vector v[i] to @@ -1418,4 +1646,7 @@ baseFuncVal = func(Q, coords); currFuncVal = baseFuncVal; + psTrace(".psLib.dataManip.psMinimizePowell", 6, + "Current function value is %f\n", currFuncVal); + diff = 0.0; biggestDiff = 0; @@ -1423,5 +1654,5 @@ for (i=0;idata.U8[i] == 0) { - dummyMin.maxIter = 100; + dummyMin.maxIter = 20; dummyMin.tol = 0.01; mul = p_psLineMin(&dummyMin, Q, v[i], myParamMask, coords, func); @@ -1433,5 +1664,6 @@ } } - + psTrace(".psLib.dataManip.psMinimizePowell", 6, + "New function value is %f\n", currFuncVal); // 4: Set the vector u = Q - P for (i=0;iiter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (true) ----\n"); @@ -1513,9 +1745,10 @@ psFree(pQP); psFree(Q); + if (paramMask == NULL) { + psFree(myParamMask); + } + min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (true) ----\n"); - if (paramMask == NULL) { - psFree(myParamMask); - } return(true); } @@ -1532,4 +1765,5 @@ psFree(pQP); psFree(Q); + min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (false) ----\n"); @@ -1549,4 +1783,12 @@ const psArray *coords) { + psTrace(".psLib.dataManip.myPowellChi2Func", 4, + "---- myPowellChi2Func() begin ----\n"); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_EMPTY(params, NAN); + PS_VECTOR_CHECK_NULL(myValue, NAN); + PS_VECTOR_CHECK_EMPTY(myValue, NAN); + PS_PTR_CHECK_NULL(coords, NAN); + float chi2 = 0.0; float d; @@ -1555,9 +1797,18 @@ tmp = Chi2PowellFunc(params, coords); - for (i=0;in;i++) { - d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i]; - chi2+= d * d; + if (myError == NULL) { + for (i=0;in;i++) { + d = (tmp->data.F32[i] - myValue->data.F32[i]); + chi2+= d * d; + } + } else { + for (i=0;in;i++) { + d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i]; + chi2+= d * d; + } } psFree(tmp); + psTrace(".psLib.dataManip.myPowellChi2Func", 4, + "---- myPowellChi2Func() end (%f) ----\n", chi2); return(chi2); } @@ -1584,4 +1835,5 @@ myValue = (psVector *) value; myError = (psVector *) error; + Chi2PowellFunc = func; Index: /trunk/psLib/src/dataManip/psMinimize.h =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.h (revision 2227) +++ /trunk/psLib/src/dataManip/psMinimize.h (revision 2228) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:43 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -81,41 +81,4 @@ const psArray *x); -/* -psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, - const psVector *params, - const psArray *coords) -{ - psVector *x = (psVector *) coords->data[0]; - psVector *y = (psVector *) coords->data[1]; - psVector *yErr = NULL; - int i; - int j; - float normalization = params->data.F32[0]; - float center = params->data.F32[1]; - float stdev = params->data.F32[2]; - psVector *out = psVectorAlloc(x->n, PS_TYPE_F32); - - - if (deriv == NULL) { - deriv = psImageAlloc(params->n, x->n, PS_TYPE_F32); - } - - if (coords->n == 2) { - yErr = (psVector *) coords->data[2]; - } - - for (i=0;in;i++) { - out->data.F32[i] = psGaussian(x->data.F32[i], center, stdev, 1); - } - - for (i=0;in;i++) { - for (j=0;jn;j++) { - deriv->data.F32[i][j] = the derivative of the gaussian w.r.t. the j-th parameter at the coord x[i] - } - } - return(out); -} -*/ - psBool psMinimizeLMChi2(psMinimization *min, psImage *covar, @@ -137,4 +100,8 @@ psMinimizePowellFunc func); +psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, + const psVector *params, + const psArray *coords); + typedef psVector *(*psMinimizeChi2PowellFunc)(const psVector *params, Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 2227) +++ /trunk/psLib/src/dataManip/psStats.c (revision 2228) @@ -9,7 +9,7 @@ * @author GLG, MHPCC * - * @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-28 00:46:04 $ -n * + * @version $Revision: 1.80 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ + * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -123,22 +123,4 @@ default: return false; - } -} - -/***************************************************************************** -p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal -function that simply prints a vector to STDOUT. Used primarily for -debugging. - *****************************************************************************/ - -void p_psVectorPrint(psVector* myVector, psVector* maskVector, psU32 maskVal, psStats* stats) -{ - psS32 i = 0; // Loop index variable. - - for (i = 0; i < myVector->n; i++) { - if (maskVector != NULL) - printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); - else - printf("Element %d is %f\n", i, myVector->data.F32[i]); } } @@ -428,8 +410,5 @@ // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; // Determine if we must only use data points within a min/max range. @@ -579,7 +558,5 @@ // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; // Determine if we must only use data points within a min/max range. @@ -817,4 +794,6 @@ p_psNormalizeVectorF32(myData): this is a private function which normalizes the elements of a vector to a range between 0.0 and 1.0. + +XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)? *****************************************************************************/ void p_psNormalizeVectorF32(psVector* myData) @@ -847,4 +826,6 @@ elements of a vector to a range between -1.0 and 1.0. XXX: 0-1 or -1:1? + +XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)? *****************************************************************************/ void p_psNormalizeVectorF64(psVector* myData) @@ -876,4 +857,5 @@ } +// XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)? void p_psNormalizeVector(psVector* myData) { @@ -885,87 +867,4 @@ psError(__func__, "Unalowable data type.\n"); } -} - -/***************************************************************************** -p_psGaussian(myData, myParams): an internal function, used by robust stats, -intended to compute the Gaussian with a specified sigma/mean, at the -specified data point. - *****************************************************************************/ -float p_psGaussian(const psVector* restrict myData, - const psVector* restrict myParams) -{ - float x = myData->data.F32[0]; - float mean = myParams->data.F32[0]; - float stdev = myParams->data.F32[1]; - float tmp = exp(-((x - mean) * (x - mean)) / (2.0 * stdev * stdev)); - - tmp /= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); - - // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); - return (tmp); -} - -/***************************************************************************** -p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which -calculates the specified partial derivative of the above Gaussian function. - *****************************************************************************/ -float p_psGaussianDeriv(const psVector* restrict myData, - const psVector* restrict myParams, psS32 whichParam) -{ - float x = myData->data.F32[0]; - float mean = myParams->data.F32[0]; - float stdev = myParams->data.F32[1]; - float tmp = 0.0; - - if (whichParam == 0) { - // Return the derivative w.r.t. the mean. - tmp = (x - mean) * p_psGaussian(myData, myParams); - tmp /= (stdev * stdev); - } else if (whichParam == 1) { - // Return the derivative w.r.t. the stdev. - tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams); - tmp /= (stdev * stdev * stdev); - } - printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp); - - return (tmp); -} - -/***************************************************************************** -p_psQuadratic(myData, myParams): an internal function, used by robust stats, -intended to compute a quadratic, with the specified parameters, at the -specified data point. - *****************************************************************************/ -float p_psQuadratic(const psVector* restrict myParams, const psVector* restrict myCoords) -{ - float x = myCoords->data.F32[0]; - float A = myParams->data.F32[0]; - float B = myParams->data.F32[1]; - float C = myParams->data.F32[2]; - float tmp = 0.0; - - tmp = (A * x * x) + (B * x) + C; - return (tmp); -} - -/***************************************************************************** -p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which -calculates the specified partial derivative of the above quadratic function. - *****************************************************************************/ -float p_psQuadraticDeriv(const psVector* restrict myParams, - const psVector* restrict myCoords, psS32 whichParamDeriv) -{ - float x = myCoords->data.F32[0]; - float tmp = 0.0; - - if (whichParamDeriv == 0) { - tmp = x * x; - } else if (whichParamDeriv == 1) { - tmp = x; - } else if (whichParamDeriv == 2) { - tmp = 1.0; - } - - return (tmp); } @@ -1022,6 +921,4 @@ and i+1 is used for x[i]). It then determines for what value x does that quadratic f(x) = yVal (the input parameter). - -// XXX: Use static variables. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1030,26 +927,13 @@ float yVal) { - // static psVector* x = NULL; - // static psVector* y = NULL; - // static psVector* yErr = NULL; - psVector* x = psVectorAlloc(3, PS_TYPE_F64); - psVector* y = psVectorAlloc(3, PS_TYPE_F64); - psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); - - /* - if (x == NULL) { - x = psVectorAlloc(3, PS_TYPE_F64); - p_psMemSetPersistent(x, true); - } - if (y == NULL) { - y = psVectorAlloc(3, PS_TYPE_F64); - p_psMemSetPersistent(y, true); - } - if (yErr == NULL) { - yErr = psVectorAlloc(3, PS_TYPE_F64); - p_psMemSetPersistent(yErr, true); - } - */ - psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); + PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); + PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); + PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); + static psPolynomial1D* myPoly = NULL; + + if (myPoly == NULL) { + myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); + p_psMemSetPersistent(myPoly, true); + } float tmpFloat; @@ -1090,8 +974,4 @@ } - psFree(x); - psFree(y); - psFree(yErr); - psFree(myPoly); return(tmpFloat); } @@ -1475,5 +1355,5 @@ /****************************************************************************** -psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure +psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure with the specifed bounds. Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 2227) +++ /trunk/psLib/src/math/psConstants.h (revision 2228) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-28 00:22:53 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -208,4 +208,10 @@ p_psMemSetPersistent(NAME->data.V, true); \ +#define PS_VECTOR_DECLARE_ALLOC_STATIC(NAME, SIZE, TYPE) \ +static psVector *NAME = NULL; \ +if (NAME == NULL) { \ + NAME = psVectorAlloc(SIZE, TYPE); \ + p_psMemSetPersistent(NAME, true); \ +} \ Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 2227) +++ /trunk/psLib/src/math/psMinimize.c (revision 2228) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.72 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:43 $ + * @version $Revision: 1.73 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -408,4 +408,54 @@ } + +/****************************************************************************** +XXX: We assume unnormalized gaussians. +XXX: Currently, yErr is ignored. + *****************************************************************************/ +psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, + const psVector *params, + const psArray *coords) +{ + psTrace(".psLib.dataManip.psMinimize", 4, + "---- psMinimizeLMChi2Gauss1D() begin ----\n"); + + PS_PTR_CHECK_NULL(coords, NULL); + PS_PTR_CHECK_NULL(params, NULL); + float x; + int i; + float mean = params->data.F32[0]; + float stdev = params->data.F32[1]; + psVector *out = psVectorAlloc(coords->n, PS_TYPE_F32); + + psTrace(".psLib.dataManip.psMinimize", 6, + "(mean, stdev) is (%f, %f)\n", mean, stdev); + + if (deriv == NULL) { + deriv = psImageAlloc(params->n, coords->n, PS_TYPE_F32); + } + + for (i=0;in;i++) { + 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]); + } + + for (i=0;in;i++) { + x = ((psVector *) (coords->data[i]))->data.F32[0]; + 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); + deriv->data.F32[i][1] = tmp / (stdev * stdev * stdev); + } + + psTrace(".psLib.dataManip.psMinimize", 4, + "---- psMinimizeLMChi2Gauss1D() end ----\n"); + return(out); +} + + /****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates @@ -429,5 +479,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. *****************************************************************************/ @@ -457,5 +507,5 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2() begin ----\n"); psS32 numData = y->n; @@ -490,8 +540,9 @@ float newChi2 = 0.0; float lamda = 0.00005; - - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + lamda = 0.05; + + psTrace(".psLib.dataManip.psMinimize", 6, "min->maxIter is %d\n", min->maxIter); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "min->tol is %f\n", min->tol); @@ -504,7 +555,7 @@ while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "------------------------------------------------------\n"); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); @@ -515,12 +566,16 @@ currValueVec = func(deriv, params, x); for (n=0;ndata.F32[n] * currValueVec->data.F32[n]); + currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) * + (currValueVec->data.F32[n] - y->data.F32[n]); + psTrace(".psLib.dataManip.psMinimize", 6, + "data[%d], chi2 calculation+= (%f * %f)\n", n, + currValueVec->data.F32[n], y->data.F32[n]); } for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "Current chi-squared is (%f)\n", currChi2); @@ -548,5 +603,5 @@ // XXX: multiple by -1 here? (beta->data.F64[p])*= -1.0; - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); } @@ -580,9 +635,9 @@ } for (j=0;jdata.F64[j][k]); - } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n"); + psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]); + } + psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n"); } @@ -597,5 +652,5 @@ // for (i=0;idata.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]); if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { @@ -607,8 +662,8 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + psTrace(".psLib.dataManip.psMinimize", 6, "Calling func() with new parameters:\n"); for (i=0;idata.F32[%d] is %f.\n", i, newParams->data.F32[i]); } @@ -621,9 +676,11 @@ newValueVec = func(deriv, newParams, x); for (n=0;ndata.F32[n] * newValueVec->data.F32[n]); + newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) * + (newValueVec->data.F32[n] - y->data.F32[n]); + } psFree(newValueVec); - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2); @@ -640,8 +697,10 @@ } lamda*= 0.1; + psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n"); } else { lamda*= 10.0; - } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n"); + } + psTrace(".psLib.dataManip.psMinimize", 4, "lamda is %f\n", lamda); min->iter++; @@ -662,5 +721,5 @@ } - psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + psTrace(".psLib.dataManip.psMinimize", 4, "---- psMinimizeLMChi2() end (false) ----\n"); return(false); @@ -942,4 +1001,14 @@ } +#define PS_VECTOR_ADD_MULTIPLE(BASE, BASEMASK, LINE, OUT, MUL) \ +for (int i=0;in;i++) { \ + if (BASEMASK->data.U8[i] == 0) { \ + OUT->data.F32[i] = BASE->data.F32[i] + (MUL * LINE->data.F32[i]); \ + } else { \ + OUT->data.F32[i] = BASE->data.F32[i]; \ + } \ +} \ + + /****************************************************************************** p_psDetermineBracket(): This routine takes as input an arbitrary function, @@ -947,18 +1016,26 @@ function produces as output a bracket [a, b, c] such that - f(param + b * line) is less than f(param + a * line) and - f(param + c * line). - -Algorithm: XXX completely ad hoc: start with the user-supplied starting +f(param + b * line) is less than f(param + a * line) and +f(param + c * line). + +Algorithm: +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. + +XXX: +output bracket vector should be an input as well. +*****************************************************************************/ psVector *p_psDetermineBracket(psVector *params, psVector *line, @@ -966,4 +1043,155 @@ const psArray *coords, psMinimizePowellFunc func) +{ + float a = 0.0; + float b = 0.0; + float c = 0.0; + float fa = 0.0; + float fb = 0.0; + float fc = 0.0; + psS32 iter = 100; + float aDir = 0.0; + float cDir = 0.0; + float new_aDir = 0.0; + float new_cDir = 0.0; + psVector *bracket = psVectorAlloc(3, PS_TYPE_F32); + float stepSize = PS_DETERMINE_BRACKET_STEP_SIZE; + psVector *tmp = NULL; + psS32 i = 0; + psS32 null = 0; + + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() begin ----\n"); + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + if (line->data.F32[i] >= FLT_EPSILON) { + null = 1; + } + } + } + + if (null == 0) { + psTrace(".psLib.dataManip.p_psDetermineBracket", 2, + "p_psDetermineBracket() called with zero line vector.\n"); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end (NULL) ----\n"); + psFree(bracket); + return(NULL); + } + + tmp = psVectorAlloc(params->n, PS_TYPE_F32); + + b = 0; + a = -stepSize; + c = stepSize; + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, a); + fa = func(tmp, coords); + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, b); + fb = func(tmp, coords); + + PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, tmp, c); + fc = func(tmp, coords); + + if (fa < fb) { + aDir = -1; + } else { + aDir = 1; + } + + if (fc < fb) { + cDir = -1; + } else { + cDir = 1; + } + + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "(a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", a, b, c, fa, fb, fc); + + while (iter > 0) { + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "psDetermineBracket(): iteration %d\n", iter); + if ((fb < fa) && (fb < fc)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + stepSize*= (1.0 + stepSize); + a =- stepSize; + c =+ stepSize; + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + tmp->data.F32[i] = params->data.F32[i] + (a * line->data.F32[i]); + } else { + tmp->data.F32[i] = params->data.F32[i]; + } + } + fa = func(tmp, coords); + + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + tmp->data.F32[i] = params->data.F32[i] + (c * line->data.F32[i]); + } else { + tmp->data.F32[i] = params->data.F32[i]; + } + } + fc = func(tmp, coords); + + psTrace(".psLib.dataManip.p_psDetermineBracket", 6, + "Iter(%d): (a, b, c) is (%f %f %f) (fa, fb, fc) is (%f %f %f)\n", iter, a, b, c, fa, fb, fc); + + if (fa < fb) { + new_aDir = -1; + } else { + new_aDir = 1; + } + + if (fc < fb) { + new_cDir = -1; + } else { + new_cDir = 1; + } + if ((new_aDir == 1) && (aDir == -1)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + + if ((new_cDir == 1) && (cDir == -1)) { + bracket->data.F32[0] = a; + bracket->data.F32[1] = b; + bracket->data.F32[2] = c; + psFree(tmp); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end ----\n"); + return(bracket); + } + aDir = new_aDir; + cDir = new_cDir; + iter--; + } + psFree(tmp); + psFree(bracket); + psTrace(".psLib.dataManip.p_psDetermineBracket", 4, + "---- p_psDetermineBracket() end (NULL) ----\n"); + return(NULL); +} + +psVector *p_psDetermineBracket2(psVector *params, + psVector *line, + const psVector *paramMask, + const psArray *coords, + psMinimizePowellFunc func) { float a = 0.0; @@ -1365,4 +1593,5 @@ float fqp = 0.0; float diff = 0.0; + int iterationNumber = 0; psTrace(".psLib.dataManip.psMinimizePowell", 4, @@ -1407,9 +1636,8 @@ } - min->iter = 0; - while (min->iter < min->maxIter) { - min->iter++; + while (iterationNumber < min->maxIter) { + iterationNumber++; psTrace(".psLib.dataManip.psMinimizePowell", 6, - "psMinimizePowell() iteration %d\n", min->iter); + "psMinimizePowell() iteration %d\n", iterationNumber); // 3: For each dimension in params, move Q only in the vector v[i] to @@ -1418,4 +1646,7 @@ baseFuncVal = func(Q, coords); currFuncVal = baseFuncVal; + psTrace(".psLib.dataManip.psMinimizePowell", 6, + "Current function value is %f\n", currFuncVal); + diff = 0.0; biggestDiff = 0; @@ -1423,5 +1654,5 @@ for (i=0;idata.U8[i] == 0) { - dummyMin.maxIter = 100; + dummyMin.maxIter = 20; dummyMin.tol = 0.01; mul = p_psLineMin(&dummyMin, Q, v[i], myParamMask, coords, func); @@ -1433,5 +1664,6 @@ } } - + psTrace(".psLib.dataManip.psMinimizePowell", 6, + "New function value is %f\n", currFuncVal); // 4: Set the vector u = Q - P for (i=0;iiter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (true) ----\n"); @@ -1513,9 +1745,10 @@ psFree(pQP); psFree(Q); + if (paramMask == NULL) { + psFree(myParamMask); + } + min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (true) ----\n"); - if (paramMask == NULL) { - psFree(myParamMask); - } return(true); } @@ -1532,4 +1765,5 @@ psFree(pQP); psFree(Q); + min->iter = iterationNumber; psTrace(".psLib.dataManip.psMinimizePowell", 4, "---- psMinimizePowell() end (false) ----\n"); @@ -1549,4 +1783,12 @@ const psArray *coords) { + psTrace(".psLib.dataManip.myPowellChi2Func", 4, + "---- myPowellChi2Func() begin ----\n"); + PS_VECTOR_CHECK_NULL(params, NAN); + PS_VECTOR_CHECK_EMPTY(params, NAN); + PS_VECTOR_CHECK_NULL(myValue, NAN); + PS_VECTOR_CHECK_EMPTY(myValue, NAN); + PS_PTR_CHECK_NULL(coords, NAN); + float chi2 = 0.0; float d; @@ -1555,9 +1797,18 @@ tmp = Chi2PowellFunc(params, coords); - for (i=0;in;i++) { - d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i]; - chi2+= d * d; + if (myError == NULL) { + for (i=0;in;i++) { + d = (tmp->data.F32[i] - myValue->data.F32[i]); + chi2+= d * d; + } + } else { + for (i=0;in;i++) { + d = (tmp->data.F32[i] - myValue->data.F32[i]) / myError->data.F32[i]; + chi2+= d * d; + } } psFree(tmp); + psTrace(".psLib.dataManip.myPowellChi2Func", 4, + "---- myPowellChi2Func() end (%f) ----\n", chi2); return(chi2); } @@ -1584,4 +1835,5 @@ myValue = (psVector *) value; myError = (psVector *) error; + Chi2PowellFunc = func; Index: /trunk/psLib/src/math/psMinimize.h =================================================================== --- /trunk/psLib/src/math/psMinimize.h (revision 2227) +++ /trunk/psLib/src/math/psMinimize.h (revision 2228) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.34 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:43 $ + * @version $Revision: 1.35 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -81,41 +81,4 @@ const psArray *x); -/* -psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, - const psVector *params, - const psArray *coords) -{ - psVector *x = (psVector *) coords->data[0]; - psVector *y = (psVector *) coords->data[1]; - psVector *yErr = NULL; - int i; - int j; - float normalization = params->data.F32[0]; - float center = params->data.F32[1]; - float stdev = params->data.F32[2]; - psVector *out = psVectorAlloc(x->n, PS_TYPE_F32); - - - if (deriv == NULL) { - deriv = psImageAlloc(params->n, x->n, PS_TYPE_F32); - } - - if (coords->n == 2) { - yErr = (psVector *) coords->data[2]; - } - - for (i=0;in;i++) { - out->data.F32[i] = psGaussian(x->data.F32[i], center, stdev, 1); - } - - for (i=0;in;i++) { - for (j=0;jn;j++) { - deriv->data.F32[i][j] = the derivative of the gaussian w.r.t. the j-th parameter at the coord x[i] - } - } - return(out); -} -*/ - psBool psMinimizeLMChi2(psMinimization *min, psImage *covar, @@ -137,4 +100,8 @@ psMinimizePowellFunc func); +psVector *psMinimizeLMChi2Gauss1D(psImage *deriv, + const psVector *params, + const psArray *coords); + typedef psVector *(*psMinimizeChi2PowellFunc)(const psVector *params, Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 2227) +++ /trunk/psLib/src/math/psStats.c (revision 2228) @@ -9,7 +9,7 @@ * @author GLG, MHPCC * - * @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-28 00:46:04 $ -n * + * @version $Revision: 1.80 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:46:57 $ + * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ @@ -123,22 +123,4 @@ default: return false; - } -} - -/***************************************************************************** -p_psVectorPrint(myVector, maskVector, maskVal, stats): a private internal -function that simply prints a vector to STDOUT. Used primarily for -debugging. - *****************************************************************************/ - -void p_psVectorPrint(psVector* myVector, psVector* maskVector, psU32 maskVal, psStats* stats) -{ - psS32 i = 0; // Loop index variable. - - for (i = 0; i < myVector->n; i++) { - if (maskVector != NULL) - printf("Element %d is %f (mask is %d)\n", i, myVector->data.F32[i], maskVector->data.U8[i]); - else - printf("Element %d is %f\n", i, myVector->data.F32[i]); } } @@ -428,8 +410,5 @@ // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; - sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; // Determine if we must only use data points within a min/max range. @@ -579,7 +558,5 @@ // Allocate temporary vectors for the data. unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - unsortedVector->n = unsortedVector->nalloc; sortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - sortedVector->n = sortedVector->nalloc; // Determine if we must only use data points within a min/max range. @@ -817,4 +794,6 @@ p_psNormalizeVectorF32(myData): this is a private function which normalizes the elements of a vector to a range between 0.0 and 1.0. + +XXX: how about an arbitrary p_psNormalizeVectorF32(myData, min, max)? *****************************************************************************/ void p_psNormalizeVectorF32(psVector* myData) @@ -847,4 +826,6 @@ elements of a vector to a range between -1.0 and 1.0. XXX: 0-1 or -1:1? + +XXX: how about an arbitrary p_psNormalizeVectorF64(myData, min, max)? *****************************************************************************/ void p_psNormalizeVectorF64(psVector* myData) @@ -876,4 +857,5 @@ } +// XXX: how about an arbitrary p_psNormalizeVector(myData, min, max)? void p_psNormalizeVector(psVector* myData) { @@ -885,87 +867,4 @@ psError(__func__, "Unalowable data type.\n"); } -} - -/***************************************************************************** -p_psGaussian(myData, myParams): an internal function, used by robust stats, -intended to compute the Gaussian with a specified sigma/mean, at the -specified data point. - *****************************************************************************/ -float p_psGaussian(const psVector* restrict myData, - const psVector* restrict myParams) -{ - float x = myData->data.F32[0]; - float mean = myParams->data.F32[0]; - float stdev = myParams->data.F32[1]; - float tmp = exp(-((x - mean) * (x - mean)) / (2.0 * stdev * stdev)); - - tmp /= ((float)sqrt(2.0 * M_PI * (stdev * stdev))); - - // printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp); - return (tmp); -} - -/***************************************************************************** -p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which -calculates the specified partial derivative of the above Gaussian function. - *****************************************************************************/ -float p_psGaussianDeriv(const psVector* restrict myData, - const psVector* restrict myParams, psS32 whichParam) -{ - float x = myData->data.F32[0]; - float mean = myParams->data.F32[0]; - float stdev = myParams->data.F32[1]; - float tmp = 0.0; - - if (whichParam == 0) { - // Return the derivative w.r.t. the mean. - tmp = (x - mean) * p_psGaussian(myData, myParams); - tmp /= (stdev * stdev); - } else if (whichParam == 1) { - // Return the derivative w.r.t. the stdev. - tmp = (x - mean) * (x - mean) * p_psGaussian(myData, myParams); - tmp /= (stdev * stdev * stdev); - } - printf("p_psGaussianDeriv((%.2f), %.2f, %.2f, (%d)) is %.2f\n", x, mean, stdev, whichParam, tmp); - - return (tmp); -} - -/***************************************************************************** -p_psQuadratic(myData, myParams): an internal function, used by robust stats, -intended to compute a quadratic, with the specified parameters, at the -specified data point. - *****************************************************************************/ -float p_psQuadratic(const psVector* restrict myParams, const psVector* restrict myCoords) -{ - float x = myCoords->data.F32[0]; - float A = myParams->data.F32[0]; - float B = myParams->data.F32[1]; - float C = myParams->data.F32[2]; - float tmp = 0.0; - - tmp = (A * x * x) + (B * x) + C; - return (tmp); -} - -/***************************************************************************** -p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which -calculates the specified partial derivative of the above quadratic function. - *****************************************************************************/ -float p_psQuadraticDeriv(const psVector* restrict myParams, - const psVector* restrict myCoords, psS32 whichParamDeriv) -{ - float x = myCoords->data.F32[0]; - float tmp = 0.0; - - if (whichParamDeriv == 0) { - tmp = x * x; - } else if (whichParamDeriv == 1) { - tmp = x; - } else if (whichParamDeriv == 2) { - tmp = 1.0; - } - - return (tmp); } @@ -1022,6 +921,4 @@ and i+1 is used for x[i]). It then determines for what value x does that quadratic f(x) = yVal (the input parameter). - -// XXX: Use static variables. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1030,26 +927,13 @@ float yVal) { - // static psVector* x = NULL; - // static psVector* y = NULL; - // static psVector* yErr = NULL; - psVector* x = psVectorAlloc(3, PS_TYPE_F64); - psVector* y = psVectorAlloc(3, PS_TYPE_F64); - psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); - - /* - if (x == NULL) { - x = psVectorAlloc(3, PS_TYPE_F64); - p_psMemSetPersistent(x, true); - } - if (y == NULL) { - y = psVectorAlloc(3, PS_TYPE_F64); - p_psMemSetPersistent(y, true); - } - if (yErr == NULL) { - yErr = psVectorAlloc(3, PS_TYPE_F64); - p_psMemSetPersistent(yErr, true); - } - */ - psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); + PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); + PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); + PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); + static psPolynomial1D* myPoly = NULL; + + if (myPoly == NULL) { + myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); + p_psMemSetPersistent(myPoly, true); + } float tmpFloat; @@ -1090,8 +974,4 @@ } - psFree(x); - psFree(y); - psFree(yErr); - psFree(myPoly); return(tmpFloat); } @@ -1475,5 +1355,5 @@ /****************************************************************************** -psHistogramAllocGenric(bounds): allocate a non-uniform histogram structure +psHistogramAllocGeneric(bounds): allocate a non-uniform histogram structure with the specifed bounds. Index: /trunk/psLib/src/sys/psTrace.h =================================================================== --- /trunk/psLib/src/sys/psTrace.h (revision 2227) +++ /trunk/psLib/src/sys/psTrace.h (revision 2228) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:44 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:47:00 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -18,5 +18,5 @@ #define PS_UNKNOWN_TRACE_LEVEL -9999 // we don't know this name's level -#define PS_DEFAULT_TRACE_LEVEL 0 +#define PS_DEFAULT_TRACE_LEVEL -1 /** \addtogroup LogTrace Index: /trunk/psLib/src/sysUtils/psTrace.h =================================================================== --- /trunk/psLib/src/sysUtils/psTrace.h (revision 2227) +++ /trunk/psLib/src/sysUtils/psTrace.h (revision 2228) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-27 23:31:44 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-28 22:47:00 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -18,5 +18,5 @@ #define PS_UNKNOWN_TRACE_LEVEL -9999 // we don't know this name's level -#define PS_DEFAULT_TRACE_LEVEL 0 +#define PS_DEFAULT_TRACE_LEVEL -1 /** \addtogroup LogTrace