Index: /branches/2dbias/ippconfig/gpc1/format_20100723.config =================================================================== --- /branches/2dbias/ippconfig/gpc1/format_20100723.config (revision 42718) +++ /branches/2dbias/ippconfig/gpc1/format_20100723.config (revision 42719) @@ -166,5 +166,6 @@ # 2D bias subtraction requires ppImage.config OVERSCAN.2D to be TRUE CELL.BIASSEC.SOURCE STR VALUE - CELL.BIASSEC STR [591:624,1:598],[1:624,599:608] +# exclude the last row and column from the bias region + CELL.BIASSEC STR [591:623,1:598],[1:623,599:607] # # a single region will be used for constant and 1D bias subtraction Index: /branches/2dbias/ippconfig/gpc1/ppImage.config =================================================================== --- /branches/2dbias/ippconfig/gpc1/ppImage.config (revision 42718) +++ /branches/2dbias/ippconfig/gpc1/ppImage.config (revision 42719) @@ -6,10 +6,13 @@ BIN2.YBIN S32 150 -OVERSCAN.SINGLE BOOL FALSE # Reduce overscan to a single value? -OVERSCAN.STAT STR MEDIAN -OVERSCAN.BOXCAR S32 3 +OVERSCAN.SINGLE BOOL FALSE # Reduce overscan to a single value? +OVERSCAN.STAT STR MEAN # MEAN | MEDIAN +OVERSCAN.BOXCAR S32 10 # Boxcar smoothing radius +OVERSCAN.2D BOOL FALSE # use a 2D model for the overscan subtraction? +OVERSCAN.2D.STAT STR MEAN # MEAN | MEDIAN +OVERSCAN.2D.BOXCAR S32 30 # Boxcar smoothing radius + OVERSCAN.2D.SUBSET METADATA - XY00 BOOL TRUE XY01 BOOL TRUE XY02 BOOL TRUE @@ -74,5 +77,4 @@ XY75 BOOL TRUE XY76 BOOL TRUE - XY77 BOOL TRUE END Index: /branches/2dbias/psLib/src/math/psStats.c =================================================================== --- /branches/2dbias/psLib/src/math/psStats.c (revision 42718) +++ /branches/2dbias/psLib/src/math/psStats.c (revision 42719) @@ -19,5 +19,4 @@ */ - // Note: choice of return values in many places is quite grey. For example, calculating the sample standard // deviation: here it's an error to get an array of size zero, but it's not an error to get an array of size @@ -63,10 +62,8 @@ #include "psString.h" - - /*****************************************************************************/ /* DEFINE STATEMENTS */ /*****************************************************************************/ -#define PS_GAUSS_WIDTH 3 // The width of the Gaussian smoothing. +#define PS_GAUSS_WIDTH 3 // The width of the Gaussian smoothing. #define PS_CLIPPED_NUM_ITER_LB 1 // This corresponds to N in the ADD. #define PS_CLIPPED_NUM_ITER_UB 10 @@ -77,66 +74,86 @@ #define TRACE "psLib.math" -#define MASK_MARK 0x80 // XXX : can we change this? bit to use internally to mark data as bad -#define PS_ROBUST_MAX_ITERATIONS 20 // Maximum number of iterations for robust statistics +#define MASK_MARK 0x80 // XXX : can we change this? bit to use internally to mark data as bad +#define PS_ROBUST_MAX_ITERATIONS 20 // Maximum number of iterations for robust statistics #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \ -(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) + (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, // 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(TRACE, 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; \ - } } - -# define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) { \ - 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]; \ - } 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)); \ +#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(TRACE, 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; \ + } \ + } + +#define PS_BIN_INTERPOLATE(RESULT, VECTOR, BOUNDS, BIN, VALUE) \ + { \ + 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]; \ + } \ + 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)); \ psTrace(TRACE, 6, "(Xo, Yo, dX, dY, Xt, Yt) is (%.2f %.2f %.2f %.2f %.2f %.2f)\n", \ - Xo, Yo, dX, dY, Xt, VALUE); \ - RESULT = Xt; } - -# define COUNT_WARNING(LIMIT, INTERVAL, ...) { \ - static int nCalls = 1; \ - if (nCalls < LIMIT) { \ - psWarning(__VA_ARGS__); \ - } \ - if (!(nCalls % INTERVAL)) { \ - psWarning(__VA_ARGS__); \ + Xo, Yo, dX, dY, Xt, VALUE); \ + RESULT = Xt; \ + } + +#define COUNT_WARNING(LIMIT, INTERVAL, ...) \ + { \ + static int nCalls = 1; \ + if (nCalls < LIMIT) \ + { \ + psWarning(__VA_ARGS__); \ + } \ + if (!(nCalls % INTERVAL)) \ + { \ + psWarning(__VA_ARGS__); \ psWarning("(warning raised %d times)", nCalls); \ - } \ - nCalls ++; \ -} + } \ + nCalls++; \ + } // Debug information @@ -195,24 +212,25 @@ To optmize this, use a macro and ifdef in or out the three states (errors, mask, range) *****************************************************************************/ - static bool vectorSampleMean(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psVectorMaskType maskVal, - psStats* stats) +static bool vectorSampleMean(const psVector *myVector, + const psVector *errors, + const psVector *maskVector, + psVectorMaskType maskVal, + psStats *stats) { - long count = 0; // Number of points contributing to this mean - psF32 mean = 0.0; // The mean + long count = 0; // Number of points contributing to this mean + psF32 mean = 0.0; // The mean psF32 weight; - psF32 *data = myVector->data.F32; // Dereference - int numData = myVector->n; // Number of data points + psF32 *data = myVector->data.F32; // Dereference + int numData = myVector->n; // Number of data points psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; bool useRange = stats->options & PS_STAT_USE_RANGE; - psF32 sumWeights = 0.0; // The sum of the weights + psF32 sumWeights = 0.0; // The sum of the weights psF32 *errorsData = (errors == NULL) ? NULL : errors->data.F32; - for (long i = 0; i < numData; i++) { + for (long i = 0; i < numData; i++) + { // Check if the data is with the specified range if (!isfinite(data[i])) @@ -224,22 +242,28 @@ if (maskData && (maskData[i] & maskVal)) continue; - if (errors) { + if (errors) + { weight = (errorsData[i] == 0) ? 0.0 : PS_SQR(errorsData[i]); mean += data[i] * weight; sumWeights += weight; - } else { + } + else + { mean += data[i]; } count++; - - } - if (errors) { + } + if (errors) + { mean = (count > 0) ? mean / sumWeights : NAN; - } else { + } + else + { mean = (count > 0) ? mean / count : NAN; } stats->sampleMean = mean; - if (!isnan(mean)) { + if (!isnan(mean)) + { stats->results |= PS_STAT_SAMPLE_MEAN; } @@ -260,9 +284,8 @@ (mask: 1, range: 0): 0.200 sec 0.244 sec *****************************************************************************/ - static long vectorMinMax(const psVector* myVector, - const psVector* maskVector, - psVectorMaskType maskVal, - psStats* stats - ) +static long vectorMinMax(const psVector *myVector, + const psVector *maskVector, + psVectorMaskType maskVal, + psStats *stats) { psF32 max = -PS_MAX_F32; // The calculated maximum @@ -270,11 +293,12 @@ psF32 *vector = myVector->data.F32; // Dereference the vector - int num = myVector->n; // Number of values - int numValid = 0; // Number of valid values + int num = myVector->n; // Number of values + int numValid = 0; // Number of valid values psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; bool useRange = stats->options & PS_STAT_USE_RANGE; - for (long i = 0; i < num; i++) { + for (long i = 0; i < num; i++) + { // Check if the data is with the specified range if (!isfinite(vector[i])) @@ -288,13 +312,16 @@ numValid++; - max = PS_MAX (vector[i], max); - min = PS_MIN (vector[i], min); + max = PS_MAX(vector[i], max); + min = PS_MIN(vector[i], min); } // XXX save numValid in psStats? - if (numValid == 0) { + if (numValid == 0) + { stats->max = NAN; stats->min = NAN; - } else { + } + else + { stats->max = max; stats->min = min; @@ -310,34 +337,40 @@ were no valid input vector elements). Expects F32 vector for input. *****************************************************************************/ -static bool vectorSampleMedian(const psVector* inVector, - const psVector* maskVector, +static bool vectorSampleMedian(const psVector *inVector, + const psVector *maskVector, psVectorMaskType maskVal, - psStats* stats) + psStats *stats) { bool useRange = stats->options & PS_STAT_USE_RANGE; psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; // Dereference the vector - psF32 *input = inVector->data.F32; // Dereference the vector + psF32 *input = inVector->data.F32; // Dereference the vector // use the temporary vector for the sorted output - stats->tmpData = psVectorRecycle (stats->tmpData, inVector->n, PS_TYPE_F32); + stats->tmpData = psVectorRecycle(stats->tmpData, inVector->n, PS_TYPE_F32); psVector *outVector = stats->tmpData; psF32 *output = outVector->data.F32; // Dereference the vector - if (maskData) psAssert (maskVector->n == inVector->n, "oops"); - - long count = 0; // Number of valid entries + if (maskData) + psAssert(maskVector->n == inVector->n, "oops"); + + long count = 0; // Number of valid entries // Store all non-masked data points within the min/max range // into the temporary vectors. - for (long i = 0; i < inVector->n; i++) { - psAssert (count >= 0, "oops"); - psAssert (count < outVector->n, "oops"); - psAssert (i >= 0, "oops"); - psAssert (i < inVector->n, "oops"); - - if (!isfinite(input[i])) continue; - if (useRange && (input[i] < stats->min)) continue; - if (useRange && (input[i] > stats->max)) continue; - if (maskData && (maskData[i] & maskVal)) continue; + for (long i = 0; i < inVector->n; i++) + { + psAssert(count >= 0, "oops"); + psAssert(count < outVector->n, "oops"); + psAssert(i >= 0, "oops"); + psAssert(i < inVector->n, "oops"); + + if (!isfinite(input[i])) + continue; + if (useRange && (input[i] < stats->min)) + continue; + if (useRange && (input[i] > stats->max)) + continue; + if (maskData && (maskData[i] & maskVal)) + continue; output[count] = input[i]; @@ -346,5 +379,6 @@ outVector->n = count; - if (count == 0) { + if (count == 0) + { COUNT_WARNING(10, 100, "No valid data in input vector.\n"); stats->sampleUQ = NAN; @@ -355,7 +389,8 @@ // Sort the temporary vector. - if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK) + if (!psVectorSort(outVector, outVector)) + { // Sort in-place (since it's a copy, it's OK) // an error in psVectorSort is a serious error: - // NULL input vector, psVectorCopy failure, invalid vector type + // NULL input vector, psVectorCopy failure, invalid vector type psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data.")); stats->sampleUQ = NAN; @@ -366,26 +401,29 @@ // Calculate the median. Use the average if the number of samples if even. - int midPt = (count/2); - psAssert (midPt >= 0, "oops"); - psAssert (midPt < outVector->n, "oops"); - if (count % 2 == 0) { - psAssert ((midPt - 1) >= 0, "oops"); - psAssert ((midPt - 1) < outVector->n, "oops"); + int midPt = (count / 2); + psAssert(midPt >= 0, "oops"); + psAssert(midPt < outVector->n, "oops"); + if (count % 2 == 0) + { + psAssert((midPt - 1) >= 0, "oops"); + psAssert((midPt - 1) < outVector->n, "oops"); stats->sampleMedian = 0.5 * (output[midPt - 1] + output[midPt]); - } else { + } + else + { stats->sampleMedian = output[midPt]; } - int Qmin = (int)(0.25*count); - int Qmax = (int)(0.75*count); - psAssert (Qmin >= 0, "oops"); - psAssert (Qmin < outVector->n, "oops"); - psAssert (Qmax >= 0, "oops"); - psAssert (Qmax < outVector->n, "oops"); + int Qmin = (int)(0.25 * count); + int Qmax = (int)(0.75 * count); + psAssert(Qmin >= 0, "oops"); + psAssert(Qmin < outVector->n, "oops"); + psAssert(Qmax >= 0, "oops"); + psAssert(Qmax < outVector->n, "oops"); // Calculate the quartile points exactly. stats->sampleUQ = output[Qmax]; stats->sampleLQ = output[Qmin]; - + stats->results |= PS_STAT_SAMPLE_MEDIAN; stats->results |= PS_STAT_SAMPLE_QUARTILE; @@ -409,18 +447,20 @@ *****************************************************************************/ -static bool vectorSampleStdev(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, +static bool vectorSampleStdev(const psVector *myVector, + const psVector *errors, + const psVector *maskVector, psVectorMaskType maskVal, - psStats* stats) + psStats *stats) { // This procedure requires the mean. If it has not been already // calculated, then call vectorSampleMean() - if (!(stats->results & PS_STAT_SAMPLE_MEAN)) { + if (!(stats->results & PS_STAT_SAMPLE_MEAN)) + { vectorSampleMean(myVector, errors, maskVector, maskVal, stats); } // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (isnan(stats->sampleMean)) { + if (isnan(stats->sampleMean)) + { COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN."); stats->sampleStdev = NAN; @@ -428,5 +468,5 @@ } - psF32 *data = myVector->data.F32; // Dereference + psF32 *data = myVector->data.F32; // Dereference psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; bool useRange = stats->options & PS_STAT_USE_RANGE; @@ -434,20 +474,24 @@ // Accumulate the sums - double mean = stats->sampleMean; // The mean - double sumSquares = 0.0; // Sum of the squares - double sumDiffs = 0.0; // Sum of the differences - double errorDivisor = 0.0; // Division for errors - long count = 0; // Number of data points being used - for (long i = 0; i < myVector->n; i++) { + double mean = stats->sampleMean; // The mean + double sumSquares = 0.0; // Sum of the squares + double sumDiffs = 0.0; // Sum of the differences + double errorDivisor = 0.0; // Division for errors + long count = 0; // Number of data points being used + for (long i = 0; i < myVector->n; i++) + { // Check if the data is with the specified range if (!isfinite(data[i])) continue; - if (useRange && (data[i] < stats->min)) { + if (useRange && (data[i] < stats->min)) + { continue; } - if (useRange && (data[i] > stats->max)) { + if (useRange && (data[i] > stats->max)) + { continue; } - if (maskData && (maskData[i] & maskVal)) { + if (maskData && (maskData[i] & maskVal)) + { continue; } @@ -457,10 +501,12 @@ sumDiffs += diff; count++; - if (errors) { + if (errors) + { errorDivisor += 1.0 / PS_SQR(errorsData[i]); } } - if (count == 0) { + if (count == 0) + { // This is an ambiguous case: error or not? // It's not an empty array (that's been asserted on previously), but everything's been masked out. @@ -470,5 +516,6 @@ return true; } - if (count == 1) { + if (count == 1) + { stats->sampleStdev = 0.0; COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): only one valid psVector elements (%ld). Setting stats->sampleStdev = 0.0.\n", count); @@ -476,7 +523,10 @@ } - if (errors) { + if (errors) + { stats->sampleStdev = (1.0 / sqrtf(errorDivisor)); - } else { + } + else + { stats->sampleStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / (float)count)) / (float)(count - 1)); } @@ -486,50 +536,58 @@ } -static bool vectorSampleMoments(const psVector* myVector, - const psVector* maskVector, +static bool vectorSampleMoments(const psVector *myVector, + const psVector *maskVector, psVectorMaskType maskVal, - psStats* stats) + psStats *stats) { // This procedure requires the mean and standard deviation - if (!(stats->results & PS_STAT_SAMPLE_MEAN)) { + if (!(stats->results & PS_STAT_SAMPLE_MEAN)) + { vectorSampleMean(myVector, NULL, maskVector, maskVal, stats); } - if (isnan(stats->sampleMean)) { + if (isnan(stats->sampleMean)) + { COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample mean is NAN.\n"); goto SAMPLE_MOMENTS_BAD; } - if (!(stats->results & PS_STAT_SAMPLE_STDEV)) { + if (!(stats->results & PS_STAT_SAMPLE_STDEV)) + { vectorSampleStdev(myVector, NULL, maskVector, maskVal, stats); } - if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) { + if (isnan(stats->sampleStdev) || stats->sampleStdev == 0.0) + { COUNT_WARNING(10, 100, "WARNING: vectorSampleMoments(): sample stdev is NAN or 0.\n"); goto SAMPLE_MOMENTS_BAD; } - psF32 *data = myVector->data.F32; // Dereference + psF32 *data = myVector->data.F32; // Dereference psVectorMaskType *maskData = (maskVector == NULL) ? NULL : maskVector->data.PS_TYPE_VECTOR_MASK_DATA; bool useRange = stats->options & PS_STAT_USE_RANGE; // Accumulate the sums - double mean = stats->sampleMean; // The mean - double sum3 = 0.0; // Sum of the cubes of the differences - double sum4 = 0.0; // Sum of the fourth powers of the differences - long count = 0; // Number of data points being used - for (long i = 0; i < myVector->n; i++) { + double mean = stats->sampleMean; // The mean + double sum3 = 0.0; // Sum of the cubes of the differences + double sum4 = 0.0; // Sum of the fourth powers of the differences + long count = 0; // Number of data points being used + for (long i = 0; i < myVector->n; i++) + { // Check if the data is with the specified range if (!isfinite(data[i])) continue; - if (useRange && (data[i] < stats->min)) { + if (useRange && (data[i] < stats->min)) + { continue; } - if (useRange && (data[i] > stats->max)) { + if (useRange && (data[i] > stats->max)) + { continue; } - if (maskData && (maskData[i] & maskVal)) { + if (maskData && (maskData[i] & maskVal)) + { continue; } - double diff = data[i] - mean; // Difference from the mean - double temp; // Temporary variable for accumulating + double diff = data[i] - mean; // Difference from the mean + double temp; // Temporary variable for accumulating sum3 += temp = PS_SQR(diff); @@ -539,8 +597,8 @@ } - psAssert(count > 1, "impossible"); // It should be, because we have a mean and standard deviation - - double stdev = stats->sampleStdev; // Standard deviation - double variance = PS_SQR(stdev); // Variance + psAssert(count > 1, "impossible"); // It should be, because we have a mean and standard deviation + + double stdev = stats->sampleStdev; // Standard deviation + double variance = PS_SQR(stdev); // Variance // Formula for skewness and kurtosis from Numerical Recipes in C, p 613. @@ -552,5 +610,5 @@ return true; - SAMPLE_MOMENTS_BAD: +SAMPLE_MOMENTS_BAD: // stats->sampleStdev has already been set stats->sampleSkewness = NAN; @@ -571,10 +629,9 @@ true for success; false otherwise *****************************************************************************/ -static bool vectorClippedStats(const psVector* myVector, - const psVector* errors, - psVector* maskInput, +static bool vectorClippedStats(const psVector *myVector, + const psVector *errors, + psVector *maskInput, psVectorMaskType maskValInput, - psStats* stats - ) + psStats *stats) { // Ensure that stats->clipIter is within the proper range. @@ -590,4 +647,5 @@ // unless we succeed, these will have NAN values stats->clippedMean = NAN; + stats->clippedMedian = NAN; stats->clippedStdev = NAN; stats->clippedNvalues = 0; @@ -597,10 +655,13 @@ // use the temporary vector for local temporary mask - stats->tmpMask = psVectorRecycle (stats->tmpMask, myVector->n, PS_TYPE_VECTOR_MASK); + stats->tmpMask = psVectorRecycle(stats->tmpMask, myVector->n, PS_TYPE_VECTOR_MASK); psVector *tmpMask = stats->tmpMask; psVectorInit(tmpMask, 0); - if (maskInput) { - for (long i = 0; i < myVector->n; i++) { - if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) { + if (maskInput) + { + for (long i = 0; i < myVector->n; i++) + { + if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) + { tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal; } @@ -610,5 +671,6 @@ // 1. Compute the sample median, which we save for output vectorSampleMedian(myVector, tmpMask, maskVal, stats); - if (isnan(stats->sampleMedian)) { + if (isnan(stats->sampleMedian)) + { COUNT_WARNING(10, 100, "Call to vectorSampleMedian returned NAN\n"); return true; @@ -618,5 +680,6 @@ // 2. Compute the sample standard deviation, which we save for output vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats); - if (isnan(stats->sampleStdev)) { + if (isnan(stats->sampleStdev)) + { COUNT_WARNING(10, 100, "Call to vectorSampleStdev returned NAN\n"); return true; @@ -625,22 +688,26 @@ // 3. Use the sample median as the first estimator of the mean X. - psF32 clippedMean = stats->sampleMedian; - + psF32 clippedMedian = stats->sampleMedian; + // 4. Use the sample stdev as the first estimator of the mean stdev. psF32 clippedStdev = stats->sampleStdev; // 5. Repeat N (stats->clipIter) times: - long numClipped = 0; // Number of values clipped - bool clipped = true; // Have we clipped anything in this iteration - for (int iter = 0; iter < stats->clipIter && clipped; iter++) { + long numClipped = 0; // Number of values clipped + bool clipped = true; // Have we clipped anything in this iteration + for (int iter = 0; iter < stats->clipIter && clipped; iter++) + { clipped = false; psTrace(TRACE, 6, "------------ Iteration %d ------------\n", iter); // a) Exclude all values x_i for which |x_i - x| > K * stdev - if (errors) { + if (errors) + { // XXXX if we convert errors to variance, this should square the other terms (A*A faster than // sqrt(A)) - for (long j = 0; j < myVector->n; j++) { + for (long j = 0; j < myVector->n; j++) + { if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] && - fabsf(myVector->data.F32[j] - clippedMean) > stats->clipSigma * errors->data.F32[j]) { + fabsf(myVector->data.F32[j] - clippedMedian) > stats->clipSigma * errors->data.F32[j]) + { tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; psTrace(TRACE, 10, "Clipped %ld: %f +/- %f\n", j, @@ -650,8 +717,12 @@ } } - } else { - for (long j = 0; j < myVector->n; j++) { + } + else + { + for (long j = 0; j < myVector->n; j++) + { if (!tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] && - fabsf(myVector->data.F32[j] - clippedMean) > (stats->clipSigma * clippedStdev)) { + fabsf(myVector->data.F32[j] - clippedMedian) > (stats->clipSigma * clippedStdev)) + { tmpMask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0xff; psTrace(TRACE, 10, "Clipped %ld: %f\n", j, myVector->data.F32[j]); @@ -662,25 +733,28 @@ } - // b) compute new mean and stdev + // b) compute new median and stdev // Allocate a psStats structure for calculating the mean and stdev. - // XXX Can we just use this psStats structure to calculate the SAMPLE MEAN and STDEV? - // psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - // vectorSampleMean(myVector, errors, tmpMask, maskVal, statsTmp); + // XXX Can we just use this psStats structure to calculate the SAMPLE MEDIAN and STDEV? + // psStats *statsTmp = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); + // vectorSampleMedian(myVector, errors, tmpMask, maskVal, statsTmp); // vectorSampleStdev(myVector, errors, tmpMask, maskVal, statsTmp); - vectorSampleMean(myVector, errors, tmpMask, maskVal, stats); + vectorSampleMedian(myVector, tmpMask, maskVal, stats); vectorSampleStdev(myVector, errors, tmpMask, maskVal, stats); - psTrace(TRACE, 6, "The new sample mean is %f\n", stats->sampleMean); + psTrace(TRACE, 6, "The new sample mean is %f\n", stats->sampleMedian); psTrace(TRACE, 6, "The new sample stdev is %f\n", stats->sampleStdev); - // If the new mean and stdev are NAN, we must exit the loop. + // If the new median and stdev are NAN, we must exit the loop. // Otherwise, use the new results and continue. - if (isnan(stats->sampleMean) || isnan(stats->sampleStdev)) { + if (isnan(stats->sampleMedian) || isnan(stats->sampleStdev)) + { iter = stats->clipIter; - COUNT_WARNING(10, 100, "vectorSampleMean() or vectorSampleStdev() returned a NAN.\n"); - clippedMean = NAN; + COUNT_WARNING(10, 100, "vectorSampleMedian() or vectorSampleStdev() returned a NAN.\n"); + clippedMedian = NAN; clippedStdev = NAN; return true; - } else { - clippedMean = stats->sampleMean; + } + else + { + clippedMedian = stats->sampleMedian; clippedStdev = stats->sampleStdev; } @@ -690,15 +764,31 @@ // Number of values used in calculation is the total number of data values, minus those we clipped stats->clippedNvalues = myVector->n - numClipped; - - // 7. The last calcuated value of x is the clipped mean. + // calculate the clipped mean + psF32 clippedMean = stats->sampleMean; + vectorSampleMean(myVector, errors, tmpMask, maskVal, stats); + if (isnan(stats->sampleMean)) + { + COUNT_WARNING(10, 100, "vectorSampleMean() returned a NAN.\n"); + clippedMean = NAN; + return true; + } + else + { + clippedMean = stats->sampleMean; + } + + // 7. The last calcuated value of x is the clipped median. // 8. The last calcuated value of stdev is the clipped stdev. - // we always return both stats even if only one was requested + // we always return all stats even if only one was requested stats->clippedMean = clippedMean; + stats->clippedMedian = clippedMedian; stats->clippedStdev = clippedStdev; stats->results |= PS_STAT_CLIPPED_MEAN; + stats->results |= PS_STAT_CLIPPED_MEDIAN; stats->results |= PS_STAT_CLIPPED_STDEV; psTrace(TRACE, 6, "The final clipped mean is %f\n", clippedMean); + psTrace(TRACE, 6, "The final clipped median is %f\n", clippedMedian); psTrace(TRACE, 6, "The final clipped stdev is %f\n", clippedStdev); @@ -725,11 +815,12 @@ *****************************************************************************/ #define INITIAL_NUM_BINS 1000.0 -static bool vectorRobustStats(const psVector* myVector, - const psVector* errors, - psVector* maskInput, +static bool vectorRobustStats(const psVector *myVector, + const psVector *errors, + psVector *maskInput, psVectorMaskType maskValInput, - psStats* stats) + psStats *stats) { - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(myVector); } @@ -741,7 +832,10 @@ psVector *mask = psVectorAlloc(myVector->n, PS_TYPE_VECTOR_MASK); // The actual mask we will use psVectorInit(mask, 0); - if (maskInput) { - for (long i = 0; i < myVector->n; i++) { - if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) { + if (maskInput) + { + for (long i = 0; i < myVector->n; i++) + { + if (maskInput->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskValInput) + { mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal; } @@ -751,11 +845,11 @@ // statsMinMax is only applied to a subset of the data points psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max - psHistogram *histogram = NULL; // Histogram of the data - psHistogram *cumulative = NULL; // Cumulative histogram of the data - float min = NAN, max = NAN; // Mimimum and maximum values - float sigma = NAN; // The robust standard deviation - long totalDataPoints = 0; // Total number of (unmasked) data points - - float binSize = 0.0; // Size of bins for histogram + psHistogram *histogram = NULL; // Histogram of the data + psHistogram *cumulative = NULL; // Cumulative histogram of the data + float min = NAN, max = NAN; // Mimimum and maximum values + float sigma = NAN; // The robust standard deviation + long totalDataPoints = 0; // Total number of (unmasked) data points + + float binSize = 0.0; // Size of bins for histogram long binLo, binHi; long binL2, binH2; @@ -764,13 +858,15 @@ // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop. - for (int iterate = 1; iterate > 0; iterate++) { + for (int iterate = 1; iterate > 0; iterate++) + { psTrace(TRACE, 6, "-------------------- Iterating on Bin size. Iteration number %d --------------------\n", iterate); - if (iterate >= PS_ROBUST_MAX_ITERATIONS) { - // This occurs when a large number of the values are identical --- a bin size cannot be found - // that will spread out the distribution. Therefore, set what we can, and fall over - // gracefully. - COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS); - goto escape; + if (iterate >= PS_ROBUST_MAX_ITERATIONS) + { + // This occurs when a large number of the values are identical --- a bin size cannot be found + // that will spread out the distribution. Therefore, set what we can, and fall over + // gracefully. + COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS); + goto escape; } @@ -779,11 +875,13 @@ min = statsMinMax->min; max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { + if (numValid == 0 || isnan(min) || isnan(max)) + { // Data range calculation failed COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); goto escape; } - if (!isfinite(max - min)) { - COUNT_WARNING(10, 100, "Range of of the input vector is too large: %lf.\n", (double)max - (double) min); + if (!isfinite(max - min)) + { + COUNT_WARNING(10, 100, "Range of of the input vector is too large: %lf.\n", (double)max - (double)min); goto escape; } @@ -791,5 +889,6 @@ // If all data points have the same value, then we set the appropriate members of stats and return. - if (fabs(max - min) <= FLT_EPSILON) { + if (fabs(max - min) <= FLT_EPSILON) + { stats->robustMedian = min; stats->robustStdev = 0.0; @@ -807,11 +906,14 @@ } - if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) { + if ((iterate == 1) && (stats->options & PS_STAT_USE_BINSIZE)) + { // Set initial bin size to the specified value. binSize = stats->binsize; psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { + } + else + { // Determine the bin size of the robust histogram, using the pre-defined number of bins - binSize = (max - min) / INITIAL_NUM_BINS; + binSize = (max - min) / INITIAL_NUM_BINS; } psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize); @@ -830,9 +932,10 @@ // Assert here so we can get more information about what is going wrong. psAssert(numBins > 0, "Invalid numBins: %ld max: %f min: %f binSize: %f", numBins, max, min, binSize); - + // Generate the histogram - histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins); + histogram = psHistogramAlloc(min - 2.0 * binSize, max + 2.0 * binSize, numBins); // XXXXX we need to consider this step if errors -> variance - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { + if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) + { // if psVectorHistogram returns false, we have a programming error psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n"); @@ -843,5 +946,6 @@ return false; } - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(histogram->bounds); PS_VECTOR_PRINT_F32(histogram->nums); @@ -854,39 +958,46 @@ int nMaxBin = histogram->nums->data.F32[0]; int iMaxBin = 0; - for (long i = 1; i < histogram->nums->n; i++) { - if (histogram->nums->data.F32[i] > nMaxBin) { + for (long i = 1; i < histogram->nums->n; i++) + { + if (histogram->nums->data.F32[i] > nMaxBin) + { nMaxBin = histogram->nums->data.F32[i]; iMaxBin = i; } } - if (nMaxBin > numValid / 2) { - float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize; - float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize; + if (nMaxBin > numValid / 2) + { + float minKeep = histogram->bounds->data.F32[iMaxBin] - 10 * binSize; + float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10 * binSize; int nInvalid = 0; - for (long i = 0; i < myVector->n; i++) { + for (long i = 0; i < myVector->n; i++) + { // skip the already-masked values - if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue; + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) + continue; bool invalid = false; invalid |= (myVector->data.F32[i] < minKeep); invalid |= (myVector->data.F32[i] > maxKeep); invalid |= (!isfinite(myVector->data.F32[i])); - if (!invalid) continue; + if (!invalid) + continue; mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal; - nInvalid ++; - } - - if (nInvalid) { - psTrace(TRACE, 6, "data is concentrated in a single bin (%d = %f - %f), masking %d extreme outliers and retrying\n", - iMaxBin, histogram->bounds->data.F32[iMaxBin], histogram->bounds->data.F32[iMaxBin+1], nInvalid); - psFree(histogram); - psFree(cumulative); - histogram = NULL; - cumulative = NULL; - continue; + nInvalid++; + } + + if (nInvalid) + { + psTrace(TRACE, 6, "data is concentrated in a single bin (%d = %f - %f), masking %d extreme outliers and retrying\n", + iMaxBin, histogram->bounds->data.F32[iMaxBin], histogram->bounds->data.F32[iMaxBin + 1], nInvalid); + psFree(histogram); + psFree(cumulative); + histogram = NULL; + cumulative = NULL; + continue; } // if we did not mask anything, give up. } - // We were causing psHistogramAlloc to assert. + // We were causing psHistogramAlloc to assert. // Assert here so we can get more information about what is going wrong. psAssert(numBins > 0, "Invalid numBins %ld max: %f min: %f binSize: %f", numBins, max, min, binSize); @@ -898,21 +1009,25 @@ cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1]; - // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins). - long Nc = 1; // track the current bin of cumulative - // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the - // upper end of the current histogram bin - for (long i = 1; i < histogram->nums->n - 1; i++) { - if (histogram->nums->data.F32[i] == 0.0) continue; - cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i]; - cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1]; - Nc ++; - } - long Nlast = Nc - 1; // last valid cumulative bin - for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled. - cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast]; - cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0; - } - - if (psTraceGetLevel("psLib.math") >= 8) { + // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins). + long Nc = 1; // track the current bin of cumulative + // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the + // upper end of the current histogram bin + for (long i = 1; i < histogram->nums->n - 1; i++) + { + if (histogram->nums->data.F32[i] == 0.0) + continue; + cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i]; + cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i + 1]; + Nc++; + } + long Nlast = Nc - 1; // last valid cumulative bin + for (long i = Nc; i < histogram->nums->n; i++) + { // Ensure the unused entries are filled. + cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast]; + cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i - 1] + 1.0; + } + + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(cumulative->bounds); PS_VECTOR_PRINT_F32(cumulative->nums); @@ -924,18 +1039,19 @@ // find bin which is the lower bound of median (value[binMedian] < median < value[binMedian+1] - PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0); - - psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian+1]); + PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints / 2.0, 0); + + psTrace(TRACE, 6, "The median bin is %ld (%.2f to %.2f)\n", binMedian, cumulative->bounds->data.F32[binMedian], cumulative->bounds->data.F32[binMedian + 1]); // 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); + // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian); - // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there. - // Instead, do a 5-point linear fit. - stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there. + // Instead, do a 5-point linear fit. + stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints / 2.0); // convert bin to bin value: this is the robust histogram median. - if (isnan(stats->robustMedian)) { + if (isnan(stats->robustMedian)) + { COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n"); goto escape; @@ -949,8 +1065,7 @@ PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0); PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0); - PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0); - PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0); - - + PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0); + PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0); + psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n", binLo, binHi); @@ -960,15 +1075,16 @@ psTrace(TRACE, 6, "binH2 midpoint is %f\n", PS_BIN_MIDPOINT(cumulative, binH2)); - if ((binLo < 0) || (binHi < 0)) { + if ((binLo < 0) || (binHi < 0)) + { COUNT_WARNING(10, 100, "Failed to calculate the 15.8655%% and 84.1345%% data points.\n"); goto escape; } - + // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma // positions. psTrace(TRACE, 6, "binLo is %ld. Nums at that bin and the next are (%.2f, %.2f)\n", - binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo+1]); + binLo, cumulative->nums->data.F32[binLo], cumulative->nums->data.F32[binLo + 1]); psTrace(TRACE, 6, "binHi is %ld. Nums at that bin and the next are (%.2f, %.2f)\n", - binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi+1]); + binHi, cumulative->nums->data.F32[binHi], cumulative->nums->data.F32[binHi + 1]); // find the +0.5 and -0.5 sigma points with linear interpolation. binLo and binHi are the bins @@ -977,24 +1093,24 @@ float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32; #if (0) - PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo, - totalDataPoints * 0.158655f); - PS_BIN_INTERPOLATE (binHiF32, cumulative->nums, cumulative->bounds, binHi, - 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 (binL4F32, cumulative->nums, cumulative->bounds, binL4, - totalDataPoints * 0.022750f); - PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, - totalDataPoints * 0.977250f); + PS_BIN_INTERPOLATE(binLoF32, cumulative->nums, cumulative->bounds, binLo, + totalDataPoints * 0.158655f); + PS_BIN_INTERPOLATE(binHiF32, cumulative->nums, cumulative->bounds, binHi, + 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(binL4F32, cumulative->nums, cumulative->bounds, binL4, + totalDataPoints * 0.022750f); + PS_BIN_INTERPOLATE(binH4F32, cumulative->nums, cumulative->bounds, binH4, + totalDataPoints * 0.977250f); #else binLoF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo, totalDataPoints * 0.158655); - binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345); + binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345); binL2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL2, totalDataPoints * 0.308538); - binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462); + binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462); binL4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL4, totalDataPoints * 0.022750); - binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250); -#endif + binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250); +#endif // report +/- 1 sigma points psTrace(TRACE, 5, @@ -1008,14 +1124,31 @@ binL4F32, binH4F32); - // If some of the fits failed, attempt to fix this - if (!isfinite(binLoF32) && isfinite(binHiF32)) { binLoF32 = -1.0 * binHiF32; } - if (!isfinite(binHiF32) && isfinite(binLoF32)) { binHiF32 = -1.0 * binLoF32; } - if (!isfinite(binL2F32) && isfinite(binH2F32)) { binL2F32 = -1.0 * binH2F32; } - if (!isfinite(binH2F32) && isfinite(binL2F32)) { binH2F32 = -1.0 * binL2F32; } - if (!isfinite(binL4F32) && isfinite(binH4F32)) { binL4F32 = -1.0 * binH4F32; } - if (!isfinite(binH4F32) && isfinite(binL4F32)) { binH4F32 = -1.0 * binL4F32; } - + // If some of the fits failed, attempt to fix this + if (!isfinite(binLoF32) && isfinite(binHiF32)) + { + binLoF32 = -1.0 * binHiF32; + } + if (!isfinite(binHiF32) && isfinite(binLoF32)) + { + binHiF32 = -1.0 * binLoF32; + } + if (!isfinite(binL2F32) && isfinite(binH2F32)) + { + binL2F32 = -1.0 * binH2F32; + } + if (!isfinite(binH2F32) && isfinite(binL2F32)) + { + binH2F32 = -1.0 * binL2F32; + } + if (!isfinite(binL4F32) && isfinite(binH4F32)) + { + binL4F32 = -1.0 * binH4F32; + } + if (!isfinite(binH4F32) && isfinite(binL4F32)) + { + binH4F32 = -1.0 * binL4F32; + } + // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma) - float sigma1 = (binH2F32 - binL2F32); @@ -1023,21 +1156,29 @@ float sigma4 = (binH4F32 - binL4F32) / 4.0; - // Fix again? - if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) { sigma1 = (sigma2 + sigma4) / 2.0; } - if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) { sigma2 = (sigma1 + sigma4) / 2.0; } - if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) { sigma4 = (sigma2 + sigma1) / 2.0; } - + // Fix again? + if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) + { + sigma1 = (sigma2 + sigma4) / 2.0; + } + if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) + { + sigma2 = (sigma1 + sigma4) / 2.0; + } + if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) + { + sigma4 = (sigma2 + sigma1) / 2.0; + } + // take the smallest of the three: if we have a clump with wide outliers, sigma2 and // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2 // will be biased high. - // sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4)); - // CZW: Instead, take the median. Taking the MIN forces a bias on unbiased data. - // It seems like occasionally getting the wrong answer on a complex distribution - // is more acceptable than always getting the wrong answer for simple ones. - - - sigma = PS_MAX( PS_MIN(sigma1,sigma2), - PS_MIN( PS_MAX(sigma1,sigma2), - sigma4)); + // sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4)); + // CZW: Instead, take the median. Taking the MIN forces a bias on unbiased data. + // It seems like occasionally getting the wrong answer on a complex distribution + // is more acceptable than always getting the wrong answer for simple ones. + + sigma = PS_MAX(PS_MIN(sigma1, sigma2), + PS_MIN(PS_MAX(sigma1, sigma2), + sigma4)); psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1); @@ -1046,57 +1187,57 @@ psTrace(TRACE, 6, "The current sigma is %f.\n", sigma); - // stats->robustStdev = sigma; - stats->robustStdev = sigma; - -#if (CZW && 0) - // Skewness check: Find least biased sample for each pair. - sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian, - stats->robustMedian - binL2F32); - sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian, - stats->robustMedian - binLoF32); - sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian, - stats->robustMedian - binL4F32); - // Kurtosis check: Take median sample as the solution. - stats->robustStdev = PS_MAX( PS_MIN(sigma1,sigma2), - PS_MIN( PS_MAX(sigma1,sigma2), - sigma4)); + // stats->robustStdev = sigma; + stats->robustStdev = sigma; + +#if (CZW && 0) + // Skewness check: Find least biased sample for each pair. + sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian, + stats->robustMedian - binL2F32); + sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian, + stats->robustMedian - binLoF32); + sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian, + stats->robustMedian - binL4F32); + // Kurtosis check: Take median sample as the solution. + stats->robustStdev = PS_MAX(PS_MIN(sigma1, sigma2), + PS_MIN(PS_MAX(sigma1, sigma2), + sigma4)); #endif - #if (CZW) - // printf("CZW: bad sigma?: %f %f %f %f %f %f %f %f %f %f\n", - // binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32, - // sigma1,sigma2,sigma4,sigma); - - printf("CZW Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n", - iterate, - stats->robustMedian,stats->robustStdev, - fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]), - - cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2], - cumulative->bounds->data.F32[binMedian-1], - cumulative->bounds->data.F32[binMedian], - cumulative->bounds->data.F32[binMedian+1], - cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3], - - cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2], - cumulative->nums->data.F32[binMedian-1], - cumulative->nums->data.F32[binMedian], - cumulative->nums->data.F32[binMedian+1], - cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]); - // PS_VECTOR_PRINT_F32(histogram->bounds); - // PS_VECTOR_PRINT_F32(histogram->nums); - // PS_VECTOR_PRINT_F32(cumulative->bounds); - // PS_VECTOR_PRINT_F32(cumulative->nums); + // printf("CZW: bad sigma?: %f %f %f %f %f %f %f %f %f %f\n", + // binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32, + // sigma1,sigma2,sigma4,sigma); + + printf("CZW Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n", + iterate, + stats->robustMedian, stats->robustStdev, + fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]), + + cumulative->bounds->data.F32[binMedian - 3], cumulative->bounds->data.F32[binMedian - 2], + cumulative->bounds->data.F32[binMedian - 1], + cumulative->bounds->data.F32[binMedian], + cumulative->bounds->data.F32[binMedian + 1], + cumulative->bounds->data.F32[binMedian + 2], cumulative->bounds->data.F32[binMedian + 3], + + cumulative->nums->data.F32[binMedian - 3], cumulative->nums->data.F32[binMedian - 2], + cumulative->nums->data.F32[binMedian - 1], + cumulative->nums->data.F32[binMedian], + cumulative->nums->data.F32[binMedian + 1], + cumulative->nums->data.F32[binMedian + 2], cumulative->nums->data.F32[binMedian + 3]); + // PS_VECTOR_PRINT_F32(histogram->bounds); + // PS_VECTOR_PRINT_F32(histogram->nums); + // PS_VECTOR_PRINT_F32(cumulative->bounds); + // PS_VECTOR_PRINT_F32(cumulative->nums); #endif // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more // than 25 bins from the median, recalculate the bin size, and perform the algorithm again. - if (sigma < (3.0 * binSize)) { + if (sigma < (3.0 * binSize)) + { psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); - // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins - psF32 medianLo = stats->robustMedian - 25*binSize; - psF32 medianHi = stats->robustMedian + 25*binSize; + // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins + psF32 medianLo = stats->robustMedian - 25 * binSize; + psF32 medianHi = stats->robustMedian + 25 * binSize; // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region @@ -1107,19 +1248,23 @@ // psTrace(TRACE, 6, "The median is at bin number %ld. We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi); psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi); - int Nmasked = 0; - for (long i = 0 ; i < myVector->n ; i++) { - if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) { - if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue; + int Nmasked = 0; + for (long i = 0; i < myVector->n; i++) + { + if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) + { + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) + continue; mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK; psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]); - Nmasked ++; + Nmasked++; } } - if (Nmasked == 0) { - // no significant change to the sigma & binsize -- we are done here - iterate = -1; - continue; - } + if (Nmasked == 0) + { + // no significant change to the sigma & binsize -- we are done here + iterate = -1; + continue; + } // Free the histograms; they will be recreated on the next iteration, with new bounds @@ -1130,5 +1275,6 @@ cumulative = NULL; - if (iterate >= PS_ROBUST_MAX_ITERATIONS) { + if (iterate >= PS_ROBUST_MAX_ITERATIONS) + { // This occurs when a large number of the values are identical --- a bin size cannot be found // that will spread out the distribution. Therefore, set what we can, and fall over @@ -1147,5 +1293,7 @@ return true; } - } else { + } + else + { // We've got the bin size correct now psTrace(TRACE, 6, "*************: No more iteration. sigma is %f\n", sigma); @@ -1153,5 +1301,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", @@ -1161,6 +1309,6 @@ // ADD step 7: Find the bins which contains the 25% and 75% data points. long binLo25, binHi25; - PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0); - PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0); + PS_BIN_FOR_VALUE(binLo25, cumulative->nums, totalDataPoints * 0.25f, 0); + PS_BIN_FOR_VALUE(binHi25, cumulative->nums, totalDataPoints * 0.75f, 0); psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25); @@ -1169,5 +1317,6 @@ psF32 binLo25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f); psF32 binHi25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f); - if (isnan(binLo25F32) || isnan(binHi25F32)) { + if (isnan(binLo25F32) || isnan(binHi25F32)) + { COUNT_WARNING(10, 100, "could not determine the robustUQ or LQ: fitLinearSearchForYThenReturnBin() returned a NAN.\n"); goto escape; @@ -1179,7 +1328,9 @@ binLo25F32, binHi25F32); long N50 = 0; - for (long i = 0 ; i < myVector->n ; i++) { + for (long i = 0; i < myVector->n; i++) + { if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[i] && - (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) { + (binLo25F32 <= myVector->data.F32[i]) && (binHi25F32 >= myVector->data.F32[i])) + { N50++; } @@ -1226,15 +1377,17 @@ * version follows the upper portion of the distribution until it passes 0.5*peak ********************/ -static bool vectorFittedStats (const psVector* myVector, - const psVector* errors, - psVector* mask, - psVectorMaskType maskVal, - psStats* stats) +static bool vectorFittedStats(const psVector *myVector, + const psVector *errors, + psVector *mask, + psVectorMaskType maskVal, + psStats *stats) { // This procedure requires the mean. If it has not been already // calculated, then call vectorSampleMean() - if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) { - if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) { + if (!(stats->results & PS_STAT_ROBUST_MEDIAN)) + { + if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); return false; @@ -1243,22 +1396,28 @@ // If the mean is NAN, then generate a warning and set the stdev to NAN. - if (isnan(stats->robustMedian)) { - stats->fittedMean = NAN; - stats->fittedStdev = NAN; - stats->results |= PS_STAT_FITTED_MEAN; - stats->results |= PS_STAT_FITTED_STDEV; + if (isnan(stats->robustMedian)) + { + stats->fittedMean = NAN; + stats->fittedStdev = NAN; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; } - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN; - stats->results |= PS_STAT_FITTED_STDEV; + if (stats->robustStdev <= FLT_EPSILON) + { + stats->fittedMean = stats->robustMedian; + stats->fittedStdev = stats->robustStdev; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; } - if (myVector->n < 1) { printf("There are no elements in this vector.\n"); abort(); } - float guessStdev = stats->robustStdev; // pass the guess sigma - float guessMean = stats->robustMedian; // pass the guess mean + if (myVector->n < 1) + { + printf("There are no elements in this vector.\n"); + abort(); + } + float guessStdev = stats->robustStdev; // pass the guess sigma + float guessMean = stats->robustMedian; // pass the guess mean psTrace(TRACE, 6, "The ** starting ** guess mean is %f.\n", guessMean); @@ -1266,27 +1425,31 @@ bool done = false; - for (int iteration = 0; !done && (iteration < 2); iteration ++) { + for (int iteration = 0; !done && (iteration < 2); iteration++) + { psStats *statsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); // Statistics for min and max psF32 binSize = 1; - if (stats->options & PS_STAT_USE_BINSIZE) { + if (stats->options & PS_STAT_USE_BINSIZE) + { // Set initial bin size to the specified value. binSize = stats->binsize; psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { + } + else + { // construct a histogram with (sigma/2 < binsize < sigma) // set roughly so that the lowest bins have about 2 cnts // Nsmallest ~ N50 / (4*dN)) - // psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); - - // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian. - // Furthermore, we only use the inner +/- 2 sigma parts. Therefore, define the - // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7). robustN50 - // contains half the total points, so 2 * robustN50 / 50 is the fraction of all - // points in the 2 sigma bin. Dance the erf() relations around, and it looks like - // there's a factor of about 1/20 to include. Keep the PS_MAX to ensure we never bin - // wider than 1 sigma when the number of points is small. - psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0)); - binSize = guessStdev / dN; + // psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); + + // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian. + // Furthermore, we only use the inner +/- 2 sigma parts. Therefore, define the + // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7). robustN50 + // contains half the total points, so 2 * robustN50 / 50 is the fraction of all + // points in the 2 sigma bin. Dance the erf() relations around, and it looks like + // there's a factor of about 1/20 to include. Keep the PS_MAX to ensure we never bin + // wider than 1 sigma when the number of points is small. + psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0)); + binSize = guessStdev / dN; } @@ -1295,5 +1458,6 @@ float min = statsMinMax->min; float max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { + if (numValid == 0 || isnan(min) || isnan(max)) + { COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); @@ -1302,5 +1466,6 @@ // If all data points have the same value, then we set the appropriate members of stats and return. - if (fabs(max - min) <= FLT_EPSILON) { + if (fabs(max - min) <= FLT_EPSILON) + { COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min); stats->fittedMean = min; @@ -1314,24 +1479,29 @@ // XXX can we calculate the binMin, binMax **before** building this histogram? // if the range is too absurd, adjust numBins & binSize - // We no longer want to reset the binSize here, as it can cause odd things. Better to select - // a number of bins, and then set the min/max values to put those bins sanely around the mean. - // long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); - // binSize = (max - min) / (float) numBins; + // We no longer want to reset the binSize here, as it can cause odd things. Better to select + // a number of bins, and then set the min/max values to put those bins sanely around the mean. + // long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); + // binSize = (max - min) / (float) numBins; psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); - // psTrace(TRACE, 6, "The numBins is %ld\n", numBins); - + // psTrace(TRACE, 6, "The numBins is %ld\n", numBins); #define FITTED_CLIPPING_NUM 5.0 - if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) { - min = guessMean - FITTED_CLIPPING_NUM * guessStdev; - } - if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) { - max = guessMean + FITTED_CLIPPING_NUM * guessStdev; - } - long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); - if (CZW) { printf("I've clipped: %f %f => %f %f ; %f %ld\n",guessMean,guessStdev,min,max,binSize,stats->robustN50); } + if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) + { + min = guessMean - FITTED_CLIPPING_NUM * guessStdev; + } + if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) + { + max = guessMean + FITTED_CLIPPING_NUM * guessStdev; + } + long numBins = PS_MAX(50, PS_MIN(10000, (max - min) / binSize)); + if (CZW) + { + printf("I've clipped: %f %f => %f %f ; %f %ld\n", guessMean, guessStdev, min, max, binSize, stats->robustN50); + } psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { + if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) + { COUNT_WARNING(10, 100, "Unable to generate histogram for fitted statistics v4.\n"); psFree(histogram); @@ -1339,5 +1509,6 @@ goto escape; } - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(histogram->nums); } @@ -1347,16 +1518,20 @@ // set the full-range upper and lower limits psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { + if (isfinite(stats->clipSigma)) + { maxFitSigma = fabs(stats->clipSigma); } - if (isfinite(stats->max)) { + if (isfinite(stats->max)) + { maxFitSigma = fabs(stats->max); } psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { + if (isfinite(stats->clipSigma)) + { minFitSigma = fabs(stats->clipSigma); } - if (isfinite(stats->min)) { + if (isfinite(stats->min)) + { minFitSigma = fabs(stats->min); } @@ -1364,8 +1539,9 @@ // 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, 0); - PS_BIN_FOR_VALUE (binMax, histogram->bounds, guessMean + maxFitSigma*guessStdev, 0); - - if (binMin == binMax) { + PS_BIN_FOR_VALUE(binMin, histogram->bounds, guessMean - minFitSigma * guessStdev, 0); + PS_BIN_FOR_VALUE(binMax, histogram->bounds, guessMean + maxFitSigma * guessStdev, 0); + + if (binMin == binMax) + { COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); @@ -1375,24 +1551,31 @@ // search for mode (peak of histogram within range mean-2sigma - mean+2sigma - long binPeak = binMin; + long binPeak = binMin; float valPeak = histogram->nums->data.F32[binPeak]; - for (int i = binMin; i < binMax; i++) { - if (histogram->nums->data.F32[i] > valPeak) { + for (int i = binMin; i < binMax; i++) + { + if (histogram->nums->data.F32[i] > valPeak) + { binPeak = i; valPeak = histogram->nums->data.F32[binPeak]; } - psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); - if (CZW) { printf("CENTERED_HISTOGRAM: %f %f\n", - PS_BIN_MIDPOINT(histogram,i), - histogram->nums->data.F32[i]); } - } - psTrace (TRACE, 6, "\n"); - - if (CZW) { printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n", - binMin,PS_BIN_MIDPOINT(histogram,binMin),histogram->nums->data.F32[binMin], - binMax,PS_BIN_MIDPOINT(histogram,binMax),histogram->nums->data.F32[binMax], - binPeak,PS_BIN_MIDPOINT(histogram,binPeak),histogram->nums->data.F32[binPeak]); - } - + psTrace(TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]); + if (CZW) + { + printf("CENTERED_HISTOGRAM: %f %f\n", + PS_BIN_MIDPOINT(histogram, i), + histogram->nums->data.F32[i]); + } + } + psTrace(TRACE, 6, "\n"); + + if (CZW) + { + printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n", + binMin, PS_BIN_MIDPOINT(histogram, binMin), histogram->nums->data.F32[binMin], + binMax, PS_BIN_MIDPOINT(histogram, binMax), histogram->nums->data.F32[binMax], + binPeak, PS_BIN_MIDPOINT(histogram, binPeak), histogram->nums->data.F32[binPeak]); + } + // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); @@ -1402,7 +1585,6 @@ psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); - - float lowfitMean = NAN; - float lowfitStdev = NAN; + float lowfitMean = NAN; + float lowfitStdev = NAN; { // fit the lower half of the distribution @@ -1410,13 +1592,17 @@ // run up until we drop below 0.50*valPeak long binS = binMin; - long binE = PS_MIN (binPeak + 3, binMax); - for (int i = binPeak - 3; i >= binMin; i--) { - if (histogram->nums->data.F32[i] < 0.25*valPeak) { + long binE = PS_MIN(binPeak + 3, binMax); + for (int i = binPeak - 3; i >= binMin; i--) + { + if (histogram->nums->data.F32[i] < 0.25 * valPeak) + { binS = i; break; } } - for (int i = binPeak + 3; i < binMax; i++) { - if (histogram->nums->data.F32[i] < 0.50*valPeak) { + for (int i = binPeak + 3; i < binMax; i++) + { + if (histogram->nums->data.F32[i] < 0.50 * valPeak) + { binE = i; break; @@ -1431,5 +1617,6 @@ psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates long j = 0; - for (long i = binS; i < binE; i++) { + for (long i = binS; i < binE; i++) + { if (histogram->nums->data.F32[i] <= 0.0) continue; @@ -1440,23 +1627,25 @@ } y->n = x->n = j; - + // 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 psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); + bool status = psVectorFitPolynomial1D(poly, NULL, 0, y, NULL, x); #if (CZW && 1) - printf("CZW: LowfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]); - for (long i = 0; i < x->n; i++) { - printf("CZW: Lowfit: %d %ld %f %f %f\n", - status,i,x->data.F32[i],y->data.F32[i], - poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + - poly->coeff[2] * pow(x->data.F32[i],2)); - } + printf("CZW: LowfitPoly: %f %f %f\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + for (long i = 0; i < x->n; i++) + { + printf("CZW: Lowfit: %d %ld %f %f %f\n", + status, i, x->data.F32[i], y->data.F32[i], + poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + + poly->coeff[2] * pow(x->data.F32[i], 2)); + } #endif psFree(x); psFree(y); - if (!status) { + if (!status) + { psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); @@ -1467,5 +1656,6 @@ } - if (poly->coeff[2] >= 0.0) { + if (poly->coeff[2] >= 0.0) + { COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); @@ -1477,6 +1667,7 @@ // tends to be found in a single bin. make one attempt to recover by dropping the guessStdev // down by a jump and trying again - if (iteration == 0) { - guessStdev = 0.25*guessStdev; + if (iteration == 0) + { + guessStdev = 0.25 * guessStdev; psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev); continue; @@ -1488,6 +1679,6 @@ // calculate lower mean & stdev from parabolic fit -- use this as the result - lowfitStdev = sqrt(-0.5/poly->coeff[2]); - lowfitMean = poly->coeff[1]*PS_SQR(lowfitStdev); + lowfitStdev = sqrt(-0.5 / poly->coeff[2]); + lowfitMean = poly->coeff[1] * PS_SQR(lowfitStdev); psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); @@ -1498,11 +1689,11 @@ } - float fullfitMean = NAN; - float fullfitStdev = NAN; - float minValueSym = NAN; - float maxValueSym = NAN; - - // try the full fit as well: - { + float fullfitMean = NAN; + float fullfitStdev = NAN; + float minValueSym = NAN; + float maxValueSym = NAN; + + // try the full fit as well: + { // fit a symmetric distribution // run up until we drop below 0.15*valPeak @@ -1510,12 +1701,16 @@ long binS = binMin; long binE = binMax; - for (int i = binPeak - 3; i >= binMin; i--) { - if (histogram->nums->data.F32[i] < 0.15*valPeak) { + for (int i = binPeak - 3; i >= binMin; i--) + { + if (histogram->nums->data.F32[i] < 0.15 * valPeak) + { binS = i; break; } } - for (int i = binPeak + 3; i < binMax; i++) { - if (histogram->nums->data.F32[i] < 0.15*valPeak) { + for (int i = binPeak + 3; i < binMax; i++) + { + if (histogram->nums->data.F32[i] < 0.15 * valPeak) + { binE = i; break; @@ -1531,5 +1726,6 @@ psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates long j = 0; - for (long i = binS; i < binE; i++) { + for (long i = binS; i < binE; i++) + { if (histogram->nums->data.F32[i] <= 0.0) continue; @@ -1543,18 +1739,20 @@ // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); + bool status = psVectorFitPolynomial1D(poly, NULL, 0, y, NULL, x); #if (CZW && 1) - printf("CZW: FullfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]); - for (long i = 0; i < x->n; i++) { - printf("CZW: Fullfit: %d %ld %f %f %f\n", - status,i,x->data.F32[i],y->data.F32[i], - poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + - poly->coeff[2] * pow(x->data.F32[i],2)); - } + printf("CZW: FullfitPoly: %f %f %f\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + for (long i = 0; i < x->n; i++) + { + printf("CZW: Fullfit: %d %ld %f %f %f\n", + status, i, x->data.F32[i], y->data.F32[i], + poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + + poly->coeff[2] * pow(x->data.F32[i], 2)); + } #endif psFree(x); psFree(y); - if (!status) { + if (!status) + { psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); @@ -1566,6 +1764,6 @@ // calculate upper mean & stdev from parabolic fit -- ignore this value - fullfitStdev = sqrt(-0.5/poly->coeff[2]); - fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev); + fullfitStdev = sqrt(-0.5 / poly->coeff[2]); + fullfitMean = poly->coeff[1] * PS_SQR(fullfitStdev); #ifndef PS_NO_TRACE @@ -1580,5 +1778,6 @@ // saturate on min or max value - if (fullfitMean < minValueSym) { + if (fullfitMean < minValueSym) + { fullfitMean = minValueSym; psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); @@ -1586,50 +1785,54 @@ // saturate on min or max value - if (fullfitMean > maxValueSym) { + if (fullfitMean > maxValueSym) + { fullfitMean = maxValueSym; psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); } - - psFree (poly); - } - - // we now have the fullfit and the lowfit mean and stdev values - // accept the fullfit unless minValueSym < lowfitMean < fullfitMean - - if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym)) { - guessMean = lowfitMean; - guessStdev = lowfitStdev; - } else { - guessMean = fullfitMean; - guessStdev = fullfitStdev; - } - - if (!isfinite(guessMean) || !isfinite(guessStdev)) { - guessMean = stats->robustMedian; - guessStdev = stats->robustStdev; - } - - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - - + psFree(poly); + } + + // we now have the fullfit and the lowfit mean and stdev values + // accept the fullfit unless minValueSym < lowfitMean < fullfitMean + + if (isfinite(lowfitMean) && isfinite(lowfitStdev) && (lowfitMean < fullfitMean) && (lowfitMean > minValueSym)) + { + guessMean = lowfitMean; + guessStdev = lowfitStdev; + } + else + { + guessMean = fullfitMean; + guessStdev = fullfitStdev; + } + + if (!isfinite(guessMean) || !isfinite(guessStdev)) + { + guessMean = stats->robustMedian; + guessStdev = stats->robustStdev; + } + + if (guessStdev > 0.75 * stats->robustStdev) + { + done = true; + } + #if (CZW && 1) - printf("CZW IN FITTED: iter %d %f \n" - " low %f %f \n" - " full %f %f \n" - " robust %f %f \n" - " final %f %f\n", - iteration,minValueSym, - lowfitMean,lowfitStdev, - fullfitMean,fullfitStdev, - stats->robustMedian,stats->robustStdev, - guessMean,guessStdev); + printf("CZW IN FITTED: iter %d %f \n" + " low %f %f \n" + " full %f %f \n" + " robust %f %f \n" + " final %f %f\n", + iteration, minValueSym, + lowfitMean, lowfitStdev, + fullfitMean, fullfitStdev, + stats->robustMedian, stats->robustStdev, + guessMean, guessStdev); #endif // Clean up after fitting - psFree (histogram); - psFree (statsMinMax); + psFree(histogram); + psFree(statsMinMax); } @@ -1655,5 +1858,4 @@ return true; } - /****************************************************************************** @@ -1668,17 +1870,19 @@ { psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); - psTrace(TRACE, 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int) histogram->nums->n, sigma); + psTrace(TRACE, 5, "(histogram->nums->n, sigma) is (%d, %.2f\n", (int)histogram->nums->n, sigma); PS_ASSERT_PTR_NON_NULL(histogram, NULL); PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(histogram->nums); } - long numBins = histogram->nums->n; // Number of histogram bins + long numBins = histogram->nums->n; // Number of histogram bins psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); // Smoothed version of histogram bins - const psVector *bounds = histogram->bounds; // The bounds for the histogram bins - - if (!histogram->uniform) { + const psVector *bounds = histogram->bounds; // The bounds for the histogram bins + + if (!histogram->uniform) + { // // We get here if the histogram is non-uniform. This code is not tested. @@ -1688,5 +1892,6 @@ "histograms has not been tested or used.\n"); - for (long i = 0; i < numBins; i++) { + for (long i = 0; i < numBins; i++) + { // Determine the midpoint of bin i. float iMid = PS_BIN_MIDPOINT(histogram, i); @@ -1708,5 +1913,6 @@ // smooth->data.F32[i] = 0.0; - for (long j = jMin ; j <= jMax ; j++) { + for (long j = jMin; j <= jMax; j++) + { float jMid = PS_BIN_MIDPOINT(histogram, j); smooth->data.F32[i] += @@ -1714,9 +1920,12 @@ } } - } else { + } + else + { // // We get here if the histogram is uniform. // - for (long i = 0; i < numBins; i++) { + for (long i = 0; i < numBins; i++) + { psF32 binSize = bounds->data.F32[1] - bounds->data.F32[0]; psS32 gaussWidth = ((PS_GAUSS_WIDTH * sigma) / binSize); @@ -1745,5 +1954,6 @@ smooth->data.F32[i] = 0.0; float iMid = PS_BIN_MIDPOINT(histogram, i); - for (long j = jMin ; j <= jMax ; j++) { + for (long j = jMin; j <= jMax; j++) + { float jMid = PS_BIN_MIDPOINT(histogram, j); smooth->data.F32[i] += @@ -1753,9 +1963,10 @@ } - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(smooth); } psTrace(TRACE, 4, "---- %s() end ----\n", __func__); - return(smooth); + return (smooth); } @@ -1769,7 +1980,8 @@ static void statsFree(psStats *stats) { - if (!stats) return; - psFree (stats->tmpData); - psFree (stats->tmpMask); + if (!stats) + return; + psFree(stats->tmpData); + psFree(stats->tmpMask); return; } @@ -1778,5 +1990,5 @@ psStatsAlloc(): This routine must create a new psStats data structure. *****************************************************************************/ -psStats* p_psStatsAlloc(const char *file, unsigned int lineno, const char *func, psStatsOptions options) +psStats *p_psStatsAlloc(const char *file, unsigned int lineno, const char *func, psStatsOptions options) { psStats *stats = p_psAlloc(file, lineno, func, sizeof(psStats)); @@ -1784,5 +1996,5 @@ // set initial, default values - psStatsInit (stats); + psStatsInit(stats); // these values are can be set as desired by the user. they are not affected by @@ -1818,5 +2030,5 @@ stats->clippedMean = NAN; stats->clippedStdev = NAN; - stats->clippedNvalues = -1; // XXX: This is never used + stats->clippedNvalues = -1; // XXX: This is never used stats->min = NAN; stats->max = NAN; @@ -1825,9 +2037,8 @@ } - bool psMemCheckStats(psPtr ptr) { PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc)statsFree ); + return (psMemGetDeallocator(ptr) == (psFreeFunc)statsFree); } @@ -1847,8 +2058,8 @@ XXX: Should we free stats if the asserts fail? NO; we don't own it (RHL). *****************************************************************************/ -bool psVectorStats(psStats* stats, - const psVector* in, - const psVector* errors, - const psVector* mask, +bool psVectorStats(psStats *stats, + const psVector *in, + const psVector *errors, + const psVector *mask, psVectorMaskType maskVal) { @@ -1856,9 +2067,11 @@ PS_ASSERT_VECTOR_NON_NULL(in, false); PS_ASSERT_VECTOR_NON_EMPTY(in, false); - if (mask) { + if (mask) + { PS_ASSERT_VECTORS_SIZE_EQUAL(mask, in, false); PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_VECTOR_MASK, false); } - if (errors) { + if (errors) + { PS_ASSERT_VECTORS_SIZE_EQUAL(errors, in, false); PS_ASSERT_VECTOR_TYPE(errors, in->type.type, false); @@ -1866,32 +2079,45 @@ // Convert types, as necessary - psVector *inF32 = NULL; // Input vector of values, F32 version - if (in->type.type == PS_TYPE_F32) { + psVector *inF32 = NULL; // Input vector of values, F32 version + if (in->type.type == PS_TYPE_F32) + { inF32 = psMemIncrRefCounter((psPtr)in); - } else { + } + else + { inF32 = psVectorCopy(NULL, in, PS_TYPE_F32); } - psVector *errorsF32 = NULL; // Input vector of errors, F32 version - if (errors) { - if (errors->type.type == PS_TYPE_F32) { + psVector *errorsF32 = NULL; // Input vector of errors, F32 version + if (errors) + { + if (errors->type.type == PS_TYPE_F32) + { errorsF32 = psMemIncrRefCounter((psPtr)errors); - } else { + } + else + { errorsF32 = psVectorCopy(NULL, errors, PS_TYPE_F32); } } - psVector *maskVector = NULL; // Input mask vector, U8 version - if (mask) { - if (mask->type.type == PS_TYPE_VECTOR_MASK) { + psVector *maskVector = NULL; // Input mask vector, U8 version + if (mask) + { + if (mask->type.type == PS_TYPE_VECTOR_MASK) + { maskVector = psMemIncrRefCounter((psPtr)mask); - } else { + } + else + { maskVector = psVectorCopy(NULL, mask, PS_TYPE_VECTOR_MASK); } } - if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) { + if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) + { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->max, stats->min, false); } - if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) { + if ((stats->options & PS_STAT_USE_BINSIZE) && (stats->min >= stats->max)) + { PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(stats->binsize, 0.0, false); } @@ -1903,7 +2129,9 @@ // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_MEAN) { - // NOTE: vectorSampleMean cannot return 'false' - if (!vectorSampleMean(inF32, errorsF32, maskVector, maskVal, stats)) { + if (stats->options & PS_STAT_SAMPLE_MEAN) + { + // NOTE: vectorSampleMean cannot return 'false' + if (!vectorSampleMean(inF32, errorsF32, maskVector, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, "Failed to calculate vector sample mean"); status &= false; @@ -1912,17 +2140,21 @@ // ************************************************************************ - if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) { - // NOTE: vectorSampleMedian only returns 'false' for very bad cases: - // NULL input vector, psVectorCopy failure, invalid vector type (likely programming errors) - if (!vectorSampleMedian(inF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median"); - status &= false; - } + if (stats->options & (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE)) + { + // NOTE: vectorSampleMedian only returns 'false' for very bad cases: + // NULL input vector, psVectorCopy failure, invalid vector type (likely programming errors) + if (!vectorSampleMedian(inF32, maskVector, maskVal, stats)) + { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample median"); + status &= false; + } } // ************************************************************************ - if (stats->options & PS_STAT_SAMPLE_STDEV) { - // NOTE: vectorSampleStdev cannot return 'false' - if (!vectorSampleStdev(inF32, errorsF32, maskVector, maskVal, stats)) { + if (stats->options & PS_STAT_SAMPLE_STDEV) + { + // NOTE: vectorSampleStdev cannot return 'false' + if (!vectorSampleStdev(inF32, errorsF32, maskVector, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample stdev"); status &= false; @@ -1930,7 +2162,9 @@ } - if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS)) { - // NOTE: vectorSampleMoments cannot return 'false' - if (!vectorSampleMoments(inF32, maskVector, maskVal, stats)) { + if (stats->options & (PS_STAT_SAMPLE_SKEWNESS | PS_STAT_SAMPLE_KURTOSIS)) + { + // NOTE: vectorSampleMoments cannot return 'false' + if (!vectorSampleMoments(inF32, maskVector, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, "Failed to calculate sample moments"); status &= false; @@ -1939,14 +2173,17 @@ // ************************************************************************ - if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) { - // NOTE: vectorMinMax returns 0 if there are no valid values, - // but this is not an error condition. stats.min,max are set to NAN. - // vectorMinMax cannot raise an error + if (stats->options & (PS_STAT_MAX | PS_STAT_MIN)) + { + // NOTE: vectorMinMax returns 0 if there are no valid values, + // but this is not an error condition. stats.min,max are set to NAN. + // vectorMinMax cannot raise an error vectorMinMax(inF32, maskVector, maskVal, stats); } // ************************************************************************ - if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) { - if (!vectorRobustStats(inF32, errorsF32, maskVector, maskVal, stats)) { + if (stats->options & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE)) + { + if (!vectorRobustStats(inF32, errorsF32, maskVector, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, _("Failed to calculate robust statistics")); status &= false; @@ -1955,6 +2192,8 @@ // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - if (!vectorFittedStats(inF32, errorsF32, maskVector, maskVal, stats)) { + if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) + { + if (!vectorFittedStats(inF32, errorsF32, maskVector, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); status &= false; @@ -1963,6 +2202,8 @@ // ************************************************************************ - if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { - if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) { + if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_MEDIAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) + { + if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) + { psError(PS_ERR_UNKNOWN, false, "Failed to calculate clipped statistics\n"); status &= false; @@ -1980,40 +2221,42 @@ PS_ASSERT_STRING_NON_EMPTY(string, PS_STAT_NONE); -#define READ_STAT(NAME, SYMBOL) \ - if (strcasecmp(string, NAME) == 0) { \ - return SYMBOL; \ - } - - READ_STAT("MEAN", PS_STAT_SAMPLE_MEAN); - READ_STAT("STDEV", PS_STAT_SAMPLE_STDEV); - READ_STAT("SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); - READ_STAT("KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); - READ_STAT("MEDIAN", PS_STAT_SAMPLE_MEDIAN); - READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE); - READ_STAT("SAMPLE_MEAN", PS_STAT_SAMPLE_MEAN); - READ_STAT("SAMPLE_STDEV", PS_STAT_SAMPLE_STDEV); - READ_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN); +#define READ_STAT(NAME, SYMBOL) \ + if (strcasecmp(string, NAME) == 0) \ + { \ + return SYMBOL; \ + } + + READ_STAT("MEAN", PS_STAT_SAMPLE_MEAN); + READ_STAT("STDEV", PS_STAT_SAMPLE_STDEV); + READ_STAT("SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); + READ_STAT("KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); + READ_STAT("MEDIAN", PS_STAT_SAMPLE_MEDIAN); + READ_STAT("QUARTILE", PS_STAT_SAMPLE_QUARTILE); + READ_STAT("SAMPLE_MEAN", PS_STAT_SAMPLE_MEAN); + READ_STAT("SAMPLE_STDEV", PS_STAT_SAMPLE_STDEV); + READ_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN); READ_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE); READ_STAT("SAMPLE_SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); READ_STAT("SAMPLE_KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); - READ_STAT("ROBUST", PS_STAT_ROBUST_MEDIAN); - READ_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN); - READ_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); + READ_STAT("ROBUST", PS_STAT_ROBUST_MEDIAN); + READ_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN); + READ_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); READ_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE); - READ_STAT("FITTED", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); - READ_STAT("FITTED_V2", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); + READ_STAT("FITTED_V2", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN); READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV); - READ_STAT("FITTED_V3", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_V3", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN); READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV); - READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN); - READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN); + READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN); READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV); - READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); - READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); - READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); + READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); + READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); + READ_STAT("CLIPPED_MEDIAN", PS_STAT_CLIPPED_MEDIAN); + READ_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to parse statistic: %s\n", string); @@ -2023,25 +2266,27 @@ psString psStatsOptionToString(psStatsOptions option) { - psString string = NULL; // String to return - -#define WRITE_STAT(NAME, SYMBOL) \ - if (option & SYMBOL) { \ + psString string = NULL; // String to return + +#define WRITE_STAT(NAME, SYMBOL) \ + if (option & SYMBOL) \ + { \ psStringAppend(&string, "%s ", NAME); \ } // Same list as above (for psStatsFromString), but with repeat symbols removed - WRITE_STAT("SAMPLE_MEAN", PS_STAT_SAMPLE_MEAN); - WRITE_STAT("SAMPLE_STDEV", PS_STAT_SAMPLE_STDEV); - WRITE_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN); + WRITE_STAT("SAMPLE_MEAN", PS_STAT_SAMPLE_MEAN); + WRITE_STAT("SAMPLE_STDEV", PS_STAT_SAMPLE_STDEV); + WRITE_STAT("SAMPLE_MEDIAN", PS_STAT_SAMPLE_MEDIAN); WRITE_STAT("SAMPLE_QUARTILE", PS_STAT_SAMPLE_QUARTILE); WRITE_STAT("SAMPLE_SKEWNESS", PS_STAT_SAMPLE_SKEWNESS); WRITE_STAT("SAMPLE_KURTOSIS", PS_STAT_SAMPLE_KURTOSIS); - WRITE_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN); - WRITE_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); + WRITE_STAT("ROBUST_MEDIAN", PS_STAT_ROBUST_MEDIAN); + WRITE_STAT("ROBUST_STDEV", PS_STAT_ROBUST_STDEV); WRITE_STAT("ROBUST_QUARTILE", PS_STAT_ROBUST_QUARTILE); - WRITE_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); - WRITE_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); - WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); - WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); + WRITE_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); + WRITE_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); + WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); + WRITE_STAT("CLIPPED_MEDIAN", PS_STAT_CLIPPED_MEDIAN); + WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); return string; @@ -2051,5 +2296,6 @@ { psList *subStrings = psStringSplit(string, " ,;", false); // List of sub-strings - if (!subStrings || psListLength(subStrings) == 0) { + if (!subStrings || psListLength(subStrings) == 0) + { // Nothing here psError(PS_ERR_BAD_PARAMETER_VALUE, false, "No string to parse for statistics: %s\n", string); @@ -2057,10 +2303,12 @@ return NULL; } - psStats *stats = psStatsAlloc(0); // Generate empty stats structure + psStats *stats = psStatsAlloc(0); // Generate empty stats structure psListIterator *iterator = psListIteratorAlloc(subStrings, PS_LIST_HEAD, false); // Iterator - psString statString; // Statistic string, from iteration - while ((statString = psListGetAndIncrement(iterator))) { + psString statString; // Statistic string, from iteration + while ((statString = psListGetAndIncrement(iterator))) + { psStatsOptions option = psStatsOptionFromString(statString); - if (option == 0) { + if (option == 0) + { psWarning("Unable to interpret statistic option: %s --- ignored.\n", statString); continue; @@ -2080,22 +2328,24 @@ psStatsOptions psStatsSingleOption(psStatsOptions option) { - switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) { - case PS_STAT_SAMPLE_MEAN: - case PS_STAT_SAMPLE_MEDIAN: - case PS_STAT_SAMPLE_STDEV: - case PS_STAT_SAMPLE_QUARTILE: - case PS_STAT_SAMPLE_SKEWNESS: - case PS_STAT_SAMPLE_KURTOSIS: - case PS_STAT_ROBUST_MEDIAN: - case PS_STAT_ROBUST_STDEV: - case PS_STAT_ROBUST_QUARTILE: - case PS_STAT_FITTED_MEAN: - case PS_STAT_FITTED_STDEV: - case PS_STAT_CLIPPED_MEAN: - case PS_STAT_CLIPPED_STDEV: - case PS_STAT_MAX: - case PS_STAT_MIN: + switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) + { + case PS_STAT_SAMPLE_MEAN: + case PS_STAT_SAMPLE_MEDIAN: + case PS_STAT_SAMPLE_STDEV: + case PS_STAT_SAMPLE_QUARTILE: + case PS_STAT_SAMPLE_SKEWNESS: + case PS_STAT_SAMPLE_KURTOSIS: + case PS_STAT_ROBUST_MEDIAN: + case PS_STAT_ROBUST_STDEV: + case PS_STAT_ROBUST_QUARTILE: + case PS_STAT_FITTED_MEAN: + case PS_STAT_FITTED_STDEV: + case PS_STAT_CLIPPED_MEAN: + case PS_STAT_CLIPPED_MEDIAN: + case PS_STAT_CLIPPED_STDEV: + case PS_STAT_MAX: + case PS_STAT_MIN: return option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE); - default: + default: return 0; } @@ -2107,5 +2357,5 @@ { return options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | - PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN); + PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_MEDIAN | PS_STAT_FITTED_MEAN); } @@ -2116,41 +2366,43 @@ } - double psStatsGetValue(const psStats *stats, psStatsOptions option) { // We could call psStatsSingle to check, but it would be a waste since we effectively do it anyway - switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) { - case PS_STAT_SAMPLE_MEAN: + switch (option & ~(PS_STAT_USE_RANGE | PS_STAT_USE_BINSIZE)) + { + case PS_STAT_SAMPLE_MEAN: return stats->sampleMean; - case PS_STAT_SAMPLE_MEDIAN: + case PS_STAT_SAMPLE_MEDIAN: return stats->sampleMedian; - case PS_STAT_SAMPLE_STDEV: + case PS_STAT_SAMPLE_STDEV: return stats->sampleStdev; - case PS_STAT_SAMPLE_SKEWNESS: + case PS_STAT_SAMPLE_SKEWNESS: return stats->sampleSkewness; - case PS_STAT_SAMPLE_KURTOSIS: + case PS_STAT_SAMPLE_KURTOSIS: return stats->sampleKurtosis; - case PS_STAT_ROBUST_MEDIAN: + case PS_STAT_ROBUST_MEDIAN: return stats->robustMedian; - case PS_STAT_ROBUST_STDEV: + case PS_STAT_ROBUST_STDEV: return stats->robustStdev; - case PS_STAT_FITTED_MEAN: + case PS_STAT_FITTED_MEAN: return stats->fittedMean; - case PS_STAT_FITTED_STDEV: + case PS_STAT_FITTED_STDEV: return stats->fittedStdev; - case PS_STAT_CLIPPED_MEAN: + case PS_STAT_CLIPPED_MEAN: return stats->clippedMean; - case PS_STAT_CLIPPED_STDEV: + case PS_STAT_CLIPPED_MEDIAN: + return stats->clippedMedian; + case PS_STAT_CLIPPED_STDEV: return stats->clippedStdev; - case PS_STAT_MAX: + case PS_STAT_MAX: return stats->max; - case PS_STAT_MIN: + case PS_STAT_MIN: return stats->min; - case PS_STAT_SAMPLE_QUARTILE: - case PS_STAT_ROBUST_QUARTILE: + case PS_STAT_SAMPLE_QUARTILE: + case PS_STAT_ROBUST_QUARTILE: psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Cannot return a single quartile value; " - "get them yourself.\n"); + "get them yourself.\n"); return NAN; - default: + default: return NAN; } @@ -2161,5 +2413,5 @@ // other private functions used above -# if (0) +#if (0) static psF32 QuadraticInverse(psF32 a, psF32 b, @@ -2167,16 +2419,17 @@ psF32 y, psF32 xLo, - psF32 xHi - ) + psF32 xHi) { - psF64 tmp = sqrt((y - c)/a + (b*b)/(4.0 * a * a)); - - psF64 x1 = -b/(2.0*a) + tmp; - psF64 x2 = -b/(2.0*a) - tmp; - - if (xLo <= x1 && x1 <= xHi) { + psF64 tmp = sqrt((y - c) / a + (b * b) / (4.0 * a * a)); + + psF64 x1 = -b / (2.0 * a) + tmp; + psF64 x2 = -b / (2.0 * a) - tmp; + + if (xLo <= x1 && x1 <= xHi) + { return x1; } - if (xLo <= x2 && x2 <= xHi) { + if (xLo <= x2 && x2 <= xHi) + { return x2; } @@ -2185,20 +2438,20 @@ static psF32 LinearInverse(psF32 a, - psF32 b, - psF32 y, - psF32 xLo, - psF32 xHi - ) + psF32 b, + psF32 y, + psF32 xLo, + psF32 xHi) { psF64 x = (y - b) / a; - if (xLo <= x && x <= xHi) { + if (xLo <= x && x <= xHi) + { return x; } return 0.5 * (xLo + xHi); } -# endif - -# if (0) +#endif + +#if (0) /****************************************************************************** fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): A general @@ -2214,10 +2467,10 @@ psVector *yVec, psS32 binNum, - psF32 yVal - ) + psF32 yVal) { psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(xVec); PS_VECTOR_PRINT_F32(yVec); @@ -2236,14 +2489,15 @@ psF32 tmpFloat = 0.0f; - if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->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])); - x->data.F64[1] = (psF64) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); - x->data.F64[2] = (psF64) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + x->data.F64[0] = (psF64)(0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (psF64)(0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1])); + x->data.F64[2] = (psF64)(0.5 * (xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2])); y->data.F64[0] = yVec->data.F32[binNum - 1]; y->data.F64[1] = yVec->data.F32[binNum]; y->data.F64[2] = yVec->data.F32[binNum + 1]; psTrace(TRACE, 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]); + xVec->data.F32[binNum], xVec->data.F32[binNum + 1], xVec->data.F32[binNum + 2]); psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); @@ -2252,6 +2506,7 @@ // Ensure that the y value lies within range of the y values. // - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || - ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[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, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2263,5 +2518,6 @@ // 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"); @@ -2274,5 +2530,6 @@ // Determine the coefficients of the polynomial. psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { + if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to fit a 1-dimensional polynomial to the three specified data points. " @@ -2287,7 +2544,7 @@ psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n", - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2])); + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[2])); psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); @@ -2296,5 +2553,6 @@ psFree(myPoly); - if (isnan(tmpFloat)) { + if (isnan(tmpFloat)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); @@ -2302,18 +2560,25 @@ psFree(y); psTrace(TRACE, 5, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - } else { + return (NAN); + } + } + else + { // These are special cases where the bin is at the beginning or end of the vector. - if (binNum == 0) { + 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]); - } else if (binNum == (xVec->n - 1)) { + } + else if (binNum == (xVec->n - 1)) + { // The special case where we have two points only at the end of // the vectors x and y. // XXX: Is this right? tmpFloat = xVec->data.F32[binNum]; - } else if (binNum == (xVec->n - 2)) { + } + else if (binNum == (xVec->n - 2)) + { // XXX: Is this right? tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); @@ -2341,10 +2606,10 @@ psVector *yVec, psS32 binNum, - psF32 yVal - ) + psF32 yVal) { psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(xVec); PS_VECTOR_PRINT_F32(yVec); @@ -2362,14 +2627,15 @@ psF32 tmpFloat = 0.0f; - if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->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] = xVec->data.F32[binNum - 1]; x->data.F64[1] = xVec->data.F32[binNum]; - x->data.F64[2] = xVec->data.F32[binNum+1]; + x->data.F64[2] = xVec->data.F32[binNum + 1]; y->data.F64[0] = yVec->data.F32[binNum - 1]; y->data.F64[1] = yVec->data.F32[binNum]; y->data.F64[2] = yVec->data.F32[binNum + 1]; psTrace(TRACE, 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]); + xVec->data.F32[binNum], xVec->data.F32[binNum + 1], xVec->data.F32[binNum + 2]); psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); @@ -2378,6 +2644,7 @@ // Ensure that the y value lies within range of the y values. // - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || - ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[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, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2389,5 +2656,6 @@ // 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"); @@ -2400,5 +2668,6 @@ // Determine the coefficients of the polynomial. psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { + if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to fit a 1-dimensional polynomial to the three specified data points. " @@ -2413,7 +2682,7 @@ psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n", - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2])); + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[2])); psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); @@ -2422,5 +2691,6 @@ psFree(myPoly); - if (isnan(tmpFloat)) { + if (isnan(tmpFloat)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); @@ -2428,19 +2698,26 @@ psFree(y); psTrace(TRACE, 5, "---- %s(NAN) end ----\n", __func__); - return(NAN); - } - } else { + return (NAN); + } + } + else + { // These are special cases where the bin is at the beginning or end of the vector. - if (binNum == 0) { + if (binNum == 0) + { // We have two points only at the beginning of the vectors x and y. // XXX this does not seem to be doing the linear interpolation / extrapolation tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); - } else if (binNum == (xVec->n - 1)) { + } + else if (binNum == (xVec->n - 1)) + { // The special case where we have two points only at the end of // the vectors x and y. // XXX: Is this right? tmpFloat = xVec->data.F32[binNum]; - } else if (binNum == (xVec->n - 2)) { + } + else if (binNum == (xVec->n - 2)) + { // XXX: Is this right? tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]); @@ -2472,9 +2749,9 @@ psVector *yVec, psS32 binNum, - psF32 yVal - ) + psF32 yVal) { psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(xVec); PS_VECTOR_PRINT_F32(yVec); @@ -2495,27 +2772,28 @@ // if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) { - if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) { + if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) + { // The general case. We have all three points. - // x->data.F64[0] = binNum - 1; - // x->data.F64[1] = binNum; - // x->data.F64[2] = binNum + 1; - x->data.F64[0] = xVec->data.F32[binNum - 2]; - x->data.F64[1] = xVec->data.F32[binNum - 1]; - x->data.F64[2] = xVec->data.F32[binNum + 0]; - x->data.F64[3] = xVec->data.F32[binNum + 1]; - x->data.F64[4] = xVec->data.F32[binNum + 2]; + // x->data.F64[0] = binNum - 1; + // x->data.F64[1] = binNum; + // x->data.F64[2] = binNum + 1; + x->data.F64[0] = xVec->data.F32[binNum - 2]; + x->data.F64[1] = xVec->data.F32[binNum - 1]; + x->data.F64[2] = xVec->data.F32[binNum + 0]; + x->data.F64[3] = xVec->data.F32[binNum + 1]; + x->data.F64[4] = xVec->data.F32[binNum + 2]; y->data.F64[0] = yVec->data.F32[binNum - 2]; y->data.F64[1] = yVec->data.F32[binNum - 1]; y->data.F64[2] = yVec->data.F32[binNum + 0]; - y->data.F64[3] = yVec->data.F32[binNum + 1]; - y->data.F64[4] = yVec->data.F32[binNum + 2]; - psTrace(TRACE, 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]); + y->data.F64[3] = yVec->data.F32[binNum + 1]; + y->data.F64[4] = yVec->data.F32[binNum + 2]; + psTrace(TRACE, 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(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); - // Ensure that the y value lies within range of the y values. - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || - ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) { + if (!(((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || + ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0])))) + { psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2526,5 +2804,6 @@ // Ensure that the y values are monotonic. 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"); @@ -2536,5 +2815,6 @@ // Determine the coefficients of the polynomial. psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { + if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to fit a 1-dimensional polynomial to the three specified data points. " @@ -2549,7 +2829,7 @@ psTrace(TRACE, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); psTrace(TRACE, 6, "Fitted y vec is (%f %f %f)\n", - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[2])); + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[2])); psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); @@ -2557,43 +2837,54 @@ psFree(myPoly); - if (isnan(binValue)) { + if (isnan(binValue)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); psFree(x); psFree(y); - return(NAN); - } - + return (NAN); + } + // 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; - tmpFloat = binValue; - - } else { + // 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; + tmpFloat = binValue; + } + else + { // These are special cases where the bin is at the beginning or end of the vector. - if (binNum == 0) { + 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]; - if (dY == 0.0) { + if (dY == 0.0) + { tmpFloat = xVec->data.F32[0]; - } else { + } + else + { tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; } - } else if (binNum == (xVec->n - 1)) { + } + 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]; - 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]; + 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]; } } @@ -2606,6 +2897,5 @@ return tmpFloat; } -# endif - +#endif /****************************************************************************** @@ -2623,48 +2913,52 @@ *****************************************************************************/ static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal - ) + psVector *yVec, + psS32 binNum, + psF32 yVal) { -# if (1) -# define HALF_SIZE 2 - double Sx = 0.0; - - double Sy = 0.0; - double Sxx = 0.0; - double Sxy = 0.0; - double deltaY = 0.0; - int N = 0; - - for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) { - if ((u >= 0)&&(u < yVec->n)) { - if (u+1 < xVec->n) { - Sx += yVec->data.F32[u]; - Sxx += PS_SQR(yVec->data.F32[u]); - - deltaY = xVec->data.F32[u]; - //deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]); - Sy += deltaY; - Sxy += yVec->data.F32[u] * deltaY; - N += 1; - } - } - } - double Det = N * Sxx - Sx * Sx; - if (Det == 0.0) return NAN; - if (N == 0) return NAN; - - double C0 = (Sy*Sxx - Sx*Sxy) / Det; - double C1 = (Sxy*N - Sx*Sy) / Det; - - double value = C0 + yVal*C1; - return value; - - -# else +#if (1) +#define HALF_SIZE 2 + double Sx = 0.0; + + double Sy = 0.0; + double Sxx = 0.0; + double Sxy = 0.0; + double deltaY = 0.0; + int N = 0; + + for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) + { + if ((u >= 0) && (u < yVec->n)) + { + if (u + 1 < xVec->n) + { + Sx += yVec->data.F32[u]; + Sxx += PS_SQR(yVec->data.F32[u]); + + deltaY = xVec->data.F32[u]; + // deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]); + Sy += deltaY; + Sxy += yVec->data.F32[u] * deltaY; + N += 1; + } + } + } + double Det = N * Sxx - Sx * Sx; + if (Det == 0.0) + return NAN; + if (N == 0) + return NAN; + + double C0 = (Sy * Sxx - Sx * Sxy) / Det; + double C1 = (Sxy * N - Sx * Sy) / Det; + + double value = C0 + yVal * C1; + return value; + +#else psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); - if (psTraceGetLevel("psLib.math") >= 8) { + if (psTraceGetLevel("psLib.math") >= 8) + { PS_VECTOR_PRINT_F32(xVec); PS_VECTOR_PRINT_F32(yVec); @@ -2684,25 +2978,27 @@ psF32 tmpFloat = 0.0f; - if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) { - x->data.F64[0] = xVec->data.F32[binNum - 2]; - x->data.F64[1] = xVec->data.F32[binNum - 1]; - x->data.F64[2] = xVec->data.F32[binNum + 0]; - x->data.F64[3] = xVec->data.F32[binNum + 1]; - x->data.F64[4] = xVec->data.F32[binNum + 2]; - - y->data.F64[0] = yVec->data.F32[binNum - 2]; - y->data.F64[1] = yVec->data.F32[binNum - 1]; - y->data.F64[2] = yVec->data.F32[binNum + 0]; - y->data.F64[3] = yVec->data.F32[binNum + 1]; - y->data.F64[4] = yVec->data.F32[binNum + 2]; - psTrace(TRACE, 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(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); - psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); - - // Ensure that the y value lies within range of the y values. - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || - ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified yVal, %g, is not within y-range, %g to %g."), + if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) + { + x->data.F64[0] = xVec->data.F32[binNum - 2]; + x->data.F64[1] = xVec->data.F32[binNum - 1]; + x->data.F64[2] = xVec->data.F32[binNum + 0]; + x->data.F64[3] = xVec->data.F32[binNum + 1]; + x->data.F64[4] = xVec->data.F32[binNum + 2]; + + y->data.F64[0] = yVec->data.F32[binNum - 2]; + y->data.F64[1] = yVec->data.F32[binNum - 1]; + y->data.F64[2] = yVec->data.F32[binNum + 0]; + y->data.F64[3] = yVec->data.F32[binNum + 1]; + y->data.F64[4] = yVec->data.F32[binNum + 2]; + psTrace(TRACE, 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(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); + psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); + + // Ensure that the y value lies within range of the y values. + if (!(((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || + ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0])))) + { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + _("Specified yVal, %g, is not within y-range, %g to %g."), (psF64)yVal, y->data.F64[0], y->data.F64[2]); return NAN; @@ -2711,5 +3007,6 @@ // Ensure that the y values are monotonic. 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"); @@ -2721,5 +3018,6 @@ // Determine the coefficients of the polynomial. psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); - if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { + if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to fit a 1-dimensional polynomial to the three specified data points. " @@ -2733,6 +3031,6 @@ psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]); psTrace(TRACE, 6, "Fitted y vec is (%f %f)\n", - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1])); + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[0]), + (psF32)psPolynomial1DEval(myPoly, (psF64)x->data.F64[1])); psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); @@ -2740,44 +3038,54 @@ psFree(myPoly); - if (isnan(binValue)) { + if (isnan(binValue)) + { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); psFree(x); psFree(y); - return(NAN); - } - + return (NAN); + } + // 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; - tmpFloat = binValue; - - - } else { + // 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; + tmpFloat = binValue; + } + else + { // These are special cases where the bin is at the beginning or end of the vector. - if (binNum == 0) { + 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]; - if (dY == 0.0) { + if (dY == 0.0) + { tmpFloat = xVec->data.F32[0]; - } else { + } + else + { tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; } - } else if (binNum == (xVec->n - 1)) { + } + 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]; - 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]; + 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]; } } @@ -2789,4 +3097,4 @@ return tmpFloat; -# endif +#endif } Index: /branches/2dbias/psLib/src/math/psStats.h =================================================================== --- /branches/2dbias/psLib/src/math/psStats.h (revision 42718) +++ /branches/2dbias/psLib/src/math/psStats.h (revision 42719) @@ -27,24 +27,26 @@ * @see psStats, psVectorStats, psImageStats */ -typedef enum { - PS_STAT_NONE = 0x000000, ///< Empty set - PS_STAT_MIN = 0x000001, ///< Maximum - PS_STAT_MAX = 0x000002, ///< Minumum - PS_STAT_SAMPLE_MEAN = 0x000004, ///< Sample Mean - PS_STAT_SAMPLE_MEDIAN = 0x000008, ///< Sample Median - PS_STAT_SAMPLE_STDEV = 0x000010, ///< Sample Standard Deviation - PS_STAT_SAMPLE_QUARTILE = 0x000020, ///< Sample Quartile - PS_STAT_SAMPLE_SKEWNESS = 0x000040, ///< Sample Skewness (third moment) - PS_STAT_SAMPLE_KURTOSIS = 0x000080, ///< Sample Kurtosis (fourth moment) - PS_STAT_ROBUST_MEDIAN = 0x000100, ///< Robust Median - PS_STAT_ROBUST_STDEV = 0x000200, ///< Robust Standarad Deviation - PS_STAT_ROBUST_QUARTILE = 0x000400, ///< Robust Quartile - PS_STAT_ROBUST_SPARE1 = 0x000800, ///< Spare 1 - PS_STAT_FITTED_MEAN = 0x001000, ///< Fitted Mean - PS_STAT_FITTED_STDEV = 0x002000, ///< Fitted Standard Deviation - PS_STAT_CLIPPED_MEAN = 0x040000, ///< Clipped Mean - PS_STAT_CLIPPED_STDEV = 0x080000, ///< Clipped Standard Deviation - PS_STAT_USE_RANGE = 0x100000, ///< Range - PS_STAT_USE_BINSIZE = 0x200000, ///< Binsize +typedef enum +{ + PS_STAT_NONE = 0x000000, ///< Empty set + PS_STAT_MIN = 0x000001, ///< Maximum + PS_STAT_MAX = 0x000002, ///< Minumum + PS_STAT_SAMPLE_MEAN = 0x000004, ///< Sample Mean + PS_STAT_SAMPLE_MEDIAN = 0x000008, ///< Sample Median + PS_STAT_SAMPLE_STDEV = 0x000010, ///< Sample Standard Deviation + PS_STAT_SAMPLE_QUARTILE = 0x000020, ///< Sample Quartile + PS_STAT_SAMPLE_SKEWNESS = 0x000040, ///< Sample Skewness (third moment) + PS_STAT_SAMPLE_KURTOSIS = 0x000080, ///< Sample Kurtosis (fourth moment) + PS_STAT_ROBUST_MEDIAN = 0x000100, ///< Robust Median + PS_STAT_ROBUST_STDEV = 0x000200, ///< Robust Standarad Deviation + PS_STAT_ROBUST_QUARTILE = 0x000400, ///< Robust Quartile + PS_STAT_ROBUST_SPARE1 = 0x000800, ///< Spare 1 + PS_STAT_FITTED_MEAN = 0x001000, ///< Fitted Mean + PS_STAT_FITTED_STDEV = 0x002000, ///< Fitted Standard Deviation + PS_STAT_CLIPPED_MEDIAN = 0x020000, ///< Clipped Median + PS_STAT_CLIPPED_MEAN = 0x040000, ///< Clipped Mean + PS_STAT_CLIPPED_STDEV = 0x080000, ///< Clipped Standard Deviation + PS_STAT_USE_RANGE = 0x100000, ///< Range + PS_STAT_USE_BINSIZE = 0x200000, ///< Binsize } psStatsOptions; @@ -54,34 +56,34 @@ typedef struct { - double sampleMean; ///< formal mean of sample - double sampleMedian; ///< formal median of sample - double sampleStdev; ///< standard deviation of sample - double sampleUQ; ///< upper quartile of sample - double sampleLQ; ///< lower quartile of sample - double sampleSkewness; ///< skewness (third moment) of sample - double sampleKurtosis; ///< kurtosis (fourth moment) of sample - double robustMedian; ///< robust median of array - double robustStdev; ///< robust standard deviation of array - double robustUQ; ///< robust upper quartile - double robustLQ; ///< robust lower quartile - long robustN50; ///< Number of points in Gaussian fit; XXX: This is currently unused. - double fittedMean; ///< robust mean of data - double fittedStdev; ///< robust standard deviation of data - long fittedNfit; ///< Number of points in Gaussian fit; XXX: This is currently unused - double clippedMean; ///< Nsigma clipped mean - double clippedStdev; ///< standard deviation after clipping - long clippedNvalues; ///< Number of data points used for clipped mean. - double clipSigma; ///< Nsigma used for clipping; user input - int clipIter; ///< Number of clipping iterations; user input - double min; ///< minimum data value in array - double max; ///< maximum data value in array - double binsize; ///< binsize for robust fit (input/ouput) - long nSubsample; ///< maxinum number of measurements (input) - psStatsOptions options; ///< bitmask of values requested - psStatsOptions results; ///< bitmask of values calculated - psVector *tmpData; ///< temporary vector so repeated calls do not have to realloc - psVector *tmpMask; ///< temporary vector so repeated calls do not have to realloc -} -psStats; + double sampleMean; ///< formal mean of sample + double sampleMedian; ///< formal median of sample + double sampleStdev; ///< standard deviation of sample + double sampleUQ; ///< upper quartile of sample + double sampleLQ; ///< lower quartile of sample + double sampleSkewness; ///< skewness (third moment) of sample + double sampleKurtosis; ///< kurtosis (fourth moment) of sample + double robustMedian; ///< robust median of array + double robustStdev; ///< robust standard deviation of array + double robustUQ; ///< robust upper quartile + double robustLQ; ///< robust lower quartile + long robustN50; ///< Number of points in Gaussian fit; XXX: This is currently unused. + double fittedMean; ///< robust mean of data + double fittedStdev; ///< robust standard deviation of data + long fittedNfit; ///< Number of points in Gaussian fit; XXX: This is currently unused + double clippedMean; ///< Nsigma clipped mean + double clippedMedian; ///< Nsigma clipped median + double clippedStdev; ///< standard deviation after clipping + long clippedNvalues; ///< Number of data points used for clipped mean. + double clipSigma; ///< Nsigma used for clipping; user input + int clipIter; ///< Number of clipping iterations; user input + double min; ///< minimum data value in array + double max; ///< maximum data value in array + double binsize; ///< binsize for robust fit (input/ouput) + long nSubsample; ///< maxinum number of measurements (input) + psStatsOptions options; ///< bitmask of values requested + psStatsOptions results; ///< bitmask of values calculated + psVector *tmpData; ///< temporary vector so repeated calls do not have to realloc + psVector *tmpMask; ///< temporary vector so repeated calls do not have to realloc +} psStats; /** Performs statistical calculations on a vector. @@ -90,8 +92,8 @@ */ bool psVectorStats( - psStats* stats, ///< stats structure defines stats to be calculated and how - const psVector* in, ///< Vector to be analysed. - const psVector* errors, ///< Errors. - const psVector* mask, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL + psStats *stats, ///< stats structure defines stats to be calculated and how + const psVector *in, ///< Vector to be analysed. + const psVector *errors, ///< Errors. + const psVector *mask, ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL psVectorMaskType maskVal ///< Only mask elements with one of these bits set in maskVector ); @@ -103,16 +105,16 @@ */ #ifdef DOXYGEN -psStats* psStatsAlloc( - psStatsOptions options ///< Statistics to calculate +psStats *psStatsAlloc( + psStatsOptions options ///< Statistics to calculate ); #else // ifdef DOXYGEN -psStats* p_psStatsAlloc( - const char *file, ///< File of caller - unsigned int lineno, ///< Line number of caller - const char *func, ///< Function name of caller - psStatsOptions options ///< Statistics to calculate -) PS_ATTR_MALLOC; +psStats *p_psStatsAlloc( + const char *file, ///< File of caller + unsigned int lineno, ///< Line number of caller + const char *func, ///< Function name of caller + psStatsOptions options ///< Statistics to calculate + ) PS_ATTR_MALLOC; #define psStatsAlloc(options) \ - p_psStatsAlloc(__FILE__, __LINE__, __func__, options) + p_psStatsAlloc(__FILE__, __LINE__, __func__, options) #endif // ifdef DOXYGEN @@ -124,5 +126,5 @@ */ bool psMemCheckStats( - psPtr ptr ///< the pointer whose type to check + psPtr ptr ///< the pointer whose type to check ); Index: /branches/2dbias/psLib/src/math/psVectorSmooth.c =================================================================== --- /branches/2dbias/psLib/src/math/psVectorSmooth.c (revision 42718) +++ /branches/2dbias/psLib/src/math/psVectorSmooth.c (revision 42719) @@ -172,6 +172,7 @@ if (robust) \ { \ - statistic = PS_STAT_CLIPPED_MEAN; \ + statistic = PS_STAT_CLIPPED_MEDIAN; \ stats = psStatsAlloc(statistic); \ + stats->clipSigma = 2.5; \ } \ else \ Index: /branches/2dbias/psModules/src/detrend/pmOverscan.c =================================================================== --- /branches/2dbias/psModules/src/detrend/pmOverscan.c (revision 42718) +++ /branches/2dbias/psModules/src/detrend/pmOverscan.c (revision 42719) @@ -415,5 +415,5 @@ // Reduce the overscans // XXX need to save 2 different chi-square values - psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels, false); + psVector *yReduced = pmOverscanVector(&chi2, overscanOpts->primary, yscanPixels, true); psFree(yscanPixels); if (!yReduced) @@ -453,9 +453,9 @@ } // Now normVector holds the sliced data; then, compute its robust mean - psStatsOptions statistic = PS_STAT_CLIPPED_MEAN; + psStatsOptions statistic = PS_STAT_ROBUST_MEDIAN; psStats *stats = psStatsAlloc(statistic); if (!psVectorStats(stats, normVector, NULL, NULL, 0)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure clipped mean as normalization of xReduced"); + psError(PS_ERR_UNKNOWN, false, "failure to measure robust median as the normalization of xReduced"); psFree(stats); psFree(normVector);