Index: /trunk/psLib/src/math/psMinimizeLMM.c =================================================================== --- /trunk/psLib/src/math/psMinimizeLMM.c (revision 6321) +++ /trunk/psLib/src/math/psMinimizeLMM.c (revision 6322) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-27 20:08:58 $ + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-03 22:05:22 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -48,5 +48,4 @@ // XXX EAM : can we use static copies of LUv, LUm, A? -// XXX: Add trace messages, check return codes. psBool p_psMinLM_GuessABP( psImage *Alpha, @@ -75,6 +74,7 @@ A = psImageCopy(NULL, alpha, PS_TYPE_F64); for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) + if ((paramMask != NULL) && (paramMask->data.U8[j])) { continue; + } A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); } @@ -83,7 +83,20 @@ // these operations do not modify the input values (creates LUm, LUv) LUm = psMatrixLUD(NULL, &LUv, A); + if (LUm == NULL) { + psError(PS_ERR_UNKNOWN, false, "psMatrixLUD() returned NULL\n"); + } Beta = psMatrixLUSolve(Beta, LUm, beta, LUv); + if (Beta == NULL) { + psError(PS_ERR_UNKNOWN, false, "psMatrixLUSolve() returned NULL\n"); + } Alpha = psMatrixInvert(Alpha, A, &det); - + psFree(LUm); + psFree(LUv); + psFree(A); + + if (Alpha == NULL) { + psError(PS_ERR_UNKNOWN, false, "psMatrixInvert() returned NULL\n"); + return(false); + } # else // gauss-jordan version @@ -94,20 +107,21 @@ Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F64); for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) + if ((paramMask != NULL) && (paramMask->data.U8[j])) { continue; + } Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); } - // XXX: Check error codes! - psGaussJordan(Alpha, Beta); - psFree(A); - psFree(LUm); - psFree(LUv); + if (false == psGaussJordan(Alpha, Beta)) { + psError(PS_ERR_UNKNOWN, false, "psMatrixInvert() returned NULL\n"); + return(false); + } # endif // apply Beta to get new Params values for (int j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) + if ((paramMask != NULL) && (paramMask->data.U8[j])) { continue; + } // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; // compare Beta to beta limits @@ -126,9 +140,4 @@ } } - # if (USE_LU_DECOMP) - psFree(A); - psFree(LUm); - psFree(LUv); - # endif return(true); @@ -136,5 +145,4 @@ -// XXX: Add trace messages, check return codes. bool psMinimizeGaussNewtonDelta( psVector *delta, @@ -146,4 +154,5 @@ psMinimizeLMChi2Func func) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); // allocate internal arrays (current vs Guess) psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); @@ -152,4 +161,5 @@ psVector *Params = psVectorAlloc(params->n, PS_TYPE_F64); psVector *dy = NULL; + psBool rc = true; // the user provides the error or NULL. we need to convert @@ -162,6 +172,17 @@ } - p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func); - p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0); + psF64 rcF64 = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func); + if (isnan(rcF64)) { + psError(PS_ERR_UNKNOWN, false, "p_psMinLM_SetABX() retruned a NAN.\n"); + rc = false; + } + psTrace(__func__, 5, "p_psMinLM_SetABX() was succesful\n", __func__); + + psBool rcBool = p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0); + if (rcBool == false) { + psError(PS_ERR_UNKNOWN, false, "p_psMinLM_GuessABP() retruned FALSE.\n"); + rc = false; + } + psTrace(__func__, 5, "p_psMinLM_GuessABP() was succesful\n", __func__); psFree(alpha); @@ -172,5 +193,6 @@ psFree(dy); } - return (true); + psTrace(__func__, 3, "---- %s() end ----\n", __func__); + return(rc); } @@ -235,10 +257,12 @@ for (psS32 j = 0; j < params->n; j++) { - if ((paramMask != NULL) && (paramMask->data.U8[j])) + if ((paramMask != NULL) && (paramMask->data.U8[j])) { continue; + } weight = deriv->data.F32[j] * dy->data.F32[i]; for (psS32 k = 0; k <= j; k++) { - if ((paramMask != NULL) && (paramMask->data.U8[k])) + if ((paramMask != NULL) && (paramMask->data.U8[k])) { continue; + } alpha->data.F64[j][k] += weight * deriv->data.F32[k]; } @@ -394,5 +418,8 @@ } if (psTraceGetLevel (__func__) >= 6) { - //XXX: p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess"); + psTrace(__func__, 6, "The current Param vector: \n"); + for (psS32 i = 0 ; i < Params->n ; i++) { + psTrace(__func__, 6, "Params[%d] is %f\n", Params->data.F32[i]); + } } @@ -416,5 +443,10 @@ } if (psTraceGetLevel(__func__) >= 6) { - //XXX: p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess"); + if (psTraceGetLevel (__func__) >= 6) { + psTrace(__func__, 6, "The current Param vector: \n"); + for (psS32 i = 0 ; i < Params->n ; i++) { + psTrace(__func__, 6, "Params[%d] is %f\n", Params->data.F32[i]); + } + } } @@ -628,4 +660,4 @@ bool psMemCheckConstrain(psPtr tmp) { - return( psMemGetDeallocator(tmp) == (psFreeFunc)constrainFree ); -} + return(psMemGetDeallocator(tmp) == (psFreeFunc) constrainFree); +} Index: /trunk/psLib/src/math/psMinimizePowell.c =================================================================== --- /trunk/psLib/src/math/psMinimizePowell.c (revision 6321) +++ /trunk/psLib/src/math/psMinimizePowell.c (revision 6322) @@ -9,11 +9,8 @@ * @author GLG, MHPCC * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-23 22:25:31 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-03 22:05:22 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii - * - * XXX: must follow coding name standards on local functions. - * XXX: put local functions in front. * */ @@ -69,9 +66,8 @@ /*****************************************************************************/ /* GLOBAL VARIABLES */ -/* XXX: Do these conform to code standard? */ /*****************************************************************************/ static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL; -static psVector *myValue; -static psVector *myError; +static psVector *PowellValue; +static psVector *PowellError; /*****************************************************************************/ /* FILE STATIC VARIABLES */ @@ -97,6 +93,5 @@ Repeat this process until a local minimum is found. -XXX: -new algorithm: +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 @@ -104,11 +99,8 @@ direction. -XXX: -This is F32 only. Must add -F64 support (actually, make the defaults F64, +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( @@ -355,10 +347,11 @@ *****************************************************************************/ #define PS_LINEMIN_MAX_ITERATIONS 30 -psF32 p_psLineMin(psMinimization *min, - psVector *params, - psVector *line, - const psVector *paramMask, - const psArray *coords, - psMinimizePowellFunc func) +psF32 p_psLineMin( + psMinimization *min, + psVector *params, + psVector *line, + const psVector *paramMask, + const psArray *coords, + psMinimizePowellFunc func) { PS_ASSERT_PTR_NON_NULL(min, NAN); @@ -502,9 +495,10 @@ #define PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE 0.01 -bool psMinimizePowell(psMinimization *min, - psVector *params, - const psVector *paramMask, - const psArray *coords, - psMinimizePowellFunc func) +bool psMinimizePowell( + psMinimization *min, + psVector *params, + const psVector *paramMask, + const psArray *coords, + psMinimizePowellFunc func) { PS_ASSERT_PTR_NON_NULL(min, NULL); @@ -696,14 +690,16 @@ This functions uses global variables to receive the function pointer, the data values, and the data errors. + XXX: This is F32 only. Must implement F64. *****************************************************************************/ -static psF32 myPowellChi2Func(const psVector *params, - const psArray *coords) +static psF32 myPowellChi2Func( + const psVector *params, + const psArray *coords) { psTrace(__func__, 4, "---- myPowellChi2Func() begin ----\n"); PS_ASSERT_VECTOR_NON_NULL(params, NAN); PS_ASSERT_VECTOR_NON_EMPTY(params, NAN); - PS_ASSERT_VECTOR_NON_NULL(myValue, NAN); - PS_ASSERT_VECTOR_NON_EMPTY(myValue, NAN); + PS_ASSERT_VECTOR_NON_NULL(PowellValue, NAN); + PS_ASSERT_VECTOR_NON_EMPTY(PowellValue, NAN); PS_ASSERT_PTR_NON_NULL(coords, NAN); @@ -714,12 +710,12 @@ tmp = Chi2PowellFunc(params, coords); - if (myError == NULL) { + if (PowellError == NULL) { for (i=0;in;i++) { - d = (tmp->data.F32[i] - myValue->data.F32[i]); + d = (tmp->data.F32[i] - PowellValue->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]; + d = (tmp->data.F32[i] - PowellValue->data.F32[i]) / PowellError->data.F32[i]; chi2+= d * d; } @@ -741,14 +737,15 @@ psMinimizePowell(). *****************************************************************************/ -bool psMinimizeChi2Powell(psMinimization *min, - psVector *params, - const psVector *paramMask, - const psArray *coords, - const psVector *value, - const psVector *error, - psMinimizeChi2PowellFunc model) +bool psMinimizeChi2Powell( + psMinimization *min, + psVector *params, + const psVector *paramMask, + const psArray *coords, + const psVector *value, + const psVector *error, + psMinimizeChi2PowellFunc model) { - myValue = (psVector *) value; - myError = (psVector *) error; + PowellValue = (psVector *) value; + PowellError = (psVector *) error; Chi2PowellFunc = model; Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 6321) +++ /trunk/psLib/src/math/psStats.c (revision 6322) @@ -16,6 +16,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.166 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-03 01:30:14 $ + * @version $Revision: 1.167 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-02-03 22:05:22 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -303,6 +303,4 @@ max of the input vector. If there was a problem with the max calculation, this routine sets stats->max to NAN. - -XXX: Do we need to factor errors into it? *****************************************************************************/ psS32 p_psVectorMax(const psVector* myVector, @@ -549,6 +547,4 @@ median of the input vector. Returns true on success (including if there were no valid input vector elements). - -XXX: Use static vectors for sort arrays. *****************************************************************************/ bool p_psVectorSampleMedian(const psVector* myVector, @@ -783,12 +779,12 @@ { psTrace(__func__, 4, "---- %s() begin ----\n", __func__); - psVector* unsortedVector = NULL; // Temporary vector - psVector* sortedVector = NULL; // Temporary vector - psS32 i = 0; // Loop index variable - psS32 count = 0; // # of points in this mean. - psS32 nValues = 0; // # data points + psVector* unsortedVector = NULL; + psVector* sortedVector = NULL; + psS32 i = 0; + psS32 count = 0; + psS32 nValues = 0; // Determine how many data points fit inside this min/max range - // and are not maxed, IF the maskVector is not NULL. + // and are not masked, if the maskVector is not NULL. nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats); @@ -796,5 +792,4 @@ unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32); - // Determine if we must only use data points within a min/max range. if (stats->options & PS_STAT_USE_RANGE) { // Store all non-masked data points within the min/max range @@ -1001,9 +996,9 @@ psTrace(__func__, 4, "---- %s() begin ----\n", __func__); psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__)); - psF32 clippedMean = 0.0; // self-explanatory - psF32 clippedStdev = 0.0; // self-explanatory - psVector* tmpMask = NULL; // Temporary vector for masks during iterations. - static psStats *statsTmp = NULL; // Temporary psStats struct. - psS32 rc = 0; // Return code. + psF32 clippedMean = 0.0; + psF32 clippedStdev = 0.0; + psVector* tmpMask = NULL; + static psStats *statsTmp = NULL; + psS32 rc = 0; // Ensure that stats->clipIter is within the proper range. @@ -1189,12 +1184,4 @@ parameter). -XXX: After you fit the polynomial, solve for X analytically. - -XXX: the vectors do not have to be the same length. Must insert the proper -tests to ensure that binNum is within acceptable ranges for both vectors. - -XXX: This currently assumes that the three points are monotonically increasing -or decreasing: so, it works for the cumulative histogram vectors, but not for -arbitrary vectors. We should probably test that condition. *****************************************************************************/ psF32 fitQuadraticSearchForYThenReturnX( @@ -1221,9 +1208,7 @@ psVector *x = psVectorAlloc(3, PS_TYPE_F64); psVector *y = psVectorAlloc(3, PS_TYPE_F64); - psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); - psF32 tmpFloat = 0.0f; - if ((binNum > 0) && (binNum < (yVec->n - 2))) { + if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) { // The general case. We have all three points. x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); @@ -1234,9 +1219,7 @@ y->data.F64[2] = yVec->data.F32[binNum + 1]; psTrace(__func__, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], - xVec->data.F32[binNum], - xVec->data.F32[binNum+1], - xVec->data.F32[binNum+2]); - psTrace(__func__, 6, "x vec is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); - psTrace(__func__, 6, "y vec is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); + xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]); + psTrace(__func__, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); + psTrace(__func__, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); // @@ -1246,50 +1229,28 @@ // so that the following checks are not necessary. // - if (y->data.F64[0] < y->data.F64[1]) { - if (!(y->data.F64[1] <= y->data.F64[2])) { - psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); - psFree(x); - psFree(y); - psFree(yErr); - psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - // Ensure that yVal is within the range of the bins we are using. - if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, - (psF64)yVal, y->data.F64[2], y->data.F64[0]); - } - } else if (y->data.F64[0] > y->data.F64[1]) { - if (!(y->data.F64[1] >= y->data.F64[2])) { - psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); - psFree(x); - psFree(y); - psFree(yErr); - psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - // Ensure that yVal is within the range of the bins we are using. - if (!((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, - (psF64)yVal, y->data.F64[2], y->data.F64[0]); - } - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; + if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || + ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, + (psF64)yVal, y->data.F64[0], y->data.F64[2]); + } + + if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || + ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { + psError(PS_ERR_UNKNOWN, true, "This routine must be called with monotonically increasing or decreasing data points.\n"); + psFree(x); + psFree(y); + psTrace(__func__, 5, "---- %s() end ----\n", __func__); + return(NAN); + } // Determine the coefficients of the polynomial. psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x); + myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x); if (myPoly == NULL) { - psError(PS_ERR_UNEXPECTED_NULL, - false, + psError(PS_ERR_UNEXPECTED_NULL, false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT); psFree(x); psFree(y); - psFree(yErr); psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__); return(NAN); @@ -1309,17 +1270,14 @@ if (isnan(tmpFloat)) { psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN); + false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN); psFree(x); psFree(y); - psFree(yErr); psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__); return(NAN); } - } else { - // The special case where we have two points only at the beginning of - // the vectors x and y. + // These are special cases where the bin is at the beginning or end of the vector. if (binNum == 0) { + // We have two points only at the beginning of the vectors x and y. tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); @@ -1339,5 +1297,4 @@ psFree(x); psFree(y); - psFree(yErr); psTrace(__func__, 5, "---- %s(%f) end ----\n", __func__, tmpFloat); @@ -1542,5 +1499,4 @@ // ADD: Step 3. // Interpolate to the exact 50% position: this is the robust histogram median. - // XXX: Check return codes. // stats->robustMedian = fitQuadraticSearchForYThenReturnX( @@ -1549,4 +1505,13 @@ binMedian, totalDataPoints/2.0); + if (isnan(stats->robustMedian)) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpMaskVec); + psTrace(__func__, 4, "---- %s(1) end ----\n", __func__); + return(1); + } psTrace(__func__, 6, "Current robust median is %f\n", stats->robustMedian); @@ -2324,8 +2289,9 @@ binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); if (errorsF32 != NULL) { - // XXX: Check return codes. - UpdateHistogramBins(binNum, out, - inF32->data.F32[i], - errorsF32->data.F32[i]); + psS32 rc = UpdateHistogramBins(binNum, out, inF32->data.F32[i], + errorsF32->data.F32[i]); + if (rc < 0) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n"); + } } else { // XXX: This if-statement really shouldn't be necessary. @@ -2348,8 +2314,10 @@ } else { if (errorsF32 != NULL) { - // XXX: Check return codes. - UpdateHistogramBins(binNum, out, - inF32->data.F32[i], - errors->data.F32[i]); + psS32 rc = UpdateHistogramBins(binNum, out, + inF32->data.F32[i], + errors->data.F32[i]); + if (rc < 0) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n"); + } } else { (out->nums->data.F32[binNum])+= 1.0;