Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 24960) +++ trunk/psLib/src/math/psStats.c (revision 24987) @@ -129,15 +129,15 @@ # define COUNT_WARNING(LIMIT, INTERVAL, ...) { \ - static int nCalls = 1; \ - if (nCalls < LIMIT) { \ - psWarning(__VA_ARGS__); \ - } \ - if (!(nCalls % INTERVAL)) { \ - psWarning(__VA_ARGS__); \ + static int nCalls = 1; \ + if (nCalls < LIMIT) { \ + psWarning(__VA_ARGS__); \ + } \ + if (!(nCalls % INTERVAL)) { \ + psWarning(__VA_ARGS__); \ psWarning("(warning raised %d times)", nCalls); \ - } \ - nCalls ++; \ + } \ + nCalls ++; \ } - + /*****************************************************************************/ @@ -217,6 +217,6 @@ for (long i = 0; i < numData; i++) { // Check if the data is with the specified range - if (!isfinite(data[i])) - continue; + if (!isfinite(data[i])) + continue; if (useRange && (data[i] < stats->min)) continue; @@ -242,5 +242,5 @@ if (!isnan(mean)) { - stats->results |= PS_STAT_SAMPLE_MEAN; + stats->results |= PS_STAT_SAMPLE_MEAN; } return true; @@ -278,6 +278,6 @@ for (long i = 0; i < num; i++) { // Check if the data is with the specified range - if (!isfinite(vector[i])) - continue; + if (!isfinite(vector[i])) + continue; if (useRange && (vector[i] < stats->min)) continue; @@ -329,6 +329,6 @@ // into the temporary vectors. for (long i = 0; i < inVector->n; i++) { - if (!isfinite(input[i])) - continue; + if (!isfinite(input[i])) + continue; if (useRange && (input[i] < stats->min)) continue; @@ -353,5 +353,5 @@ // Sort the temporary vector. if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK) - // an error in psVectorSort is a serious error + // an error in psVectorSort is a serious error psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data.")); stats->sampleUQ = NAN; @@ -407,5 +407,5 @@ // If the mean is NAN, then generate a warning and set the stdev to NAN. if (isnan(stats->sampleMean)) { - COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); + COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); stats->sampleStdev = NAN; return true; @@ -425,6 +425,6 @@ for (long i = 0; i < myVector->n; i++) { // Check if the data is with the specified range - if (!isfinite(data[i])) - continue; + if (!isfinite(data[i])) + continue; if (useRange && (data[i] < stats->min)) { continue; @@ -502,6 +502,6 @@ for (long i = 0; i < myVector->n; i++) { // Check if the data is with the specified range - if (!isfinite(data[i])) - continue; + if (!isfinite(data[i])) + continue; if (useRange && (data[i] < stats->min)) { continue; @@ -771,5 +771,5 @@ stats->robustLQ = min; stats->robustN50 = numValid; - // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits: + // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits: stats->results |= PS_STAT_ROBUST_MEDIAN; stats->results |= PS_STAT_ROBUST_STDEV; @@ -801,12 +801,12 @@ // XXXXX we need to consider this step if errors -> variance if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - // if psVectorHistogram returns false, we have a programming error + // if psVectorHistogram returns false, we have a programming error psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n"); - psFree(histogram); - psFree(cumulative); - psFree(statsMinMax); - psFree(mask); - - return false; + psFree(histogram); + psFree(cumulative); + psFree(statsMinMax); + psFree(mask); + + return false; } if (psTraceGetLevel("psLib.math") >= 8) { @@ -838,9 +838,9 @@ // ADD step 3: Interpolate to the exact 50% position in bin units - stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); - // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); - // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian); - - // convert bin to bin value: this is the robust histogram median. + stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian); + + // convert bin to bin value: this is the robust histogram median. if (isnan(stats->robustMedian)) { COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n"); @@ -1054,6 +1054,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1107,5 +1107,5 @@ psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - // if psVectorHistogram returns false, we have a programming error + // if psVectorHistogram returns false, we have a programming error psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n"); psFree(histogram); @@ -1172,6 +1172,6 @@ PS_VECTOR_PRINT_F32(y); } - - // psMinimizeLMChi2 can return false for bad data as well as for serious failures + + // psMinimizeLMChi2 can return false for bad data as well as for serious failures if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); @@ -1235,6 +1235,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1355,5 +1355,5 @@ // psVectorInit (fitMask, 0); - // XXX not sure if these should result in errors or not... + // XXX not sure if these should result in errors or not... if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) { psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); @@ -1432,6 +1432,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1729,6 +1729,6 @@ if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } + return false; + } } @@ -1766,5 +1766,5 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - goto escape; + goto escape; } @@ -1774,6 +1774,6 @@ stats->fittedMean = min; stats->fittedStdev = 0.0; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN_V4; + stats->results |= PS_STAT_FITTED_STDEV_V4; return true; } @@ -1793,5 +1793,5 @@ psFree(histogram); psFree(statsMinMax); - goto escape; + goto escape; } if (psTraceGetLevel("psLib.math") >= 8) { @@ -1825,5 +1825,5 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - goto escape; + goto escape; } @@ -1886,6 +1886,6 @@ // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 - // XXX this fit may fail with an error for an ill-conditioned matrix (bad data) - // we probably should be able to handle the data errors gracefully + // XXX this fit may fail with an error for an ill-conditioned matrix (bad data) + // we probably should be able to handle the data errors gracefully psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); @@ -1898,5 +1898,5 @@ psFree(histogram); psFree(statsMinMax); - goto escape; + goto escape; } @@ -1917,5 +1917,5 @@ COUNT_WARNING(10, 100, "fit did not converge\n"); - goto escape; + goto escape; } @@ -1987,5 +1987,5 @@ psFree(histogram); psFree(statsMinMax); - goto escape; + goto escape; } @@ -2898,7 +2898,7 @@ *****************************************************************************/ static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal + psVector *yVec, + psS32 binNum, + psF32 yVal ) { @@ -2982,27 +2982,35 @@ } - // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1 - assert (binValue >= binNum - 1); - assert (binValue <= binNum + 1); - - int fitBin = binValue; - float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; - float dY = binValue - fitBin; - tmpFloat = xVec->data.F32[fitBin] + dY * dX; + // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1 + assert (binValue >= binNum - 1); + assert (binValue <= binNum + 1); + + int fitBin = binValue; + float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; + float dY = binValue - fitBin; + tmpFloat = xVec->data.F32[fitBin] + dY * dX; } else { // 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. - // X = (dX/dY)(Y - Yo) + Xo - float dX = xVec->data.F32[1] - xVec->data.F32[0]; - float dY = yVec->data.F32[1] - yVec->data.F32[0]; - tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; + // X = (dX/dY)(Y - Yo) + Xo + float dX = xVec->data.F32[1] - xVec->data.F32[0]; + float dY = yVec->data.F32[1] - yVec->data.F32[0]; + if (dY == 0.0) { + tmpFloat = xVec->data.F32[0]; + } else { + tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; + } } else if (binNum == (xVec->n - 1)) { // We have two points only at the end of the vectors x and y. - // X = (dX/dY)(Y - Yo) + Xo - float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1]; - float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1]; - tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1]; - } + // X = (dX/dY)(Y - Yo) + Xo + float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1]; + float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1]; + if (dY == 0.0) { + tmpFloat = xVec->data.F32[binNum-1]; + } else { + tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1]; + } + } }