Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 12494) +++ /trunk/psLib/src/math/psStats.c (revision 12495) @@ -13,6 +13,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.203 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-03-16 00:33:55 $ + * @version $Revision: 1.204 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-03-19 22:48:11 $ * * Copyright 2006 IfA, University of Hawaii @@ -67,24 +67,24 @@ (0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) -// set the bin closest to the corresponding value. if USE_END is +/- 1, +// set the bin closest to the corresponding value. if USE_END is +/- 1, // out-of-range saturates on lower/upper bin REGARDLESS of actual value #define PS_BIN_FOR_VALUE(RESULT, VECTOR, VALUE, USE_END) { \ - psVectorBinaryDisectResult result; \ - psScalar tmpScalar; \ - tmpScalar.type.type = PS_TYPE_F32; \ - tmpScalar.data.F32 = (VALUE); \ - RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \ - switch (result) { \ - case PS_BINARY_DISECT_PASS: \ - break; \ - case PS_BINARY_DISECT_OUTSIDE_RANGE: \ - psTrace ("psLib.math", 6, "selected bin outside range"); \ - if (USE_END == -1) { RESULT = 0; } \ - if (USE_END == +1) { RESULT = VECTOR->n - 1; } \ - break; \ - case PS_BINARY_DISECT_INVALID_INPUT: \ - case PS_BINARY_DISECT_INVALID_TYPE: \ - psAbort ("programming error"); \ - break; \ + psVectorBinaryDisectResult result; \ + psScalar tmpScalar; \ + tmpScalar.type.type = PS_TYPE_F32; \ + tmpScalar.data.F32 = (VALUE); \ + RESULT = psVectorBinaryDisect (&result, VECTOR, &tmpScalar); \ + switch (result) { \ + case PS_BINARY_DISECT_PASS: \ + break; \ + case PS_BINARY_DISECT_OUTSIDE_RANGE: \ + psTrace ("psLib.math", 6, "selected bin outside range"); \ + if (USE_END == -1) { RESULT = 0; } \ + if (USE_END == +1) { RESULT = VECTOR->n - 1; } \ + break; \ + case PS_BINARY_DISECT_INVALID_INPUT: \ + case PS_BINARY_DISECT_INVALID_TYPE: \ + psAbort ("programming error"); \ + break; \ } } @@ -92,20 +92,20 @@ float dX, dY, Xo, Yo, Xt; \ if (BIN == BOUNDS->n - 1) { \ - dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \ - dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \ - Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ - Yo = VECTOR->data.F32[BIN]; \ + dX = 0.5*(BOUNDS->data.F32[BIN+1] - BOUNDS->data.F32[BIN-1]); \ + dY = VECTOR->data.F32[BIN] - VECTOR->data.F32[BIN-1]; \ + Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ + Yo = VECTOR->data.F32[BIN]; \ } else { \ - dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \ - dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \ - Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ - Yo = VECTOR->data.F32[BIN]; \ - } \ - if (dY != 0) { \ - Xt = (VALUE - Yo)*dX/dY + Xo; \ - } else { \ - Xt = Xo; \ - } \ - Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \ + dX = 0.5*(BOUNDS->data.F32[BIN+2] - BOUNDS->data.F32[BIN]); \ + dY = VECTOR->data.F32[BIN+1] - VECTOR->data.F32[BIN]; \ + Xo = 0.5*(BOUNDS->data.F32[BIN+1] + BOUNDS->data.F32[BIN]); \ + Yo = VECTOR->data.F32[BIN]; \ + } \ + if (dY != 0) { \ + Xt = (VALUE - Yo)*dX/dY + Xo; \ + } else { \ + Xt = Xo; \ + } \ + Xt = PS_MIN (BOUNDS->data.F32[BIN+1], PS_MAX(BOUNDS->data.F32[BIN], Xt)); \ psTrace("psLib.math", 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \ Xo, Yo, dX, dY, Xt, VALUE); \ @@ -166,8 +166,8 @@ *****************************************************************************/ static bool vectorSampleMean(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -234,8 +234,8 @@ *****************************************************************************/ static long vectorMinMax(const psVector* myVector, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats - ) + const psVector* maskVector, + psMaskType maskVal, + psStats* stats + ) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -377,8 +377,8 @@ *****************************************************************************/ static bool vectorSampleStdev(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psMaskType maskVal, - psStats* stats) + const psVector* errors, + const psVector* maskVector, + psMaskType maskVal, + psStats* stats) { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); @@ -536,5 +536,5 @@ for (long j = 0; j < myVector->n; j++) { if (!tmpMask->data.U8[j] && - fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { + fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { tmpMask->data.U8[j] = 0xff; psTrace("psLib.math", 10, "Clipped %ld: %f +/- %f\n", j, @@ -547,5 +547,5 @@ for (long j = 0; j < myVector->n; j++) { if (!tmpMask->data.U8[j] && - fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { + fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { tmpMask->data.U8[j] = 0xff; psTrace("psLib.math", 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]); @@ -667,5 +667,5 @@ // Data range calculation failed psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); - goto escape; + goto escape; } psTrace("psLib.math", 6, "Data min/max is (%.2f, %.2f)\n", min, max); @@ -731,5 +731,5 @@ psTrace("psLib.math", 6, "Total data points is %ld\n", totalDataPoints); - PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0); + PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0); psTrace("psLib.math", 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]); @@ -737,18 +737,18 @@ // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median. stats->robustMedian = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, - binMedian, totalDataPoints/2.0); + 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"); - goto escape; + goto escape; } psTrace("psLib.math", 6, "Current robust median is %f\n", stats->robustMedian); // ADD step 4: Find the bins which contains the 15.8655% (-1 sigma) and 84.1345% (+1 sigma) data points - // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points - PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, -1); - PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, +1); - - PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, -1); - PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, +1); + // also find the 30.8538% (-0.5 sigma) and 69.1462% (+0.5 sigma) points + PS_BIN_FOR_VALUE(binLo, cumulative->nums, totalDataPoints * 0.158655f, -1); + PS_BIN_FOR_VALUE(binHi, cumulative->nums, totalDataPoints * 0.841345f, +1); + + PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, -1); + PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, +1); psTrace("psLib.math", 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n", @@ -759,5 +759,5 @@ if ((binLo < 0) || (binHi < 0)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); - goto escape; + goto escape; } @@ -769,11 +769,11 @@ // find the +0.5 and -0.5 sigma points with linear interpolation. binLo and binHi are the bins - // containing the value of interest. constrain the answer to be within the current bin. - // (extrapolation should not be needed and will result in errors) + // containing the value of interest. constrain the answer to be within the current bin. + // (extrapolation should not be needed and will result in errors) float binLoF32, binHiF32, binL2F32, binH2F32; - PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f); - PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f); - PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f); - PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f); + PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.158655f); + PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.841345f); + PS_BIN_INTERPOLATE (binL2F32, cumulative->nums, cumulative->bounds, binL2, totalDataPoints * 0.308538f); + PS_BIN_INTERPOLATE (binH2F32, cumulative->nums, cumulative->bounds, binH2, totalDataPoints * 0.691462f); // report +/- 1 sigma points @@ -810,8 +810,8 @@ // Free the histograms; they will be recreated on the next iteration, with new bounds psFree(histogram); - histogram = NULL; + histogram = NULL; psFree(cumulative); - cumulative = NULL; + cumulative = NULL; } else { // We've got the bin size correct now @@ -823,5 +823,5 @@ // XXX test lines while studying algorithm errors // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld : %f %f %f\n", - // stats->robustMedian, stats->robustStdev, + // stats->robustMedian, stats->robustStdev, // binLo, binL2, binMedian, binH2, binHi, binSize, max, min); @@ -835,11 +835,11 @@ // positions. psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, - binLo25, totalDataPoints * 0.25f); + binLo25, totalDataPoints * 0.25f); psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(cumulative->bounds, cumulative->nums, - binHi25, totalDataPoints * 0.75f); + binHi25, totalDataPoints * 0.75f); if (isnan(binLo25F32) || isnan(binHi25F32)) { psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n"); - goto escape; + goto escape; } @@ -851,5 +851,5 @@ for (long i = 0 ; i < myVector->n ; i++) { if (!mask->data.U8[i] && - (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { + (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { N50++; } @@ -986,8 +986,8 @@ } - // select the min and max bins, saturating on the lower and upper end-points - long binMin, binMax; - PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); // Generate the variables that will be used in the Gaussian fitting @@ -1164,8 +1164,8 @@ } - // select the min and max bins, saturating on the lower and upper end-points - long binMin, binMax; - PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); // Generate the variables that will be used in the Gaussian fitting @@ -1357,8 +1357,8 @@ } - // select the min and max bins, saturating on the lower and upper end-points - long binMin, binMax; - PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); + // select the min and max bins, saturating on the lower and upper end-points + long binMin, binMax; + PS_BIN_FOR_VALUE (binMin, histogram->bounds, guessMean - minFitSigma*guessStdev, -1); + PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, +1); // search for mode (peak of histogram within range mean-2sigma - mean+2sigma @@ -1494,4 +1494,5 @@ // calculate upper mean & stdev from parabolic fit -- use this as the result +#ifndef PS_NO_TRACE float upperStdev = sqrt(-0.5/poly->coeff[2]); float upperMean = poly->coeff[1]*PS_SQR(upperStdev); @@ -1499,4 +1500,5 @@ psTrace("psLib.math", 6, "The upper mean is %f.\n", upperMean); psTrace("psLib.math", 6, "The upper stdev is %f.\n", upperStdev); +#endif psFree (poly); @@ -2043,7 +2045,7 @@ *****************************************************************************/ static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal + psVector *yVec, + psS32 binNum, + psF32 yVal ) { @@ -2084,5 +2086,5 @@ // if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || - ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { + ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2094,5 +2096,5 @@ // 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]))) { + ((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");