Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 6315) +++ 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;