Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 5090) +++ trunk/psLib/src/math/psStats.c (revision 5113) @@ -7,19 +7,19 @@ * on those data structures. * - * @author GLG, MHPCC - * - * XXX: The following stats members are never used, or set in this code. - * stats->robustN50 - * stats->clippedNvalues - * stats->binsize - * - * - * - * - * @version $Revision: 1.146 $ $Name: not supported by cvs2svn $ - * @date $Date: 2005-09-22 02:47:16 $ - * - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - */ +* @author GLG, MHPCC +* +* XXX: The following stats members are never used, or set in this code. +* stats->robustN50 +* stats->clippedNvalues +* stats->binsize +* +* +* +* +* @version $Revision: 1.147 $ $Name: not supported by cvs2svn $ +* @date $Date: 2005-09-23 23:01:30 $ +* +* Copyright 2004 Maui High Performance Computing Center, University of Hawaii +*/ #include @@ -631,4 +631,6 @@ psF32 sigma) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psTrace(__func__, 6, "(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); @@ -666,4 +668,5 @@ false, PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); } @@ -681,4 +684,5 @@ false, PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); + psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); return(NULL); } @@ -703,9 +707,15 @@ // We get here if the histogram is uniform. // - for (psS32 i = 0; i < numBins; i++) { psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); + // + // XXX: The following is wrong. However, in practice, the sigma was too + // large compared to the binSize. This meant that the smoothing was done + // over 500 bins in the robust stats algorithm. This mean that the smoothing + // took way too long. + // + gaussWidth = 10; // // We determine the bin numbers (jMin:jMax) corresponding to a @@ -735,129 +745,7 @@ } + psTrace(__func__, 3, "---- %s(psVector) end ----\n", __func__); return(smooth); } -/****************************************************************************** -p_psVectorSmoothHistGaussianNEW(): This routine smoothes the data in the input -robustHistogram with a Gaussian of width sigma. It returns a psVector of the -smoothed data. - -XXX: Only PS_TYPE_F32 is supported. - -XXX: Write a general routine which smoothes a psVector. This routine should -call that. Is that possible? - -XXX: This is, or will be, prettier than the previous -psVectorSmoothHistGaussian(). However, it is not being used, yet. - *****************************************************************************/ -psVector *p_psVectorSmoothHistGaussianNEW(psHistogram *histogram, - psF32 sigma) -{ - PS_ASSERT_PTR_NON_NULL(histogram, NULL); - PS_ASSERT_PTR_NON_NULL(histogram->bounds, NULL); - PS_ASSERT_PTR_NON_NULL(histogram->nums, NULL); - - psS32 numBins = histogram->nums->n; - psS32 numBounds = histogram->bounds->n; - psVector *smooth = psVectorAlloc(numBins, PS_TYPE_F32); - psF32 firstBound = histogram->bounds->data.F32[0]; - psF32 lastBound = histogram->bounds->data.F32[numBounds-1]; - psScalar x; - x.type.type = PS_TYPE_F32; - psS32 jMin = 0; - psS32 jMax = 0; - - if (histogram->uniform == false) { - // - // We get here if the histogram is non-uniform. - // - - for (psS32 i = 0; i < numBins; i++) { - // Determine the midpoint of bin i. - psS32 iMid = PS_BIN_MIDPOINT(histogram, i); - - // - // We determine the bin numbers (jMin:jMax) corresponding to a - // range of data values surrounding iMid. The range is of size: - // 2*PS_GAUSS_WIDTH*sigma - // - x.data.F32 = iMid - (PS_GAUSS_WIDTH * sigma); - if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMin = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); - if (jMin < 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); - return(NULL); - } - } else if (x.data.F32 <= firstBound) { - jMin = 0; - } else if (x.data.F32 >= lastBound) { - jMin = histogram->bounds->n - 1; - } - - x.data.F32 = iMid + (PS_GAUSS_WIDTH * sigma); - if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { - jMax = p_psVectorBinDisect( *(psVector* *)&histogram->bounds, &x); - if (jMax < 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_STATS_VECTOR_BIN_DISECT_PROBLEM); - return(NULL); - } - } else if (x.data.F32 <= firstBound) { - jMax = 0; - } else if (x.data.F32 >= lastBound) { - jMax = histogram->bounds->n - 1; - } - - // - // Loop from jMin to jMax, computing the gaussian of data i. - // - smooth->data.F32[i] = 0.0; - for (psS32 j = jMin ; j <= jMax ; j++) { - psS32 jMid = PS_BIN_MIDPOINT(histogram, j); - smooth->data.F32[i] += - histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); - } - } - } else { - // - // We get here if the histogram is uniform. - // - - for (psS32 i = 0; i < numBins; i++) { - psF32 binSize = histogram->bounds->data.F32[1] - histogram->bounds->data.F32[0]; - psS32 gaussWidth = (psS32) ((PS_GAUSS_WIDTH * sigma) / binSize); - - // - // We determine the bin numbers (jMin:jMax) corresponding to a - // range of data values surrounding iMid. The range is of size: - // 2*PS_GAUSS_WIDTH*sigma - // - psS32 jMin = i - gaussWidth; - if (jMin < 0 ) { - jMin = 0; - } - psS32 jMax = i + gaussWidth; - if (jMax > (histogram->bounds->n - 1)) { - jMax = (histogram->bounds->n - 1); - } - - // - // Loop from jMin to jMax, computing the gaussian of data i. - // - smooth->data.F32[i] = 0.0; - psS32 iMid = PS_BIN_MIDPOINT(histogram, i); - for (psS32 j = jMin ; j <= jMax ; j++) { - psS32 jMid = PS_BIN_MIDPOINT(histogram, j); - smooth->data.F32[i] += - histogram->nums->data.F32[j] * psGaussian(jMid, iMid, sigma, true); - } - } - } - - return(smooth); -} - /****************************************************************************** p_psVectorSampleQuartiles(myVector, maskVector, maskVal, stats): calculates @@ -1435,7 +1323,5 @@ psF32 fHi = psPolynomial1DEval(myPoly, rangeHigh); if (!((fLo <= getThisValue) && (fHi >= getThisValue))) { - psError(PS_ERR_UNKNOWN, - true, - PS_ERRORTEXT_psStats_STATS_POLY_MEDIAN_OUT_OF_RANGE); + psError(PS_ERR_UNKNOWN, true, "The requested y value (%f) does not fall within the specified range (%f to %f)\n", getThisValue, fLo, fHi); return(NAN); } @@ -1491,4 +1377,11 @@ psF32 yVal) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + psTrace(__func__, 6, "binNum, yVal is (%d, %f)\n", binNum, yVal); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(xVec); + PS_VECTOR_PRINT_F32(yVec); + } + PS_ASSERT_VECTOR_NON_NULL(xVec, NAN); PS_ASSERT_VECTOR_NON_NULL(yVec, NAN); @@ -1502,5 +1395,4 @@ psVector *y = psVectorAlloc(3, PS_TYPE_F64); psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); - // XXX: Why was this 2 when the alloc function specified number of terms? Note: it's correct now. psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); @@ -1515,4 +1407,7 @@ y->data.F64[1] = yVec->data.F32[binNum]; y->data.F64[2] = yVec->data.F32[binNum + 1]; + psTrace(__func__, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]); + psTrace(__func__, 6, "x vec is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); + psTrace(__func__, 6, "y vec is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); // @@ -1529,7 +1424,14 @@ psFree(y); psFree(yErr); + psTrace(__func__, 3, "---- %s(NAN) end ----\n", __func__); return(NAN); } - } else { + // Ensure that yVal is within the range of the bins we are using. + if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, + (psF64)yVal, y->data.F64[2], y->data.F64[0]); + } + } else if (y->data.F64[0] > y->data.F64[1]) { if (!(y->data.F64[1] >= y->data.F64[2])) { psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n"); @@ -1538,14 +1440,15 @@ psFree(y); psFree(yErr); + psTrace(__func__, 3, "---- %s(NAN) end ----\n", __func__); return(NAN); } - } - - // Ensure that yVal is within the range of the bins we are using. - if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, - (psF64)yVal,y->data.F64[2],y->data.F64[0]); - } + // Ensure that yVal is within the range of the bins we are using. + if (!((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE, + (psF64)yVal, y->data.F64[2], y->data.F64[0]); + } + } + yErr->data.F64[0] = 1.0; yErr->data.F64[1] = 1.0; @@ -1553,4 +1456,5 @@ // Determine the coefficients of the polynomial. + // myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x); if (myPoly == NULL) { @@ -1562,9 +1466,19 @@ psFree(y); psFree(yErr); + psTrace(__func__, 3, "---- %s(NAN) end ----\n", __func__); return(NAN); } + psTrace(__func__, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]); + psTrace(__func__, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]); + psTrace(__func__, 6, "myPoly->coeff[2] is %f\n", myPoly->coeff[2]); + psTrace(__func__, 6, "Fitted y vec is (%.2f %.2f %.2f)\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])); + // Call p_ps1DPolyMedian(), which does a binary search on the // polynomial, looking for the value x such that f(x) = yVal + psTrace(__func__, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); + if (isnan(tmpFloat)) { psError(PS_ERR_UNEXPECTED_NULL, @@ -1575,4 +1489,5 @@ psFree(y); psFree(yErr); + psTrace(__func__, 3, "---- %s(NAN) end ----\n", __func__); return(NAN); } @@ -1596,10 +1511,126 @@ } + psTrace(__func__, 6, "FIT: return %f\n", tmpFloat); psFree(myPoly); psFree(x); psFree(y); psFree(yErr); + + psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, tmpFloat); return(tmpFloat); } + +/***************************************************************************** +XXX: Is there a psLib function for this? + *****************************************************************************/ +psVector *PsVectorDup(psVector *in) +{ + PS_ASSERT_VECTOR_NON_NULL(in, NULL); + psVector *out = psVectorAlloc(in->n, in->type.type); + + if (in->type.type == PS_TYPE_F32) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F32[i] = in->data.F32[i]; + } + } else if (in->type.type == PS_TYPE_F64) { + for (psS32 i = 0 ; i < in->n ; i++) { + out->data.F64[i] = in->data.F64[i]; + } + } else { + psError(PS_ERR_UNKNOWN, true, "Unallowable vector type.\n"); + return(NULL); + } + return(out); +} + +/****************************************************************************** +XXX: This function need to be written. Actually, it simply needs to be +retrieved from the CVS repository, since it was written earlier, then +discarded. At present, it was deleted from the CVS repository, so we might +have to retrieve it from tape. +*****************************************************************************/ +psVector *Fit1DGaussian(psVector *x, psVector*y) +{ + psError(PS_ERR_UNKNOWN, true, "This code has not been implemented yet.\n"); + // XXX: This function was previously part of psStats.c, was removed, was + // purged from CVS, and now needs to be retrieved from tape. + + return(NULL); +} + + + +/****************************************************************************** +XXX: We assume unnormalized gaussians. + *****************************************************************************/ +psF32 psMinimizeLMChi2Gauss1D(psVector *deriv, + const psVector *params, + const psVector *coords) +{ + PS_ASSERT_VECTOR_NON_NULL(params, NAN); + PS_ASSERT_VECTOR_SIZE(params, 2, NAN); + PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN); + PS_ASSERT_VECTOR_NON_NULL(coords, NAN); + PS_ASSERT_VECTOR_SIZE(coords, 1, NAN); + PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN); + + psF32 x = coords->data.F32[0]; + psF32 mean = params->data.F32[0]; + psF32 stdev = params->data.F32[1]; + + if (deriv == NULL) { + deriv = psVectorAlloc(2, PS_TYPE_F32); + } else { + PS_ASSERT_VECTOR_SIZE(deriv, 2, NAN); + PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN); + } + + psF32 tmp = (x - mean) * psGaussian(x, mean, stdev, false); + deriv->data.F32[0] = tmp / (stdev * stdev); + tmp = (x - mean) * (x - mean) * psGaussian(x, mean, stdev, 0); + deriv->data.F32[1] = tmp / (stdev * stdev * stdev); + + // printf("f(x, mean, stdev) = f(%.2f %.2f %.2f) is %.2f\n", x, mean, stdev, psGaussian(x, mean, stdev, false)); + + return(psGaussian(x, mean, stdev, false)); +} + +psF32 LinInterpolate(psF32 x0, + psF32 x1, + psF32 y0, + psF32 y1, + psF32 y) +{ + // printf("HEY: LinInterpolate(%.2f %.2f %.2f %.2f %.2f)\n", x0, x1, y0, y1, y); + if (y0 == y1) { + psLogMsg(__func__, PS_LOG_WARN, "WARNING: y0 == y1. Cannot interpolate. Returning NaN.\n"); + return(NAN); + } + // printf("HMMM: LinInterpolate(%.2f %.2f %.2f %.2f %.2f %.2f)\n", x0, x1, y0, y1, y, 0 + y * (x1 - x0) / (y1 - y0)); + + return(x0 + y * (x1 - x0) / (y1 - y0)); +} + +void PRINT_LEAKS() +{ + psScalar *tmpScalar = psScalarAlloc(0.0, PS_TYPE_F32); + psMemCheckCorruption( 1 ); + psS32 currentId = psMemGetId(); + // psMemBlock** blks; + // psS32 nLeaks = psMemCheckLeaks( currentId, &blks, stderr, false ); + psS32 nLeaks = psMemCheckLeaks(currentId,NULL,stderr,false); + printf("COOL: %d leaks.\n", nLeaks); + tmpScalar->data.F32 = 0.0; +} + +/* +#define PRINT_MEMLEAKS(NUM) \ + printf("A: ---------------------------------------- (ID: %d)\n", NUM); \ + memLeaks = psMemCheckLeaks(currentId,NULL,stderr,false); \ + printf("B: ---------------------------------------- (%d)\n", memLeaks); \ +*/ +#define PRINT_MEMLEAKS(NUM) \ +memLeaks = currentId; + /****************************************************************************** @@ -1625,339 +1656,5 @@ XXX: Check for errors in psLib routines that we call. *****************************************************************************/ -psS32 p_psVectorRobustStats(const psVector* myVector, - const psVector* errors, - const psVector* maskVector, - psU32 maskVal, - psStats* stats) -{ - psHistogram* robustHistogram = NULL; - psVector* robustHistogramVector = NULL; - psF32 binSize = 0.0; // Size of the histogram bins - psS32 LQBinNum = -1; // Bin num for lower quartile - psS32 UQBinNum = -1; // Bin num for upper quartile - psS32 medianBinNum = -1; - psS32 i = 0; // Loop index variable - psS32 modeBinNum = 0; - psF32 modeBinCount = 0.0; - psF32 dL = 0.0; - psS32 numBins = 0; - psF32 myMean = 0.0; - psF32 myStdev = 0.0; - psF32 countFloat = 0.0; - psF32 diff = 0.0; - psF32 sumSquares = 0.0; - psF32 sumDiffs = 0.0; - psVector* cumulativeRobustSums = NULL; - psF32 sumRobust = 0.0; - psF32 sumN50 = 0.0; - psF32 sumNfit = 0.0; - psScalar tmpScalar; - tmpScalar.type.type = PS_TYPE_F32; - psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); - - // Compute the initial bin size of the robust histogram. This is done - // by computing the clipped standard deviation of the vector, and dividing - // that by 10.0; - //XXX: add errors - psS32 rc = p_psVectorClippedStats(myVector, NULL, maskVector, maskVal, tmpStats); - if (rc != 0) { - psError(PS_ERR_UNEXPECTED_NULL, - false, - PS_ERRORTEXT_psStats_ROBUST_STATS_CLIPPED_STATS); - return(1); - } - binSize = tmpStats->clippedStdev / 10.0f; - - // If stats->clippedStdev == 0.0, then all data elements have the same - // value. Therefore, we can set the appropiate results and return. - if (fabs(binSize) <= FLT_EPSILON) { - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = tmpStats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - stats->robustMedian = tmpStats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = tmpStats->clippedMean; - } - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = 0.0; - } - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - stats->robustUQ = tmpStats->clippedMean; - stats->robustLQ = tmpStats->clippedMean; - } - // XXX: Set these to the number of unmasked data points? - stats->robustNfit = 0.0; - stats->robustN50 = 0.0; - psFree(tmpStats); - return(0); - } - - // Determine minimum and maximum values in the data vector. - // XXX: remove this conditional? - if (isnan(tmpStats->min)) { - if (0 != p_psVectorMin(myVector, maskVector, maskVal, tmpStats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: p_psVectorMin(): p_psVectorMin() reported a NAN mean.\n"); - return(1); - } - } - // XXX: remove this conditional? - if (isnan(tmpStats->max)) { - if (0 != p_psVectorMax(myVector, maskVector, maskVal, tmpStats)) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: p_psVectorMin(): p_psVectorMax() reported a NAN mean.\n"); - return(1); - } - } - - // Create the histogram structure. NOTE: we can not specify the bin size - // precisely since the argument to psHistogramAlloc() is the number of - // bins, not the binSize. Also, if we get here, we know that - // binSize != 0.0. - numBins = (psS32)((tmpStats->max - tmpStats->min) / binSize); - robustHistogram = psHistogramAlloc(tmpStats->min, tmpStats->max, numBins); - - // Populate the histogram array. - psVectorHistogram(robustHistogram, myVector, errors, maskVector, maskVal); - - // Smooth the histogram, Gaussian-style. - robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram, - tmpStats->clippedStdev / 4.0f); - - /************************************************************************** - Determine the median/lower/upper quartile bin numbers. - - We define a vector called "cumulativeRobustSums" where the value at - index position i is equal to the sum of bins 0:i. This will be used in - determining the median and lower/upper quartiles. - **************************************************************************/ - cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0]; - for (i = 1; i < robustHistogramVector->n; i++) { - cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] + - robustHistogramVector->data.F32[i]; - } - sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; - - tmpScalar.data.F32 = sumRobust / 4.0; - LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); - tmpScalar.data.F32 = 3.0 * sumRobust / 4.0; - UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); - tmpScalar.data.F32 = sumRobust / 2.0; - medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); - - if ((LQBinNum < 0) || (UQBinNum < 0)) { - psError(PS_ERR_UNKNOWN, true, - PS_ERRORTEXT_psStats_ROBUST_QUARTILE_BINS_FAILED); - return(1); - } - if (medianBinNum < 0) { - psError(PS_ERR_UNKNOWN, true, - PS_ERRORTEXT_psStats_ROBUST_QUARTILE_BINS_FAILED); - return(1); - } - /************************************************************************** - Determine the mode in the range LQ:UQ. - **************************************************************************/ - // Determine the bin with the peak value in the range LQ to UQ. - modeBinNum = LQBinNum; - modeBinCount = robustHistogramVector->data.F32[LQBinNum]; - sumN50 = robustHistogram->nums->data.F32[LQBinNum]; - for (i = LQBinNum + 1; i <= UQBinNum; i++) { - if (robustHistogramVector->data.F32[i] > modeBinCount) { - modeBinNum = i; - modeBinCount = robustHistogramVector->data.F32[i]; - } - sumN50 += robustHistogram->nums->data.F32[i]; - } - - dL = (UQBinNum - LQBinNum) / 4; - /************************************************************************** - Determine the mean/stdev for the bins in the range mode-dL to mode+dL - **************************************************************************/ - // Calculate the mean of the smoothed robust histogram in the range - // mode-dL to mode+dL. We use the midpoint of each bin as the mean for - // that bin (this is a non-exact approximation). - sumNfit = 0.0; - myMean = 0.0; - for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i); - countFloat += robustHistogramVector->data.F32[i]; - } - - sumNfit += robustHistogram->nums->data.F32[i]; - } - // XXX: divide by zero? - myMean /= countFloat; - - // Calculate the stdev of the smoothed robust histogram in the range - // mode-dL to mode+dL. We use the midpoint of each bin as the mean for - // that bin. - for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean; - sumSquares += diff * diff * robustHistogramVector->data.F32[i]; - sumDiffs += diff * robustHistogramVector->data.F32[i]; - } - } - myStdev = sqrtf((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - - p_psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0); - - if ((stats->options & PS_STAT_ROBUST_MEAN) || - (stats->options & PS_STAT_ROBUST_STDEV)) { - // We fit a 1-D polynomial to the data. - // XXX: I implement the changes requested in bug 366. However, this - // code no longer produces sensible results. - // XXX: Since we are no longer fitting a 1-D Gaussian, we can probably - // remove some of the above code that calculated the initial estimate - // for the mean and sigma. - - psVector *y = psVectorAlloc(2 * dL + 1, PS_TYPE_F32); - psVector *x = psVectorAlloc(2 * dL + 1, PS_TYPE_F32); - for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { - int index = i - modeBinNum + dL; - // XXX: Should this be the natural log? - // y->data.F32[index] = robustHistogramVector->data.F32[i]; - y->data.F32[index] = logf(robustHistogramVector->data.F32[i]); - x->data.F32[index] = (psF32) index; - } - - psPolynomial1D *tmpPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - // XXX: What about the NULL x argument? - tmpPoly = psVectorFitPolynomial1D(tmpPoly, NULL, 0, y, NULL, NULL); - if (tmpPoly == NULL) { - psLogMsg(__func__, PS_LOG_WARN, - "WARNING: failed fit a 1D polynomial.\n"); - } - psF32 polyFitSigma = sqrtf(-0.5 / tmpPoly->coeff[2]); - psF32 polyFitMean = tmpPoly->coeff[1] * PS_SQR(polyFitSigma); - // psF32 polyFitNorm = exp(tmpPoly->coedd[0] + PS_SQR(polyFitMean) / (2.0 * PS_SQR(polyFitSigma))); - - if (stats->options & PS_STAT_ROBUST_MEAN) { - stats->robustMean = polyFitMean; - } - - if (stats->options & PS_STAT_ROBUST_STDEV) { - stats->robustStdev = polyFitSigma; - } - // For testing only. This shows that the polynomial fit is successful. - if (0) { - printf("--- Results of the 1-D polynomial fit:\n"); - for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { - int index = i - modeBinNum + dL; - printf("(x, y, poly(x)) is (%.2f, %.2f, %.2f)\n", - x->data.F32[index], - y->data.F32[index], - psPolynomial1DEval(tmpPoly, x->data.F32[index])); - } - } - - psFree(x); - psFree(y); - psFree(tmpPoly); - } - - - /************************************************************************** - Set the appropriate members in the output stats struct. - **************************************************************************/ - - if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum); - } - - - // To determine the median, we fit a quadratic y=f(x) to the three bins - // surrounding the bin containing the median (x is the midpoint of each - // bin and y is the value of each bin). Then we figure out what value - // of x corresponds to f(x) being the median (half of all points). - if (stats->options & PS_STAT_ROBUST_MEDIAN) { - // Take a psVector. Fit a polynomial y = f(x) to the 3 data elements - // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0. - - //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. - //Determine id that is okay. - stats->robustMedian = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - cumulativeRobustSums, - medianBinNum, - sumRobust/2.0); - } - - // To determine the quartiles, we fit a quadratic y=f(x) to the three bins - // surrounding the bin containing LQ/UQ (x is the midpoint of each - // bin and y is the value of each bin). Then we figure out what value - // of x corresponds to f(x) being the LQ/UQ. - if (stats->options & PS_STAT_ROBUST_QUARTILE) { - countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; - - //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. - //Determine id that is okay. - stats->robustLQ = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - cumulativeRobustSums, - LQBinNum, - countFloat/4.0); - //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. - //Determine id that is okay. - stats->robustUQ = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - cumulativeRobustSums, - UQBinNum, - 3.0 * countFloat/4.0); - } - // XXX: I think sumNfit == sumN50 here. - stats->robustNfit = sumNfit; - stats->robustN50 = sumN50; - - psFree(tmpStats); - psFree(robustHistogram); - psFree(robustHistogramVector); - psFree(cumulativeRobustSums); - - return(0); -} - - -/***************************************************************************** -XXX: Is there a psLib function for this? - *****************************************************************************/ -psVector *PsVectorDup(psVector *in) -{ - psVector *out = psVectorAlloc(in->n, in->type.type); - - if (in->type.type == PS_TYPE_F32) { - for (psS32 i = 0 ; i < in->n ; i++) { - out->data.F32[i] = in->data.F32[i]; - } - } else if (in->type.type == PS_TYPE_F64) { - for (psS32 i = 0 ; i < in->n ; i++) { - out->data.F64[i] = in->data.F64[i]; - } - } else { - psError(PS_ERR_UNKNOWN, true, "Unallowable vector type.\n"); - return(NULL); - } - return(out); -} - -/****************************************************************************** -XXX: This function need to be written. Actually, it simply needs to be -retrieved from the CVS repository, since it was written earlier, then -discarded. At present, it was deleted from the CVS repository, so we might -have to retrieve it from tape. -*****************************************************************************/ -psVector *Fit1DGaussian(psVector *x, psVector*y) -{ - psError(PS_ERR_UNKNOWN, true, "This code has not been implemented yet.\n"); - // XXX: This function was previously part of psStats.c, was removed, was - // purged from CVS, and now needs to be retrieved from tape. - - return(NULL); -} + /****************************************************************************** @@ -1984,4 +1681,10 @@ psStats* stats) { + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(myVector); + } + // psS32 currentId = psMemGetId(); + // psS32 memLeaks = 0; psHistogram *robustHistogram = NULL; psHistogram *cumulativeRobustHistogram = NULL; @@ -1991,14 +1694,27 @@ psS32 totalDataPoints = 0; psS32 rc = 0; - psVector *tmpMaskVec = PsVectorDup((psVector *) maskVector); - - while (1) { + psS32 rcBool = false; + psVector *tmpMaskVec = psVectorAlloc(myVector->n, PS_TYPE_U8); + if (maskVector != NULL) { + for (psS32 i = 0 ; i < myVector->n ; i++) { + tmpMaskVec->data.U8[i] = maskVector->data.U8[i]; + } + } else { + for (psS32 i = 0 ; i < myVector->n ; i++) { + tmpMaskVec->data.U8[i] = 0; + } + } + psBool iterate = true; + psF32 sigma; + psStats *tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); + + while (iterate) { + psTrace(__func__, 6, "Iterating on Bin size.\n"); // // Determine the bin size of the robust histogram. This is done // by computing the total range of data values and dividing by 1000.0. // - psStats* tmpStatsMinMax = psStatsAlloc(PS_STAT_MIN | PS_STAT_MAX); - rc = p_psVectorMin(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); - rc|= p_psVectorMax(myVector, tmpMaskVec, maskVal, tmpStatsMinMax); + rc = p_psVectorMin(myVector, tmpMaskVec, 1, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, 1, tmpStatsMinMax); if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); @@ -2006,7 +1722,10 @@ psFree(tmpMaskVec); psFree(tmpScalar); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); return(1); } + psTrace(__func__, 6, "Data min/man is (%.2f, %.2f)\n", tmpStatsMinMax->min, tmpStatsMinMax->max); psF32 binSize = (tmpStatsMinMax->max - tmpStatsMinMax->min) / 1000.0f; + psTrace(__func__, 6, "Robust bin size is %.2f\n", binSize); // @@ -2015,4 +1734,5 @@ // if (fabs(tmpStatsMinMax->max - tmpStatsMinMax->min) <= FLT_EPSILON) { + psTrace(__func__, 5, "All data points have the same value.\n", binSize); if (stats->options & PS_STAT_ROBUST_MEDIAN) { stats->robustMedian = tmpStatsMinMax->min; @@ -2029,4 +1749,5 @@ psFree(tmpScalar); + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); return(0); } @@ -2041,10 +1762,14 @@ // numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / binSize); + psTrace(__func__, 6, "Numbins is %d\n", numBins); robustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); cumulativeRobustHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); // Populate the histogram array. - psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal); - + robustHistogram = psVectorHistogram(robustHistogram, myVector, errors, tmpMaskVec, maskVal); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(robustHistogram->bounds); + PS_VECTOR_PRINT_F32(robustHistogram->nums); + } // // ADD: Step 1. @@ -2056,4 +1781,8 @@ robustHistogram->nums->data.F32[i]; } + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->bounds); + PS_VECTOR_PRINT_F32(cumulativeRobustHistogram->nums); + } // @@ -2062,14 +1791,25 @@ // totalDataPoints = cumulativeRobustHistogram->nums->data.F32[numBins - 1]; - tmpScalar->data.F32 = totalDataPoints/2.0; - psS32 binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); - if (binMedian != 0) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50% data point.\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpScalar); - return(1); - } + psTrace(__func__, 6, "Total data points is %d\n", totalDataPoints); + psS32 binMedian; + if (totalDataPoints/2.0 < cumulativeRobustHistogram->nums->data.F32[0]) { + // Special case: the median is in the first bin. Perhaps we should recode this. + // XXX: Must try a special case where the median in in the last bin. + binMedian = 0; + } else { + tmpScalar->data.F32 = totalDataPoints/2.0; + binMedian = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if (binMedian < 0) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 50 precent data point (%d).\n", binMedian); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(tmpMaskVec); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); + return(1); + } + } + psTrace(__func__, 6, "The median bin is %d\n", binMedian); // @@ -2079,8 +1819,9 @@ // stats->robustMedian = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - *(psVector* *)&robustHistogram->nums, + *(psVector* *)&cumulativeRobustHistogram->bounds, + *(psVector* *)&cumulativeRobustHistogram->nums, binMedian, totalDataPoints/2.0); + psTrace(__func__, 6, "Current robust median is %f\n", stats->robustMedian); // @@ -2088,14 +1829,38 @@ // Find the bins which contains the 15.8655% and 84.1345% data points. // + psS32 binLo = 0; + psS32 binHi = 0; tmpScalar->data.F32 = totalDataPoints * 0.158655f; - psS32 binLo = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + binLo = 0; + } else { + // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. + binLo = 1 + p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if (binLo > cumulativeRobustHistogram->nums->n-1) { + binLo = cumulativeRobustHistogram->nums->n-1; + } + } tmpScalar->data.F32 = totalDataPoints * 0.841345f; - psS32 binHi = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); - if ((binLo != 0) || (binHi != 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the15.8655% and 84.1345% data point\n"); + if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + binHi = 0; + } else { + // NOTE: the (+1) here is because of the way p_psVectorBinDisect is defined. + binHi = 1+ p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + if (binHi > cumulativeRobustHistogram->nums->n-1) { + binHi = cumulativeRobustHistogram->nums->n-1; + } + } + psTrace(__func__, 6, "The 15.8655% and 84.1345% data point bins are (%d, %d).\n", binLo, binHi); + psTrace(__func__, 6, "binLo midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binLo)); + psTrace(__func__, 6, "binHi midpoint is %f\n", PS_BIN_MIDPOINT(cumulativeRobustHistogram, binHi)); + + if ((binLo < 0) || (binHi < 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 15.8655 and 84.1345 percent data point\n"); psFree(tmpStatsMinMax); psFree(robustHistogram); psFree(cumulativeRobustHistogram); psFree(tmpScalar); + psFree(tmpMaskVec); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); return(1); } @@ -2105,14 +1870,44 @@ // Interpolate Sigma to find these two positions exactly: these are the 1sigma positions. // - psF32 binLoF32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - *(psVector* *)&robustHistogram->nums, - binLo, - totalDataPoints * 0.158655f); - psF32 binHiF32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - *(psVector* *)&robustHistogram->nums, - binHi, - totalDataPoints * 0.841345f); + // XXX: I am overriding the ADD for now. The following code implements + // the ADD exactly and is having problems fitting a 2nd-order polynomial + // to data y-values suchs as (1, 1, 100). Therefore, I commented out the + // code and am doing simply linear interpolation. + // + psF32 binLoF32; + psF32 binHiF32; + if (0) { + binLoF32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&cumulativeRobustHistogram->bounds, + *(psVector* *)&cumulativeRobustHistogram->nums, + binLo, + totalDataPoints * 0.158655f); + binHiF32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&cumulativeRobustHistogram->bounds, + *(psVector* *)&cumulativeRobustHistogram->nums, + binHi, + totalDataPoints * 0.841345f); + psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32); + } + + if (1) { + psTrace(__func__, 6, "binLo is %d. Nums at that bin and the next are (%.2f, %.2f)\n", binLo, cumulativeRobustHistogram->nums->data.F32[binLo], cumulativeRobustHistogram->nums->data.F32[binLo+1]); + psTrace(__func__, 6, "binHi is %d. Nums at that bin and the next are (%.2f, %.2f)\n", binHi, cumulativeRobustHistogram->nums->data.F32[binHi], cumulativeRobustHistogram->nums->data.F32[binHi+1]); + // XXX: Add checks to ensure that these bins are not the first (seg fault). + binLoF32 = LinInterpolate( + cumulativeRobustHistogram->bounds->data.F32[binLo], + cumulativeRobustHistogram->bounds->data.F32[binLo+1], + cumulativeRobustHistogram->nums->data.F32[binLo-1], + cumulativeRobustHistogram->nums->data.F32[binLo], + totalDataPoints * 0.158655f); + binHiF32 = LinInterpolate( + cumulativeRobustHistogram->bounds->data.F32[binHi], + cumulativeRobustHistogram->bounds->data.F32[binHi+1], + cumulativeRobustHistogram->nums->data.F32[binHi-1], + cumulativeRobustHistogram->nums->data.F32[binHi], + totalDataPoints * 0.841345f); + // XXX: Check for NANs + psTrace(__func__, 6, "The exact 15.8655 and 84.1345 percent data point positions are: (%f, %f)\n", binLoF32, binHiF32); + } // @@ -2120,5 +1915,6 @@ // Determine SIGMA as 1/2 of the distance between these positions. // - psF32 sigma = (binHiF32 - binLoF32) / 2.0; + sigma = (binHiF32 - binLoF32) / 2.0; + psTrace(__func__, 6, "The current sigma is %f.\n", sigma); // @@ -2129,144 +1925,294 @@ // if (sigma < (2 * binSize)) { + psTrace(__func__, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); psF32 medianLo = robustHistogram->bounds->data.F32[binMedian - 25]; psF32 medianHi = robustHistogram->bounds->data.F32[binMedian + 25]; + psTrace(__func__, 6, "Masking data outside (%f %f)\n", medianLo, medianHi); for (psS32 i = 0 ; i < myVector->n ; i++) { if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) { tmpMaskVec->data.U8[i] = 1; - } - } - } else { - // - // ADD: Step 7. - // Find the bins which contains the 25% and 75% data points. - // - tmpScalar->data.F32 = totalDataPoints * 0.25f; - psS32 binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); - tmpScalar->data.F32 = totalDataPoints * 0.75f; - psS32 binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); - if ((binLo25 != 0) || (binHi25 != 0)) { - psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25% and 75% data points\n"); - psFree(tmpStatsMinMax); - psFree(robustHistogram); - psFree(cumulativeRobustHistogram); - psFree(tmpScalar); - return(1); - } - - // - // ADD: Step 8. - // Interpolate to find these two positions exactly: these are the upper - // and lower quartile positions. - // XXX: Check for errors. - // - psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - *(psVector* *)&robustHistogram->nums, - binLo25, - totalDataPoints * 0.25f); - psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( - *(psVector* *)&robustHistogram->bounds, - *(psVector* *)&robustHistogram->nums, - binHi25, - totalDataPoints * 0.75f); - - stats->robustLQ = binLo25F32; - stats->robustUQ = binHi25F32; - - // PS_BIN_MIDPOINT(robustHistogram, modeBinNum); - // XXX: I think sumNfit == sumN50 here. - stats->robustNfit = -1; - stats->robustN50 = -1; - - // - // Perform the Robust Histogram Statistics algorithm above - // - - // - // Smooth the resulting histogram with a Gaussian with SIGMA_s = 1 - // bin. - // - // XXX: SIGMA_s is defined nowhere in the ADD. - // - psF32 SIGMA_S = 1.0; - p_psVectorSmoothHistGaussian(robustHistogram, SIGMA_S); - - // - // Find the bin with the peak value in the range 2 SIGMA of the - // robust histogram median. - // - // XXX: SIGMA is defined nowhere in the ADD. - // - psF32 SIGMA = 2.0; - psS32 binMin = binMedian - (SIGMA * PS_GAUSS_WIDTH); - if (binMin < 0) { - binMin = 0; - } - psS32 binMax = binMedian + (2 * PS_GAUSS_WIDTH); - if (binMin > (robustHistogram->nums->n - 1)) { - binMin = (robustHistogram->nums->n - 1); - } - psS32 binNum = binNum; - psF32 binMaxNums = robustHistogram->nums->data.F32[binNum]; - for (psS32 i = binNum+1 ; i <= binMax ; i++) { - if (robustHistogram->nums->data.F32[i] > binMaxNums) { - binNum = i; - binMaxNums = robustHistogram->nums->data.F32[i]; - } - } - - // - // Fit a Gaussian to the bins in the range 2 SIGMA of the robust - // histogram median. - // - // XXX: SIGMA is defined nowhere in the ADD. - // - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psVector *x = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); - psS32 j = 0; - for (psS32 i = binNum ; i <= binMax ; i++) { - y->data.F32[j] = robustHistogram->nums->data.F32[i]; - x->data.F32[j] = PS_BIN_MIDPOINT(robustHistogram, i); - } - // - // XXX: This function need to be written. Actually, it simply - // needs to be retrieved from the CVS repository, since it was - // written earlier, then discarded. At present, it was deleted - // from the CVS repository, so we might have to retrieve it from - // tape. - // - psVector *params = Fit1DGaussian(x, y); - - // - // The robust mean mean_r is derived directly from the fitted - // Gaussian mean. - // - stats->robustMean = params->data.F32[0]; - - // - // The robust standard deviation, SIGMA_r is determined by - // subtracting the smoothing scale in quadrature: SIGMA_r^2 = SIGMA^2 - // - SIGMA_s^2 - // - // XXX: SIGMA and SIGMA_s are defined nowhere in the ADD. We must figure - // out what they are. - // - stats->robustStdev = sqrt(PS_SQR(SIGMA) - PS_SQR(SIGMA_S)); - - psFree(tmpStatsMinMax); + psTrace(__func__, 6, "Masking element %d is %f\n", i, myVector->data.F32[i]); + } + } psFree(robustHistogram); psFree(cumulativeRobustHistogram); - psFree(tmpScalar); - psFree(params); - - return(0); - } - + + } else { + psTrace(__func__, 6, "*************: No more iteration. sigma is %f\n", sigma); + iterate = false; + } + } + + // + // ADD: Step 7. + // Find the bins which contains the 25% and 75% data points. + // + psS32 binLo25 = 0; + psS32 binHi25 = 0; + + tmpScalar->data.F32 = totalDataPoints * 0.25f; + if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + // XXX: Special case where its in first bin. Must code last bin. + binLo25 = 0; + } else { + binLo25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + } + tmpScalar->data.F32 = totalDataPoints * 0.75f; + if (tmpScalar->data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) { + // XXX: Special case where its in first bin. Must code last bin. + binHi25 = 0; + } else { + binHi25 = p_psVectorBinDisect(cumulativeRobustHistogram->nums, tmpScalar); + } + if ((binLo25 < 0) || (binHi25 < 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the 25 and 75 percent data points\n"); psFree(tmpStatsMinMax); psFree(robustHistogram); psFree(cumulativeRobustHistogram); - } - return(1); + psFree(tmpScalar); + psFree(tmpMaskVec); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); + return(1); + } + psTrace(__func__, 6, "The 25% and 75% data point bins are (%d, %d).\n", binLo25, binHi25); + + // + // ADD: Step 8. + // Interpolate to find these two positions exactly: these are the upper + // and lower quartile positions. + // XXX: Check for errors. + // + psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&cumulativeRobustHistogram->bounds, + *(psVector* *)&cumulativeRobustHistogram->nums, + binLo25, + totalDataPoints * 0.25f); + psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX( + *(psVector* *)&cumulativeRobustHistogram->bounds, + *(psVector* *)&cumulativeRobustHistogram->nums, + binHi25, + totalDataPoints * 0.75f); + stats->robustLQ = binLo25F32; + stats->robustUQ = binHi25F32; + psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32); + + // + // + // New algorithm for calculated mean and stdev. + // + // + psS32 N50 = 0; + for (psS32 i = 0 ; i < myVector->n ; i++) { + if ((0 == tmpMaskVec->data.U8[i]) && + (binLo25F32 <= myVector->data.F32[i]) && + (binHi25F32 >= myVector->data.F32[i])) { + N50++; + } + } + stats->robustN50 = N50; + psTrace(__func__, 6, "The robustN50 is %d.\n", N50); + + psF32 dN = (psF32) (0.17 * N50); + if (dN < 1.0) { + dN = 1.0; + } else if (dN > 4.0) { + dN = 4.0; + } + psF32 newBinSize = sigma / dN; + + rc = p_psVectorMin(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); + rc|= p_psVectorMax(myVector, tmpMaskVec, 1 | maskVal, tmpStatsMinMax); + if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(tmpMaskVec); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); + return(1); + } + + numBins = (psS32)((tmpStatsMinMax->max - tmpStatsMinMax->min) / newBinSize); + psTrace(__func__, 6, "The new min/max values are (%f, %f).\n", tmpStatsMinMax->min, tmpStatsMinMax->max); + psTrace(__func__, 6, "The new bin size is %f.\n", newBinSize); + psTrace(__func__, 6, "The numBins is %d\n", numBins); + + psHistogram *newHistogram = psHistogramAlloc(tmpStatsMinMax->min, tmpStatsMinMax->max, numBins); + newHistogram = psVectorHistogram(newHistogram, myVector, errors, tmpMaskVec, maskVal|1); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(newHistogram->nums); + } + + // + // Smooth the resulting histogram with a Gaussian with sigma_s = 1 + // bin. + // + psF32 sigma_s = newBinSize; + + psVector *newHistogramSmoothed = p_psVectorSmoothHistGaussian(newHistogram, sigma_s); + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(newHistogramSmoothed); + } + + // + // Find the bin with the peak value in the range 2 sigma of the + // robust histogram median. + // + + psS32 binMin = 0; + psS32 binMax = 0; + tmpScalar->data.F32 = stats->robustMedian - (2.0 * sigma); + if (tmpScalar->data.F32 <= newHistogram->bounds->data.F32[0]) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + } + + tmpScalar->data.F32 = stats->robustMedian + (2.0 + sigma); + if (tmpScalar->data.F32 >= newHistogram->bounds->data.F32[newHistogram->bounds->n-1]) { + binMax = newHistogram->bounds->n-1; + } else { + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + } + if ((binMin < 0) || (binMax < 0)) { + psError(PS_ERR_UNKNOWN, false, "Failed to calculate the +- 2.0 * sigma bins\n"); + psFree(tmpMaskVec); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpStatsMinMax); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); + return(1); + } + + psS32 binNum = binMin; + psF32 binMaxNums = newHistogramSmoothed->data.F32[binNum]; + for (psS32 i = binMin+1 ; i <= binMax ; i++) { + if (newHistogramSmoothed->data.F32[i] > binMaxNums) { + binNum = i; + binMaxNums = newHistogramSmoothed->data.F32[i]; + } + } + psTrace(__func__, 6, "The peak bin is %d, with %f data.n", binNum, binMaxNums); + + // + // Fit a Gaussian to the bins in the range 20 sigma of the robust + // histogram median. + // + tmpScalar->data.F32 = stats->robustMedian - (20.0 * sigma); + if (tmpScalar->data.F32 < tmpStatsMinMax->min) { + binMin = 0; + } else { + binMin = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + // XXX: check for errors here. + } + tmpScalar->data.F32 = stats->robustMedian + (20.0 * sigma); + if (tmpScalar->data.F32 > tmpStatsMinMax->max) { + binMax = newHistogramSmoothed->n - 1; + } else { + binMax = p_psVectorBinDisect((psVector *) newHistogram->bounds, tmpScalar); + // XXX: check for errors here. + } + psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psVector *xTmp = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); + psArray *x = psArrayAlloc((1 + (binMax - binMin))); + stats->robustNfit = 0; + psS32 j = 0; + + for (psS32 i = binMin ; i <= binMax ; i++) { + y->data.F32[j] = newHistogramSmoothed->data.F32[i]; + x->data[j] = (psPtr *) psVectorAlloc(1, PS_TYPE_F32); + ((psVector *) x->data[j])->data.F32[0] = PS_BIN_MIDPOINT(newHistogram, i); + xTmp->data.F32[j] = PS_BIN_MIDPOINT(newHistogram, i); + + stats->robustNfit+= newHistogramSmoothed->data.F32[i]; + j++; + } + if (psTraceGetLevel(__func__) >= 8) { + // XXX: Print the x array somehow. + PS_VECTOR_PRINT_F32(y); + } + + // XXX: Use the min/max routines for this + psF32 minY = FLT_MAX; + psF32 maxY = -FLT_MAX; + for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + if (y->data.F32[i] > maxY) { + maxY = y->data.F32[i]; + } + if (y->data.F32[i] < minY) { + minY = y->data.F32[i]; + } + } + // + // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) + // XXX: Use the normalize routines for this. + // + for (psS32 i = 0 ; i < 1 + (binMax - binMin) ; i++) { + y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY); + } + + // + psMinimization *min = psMinimizationAlloc(100, 0.01); + psVector *params = psVectorAlloc(2, PS_TYPE_F32); + // Initial guess for the mean ([0]) and standard dev. + params->data.F32[0] = stats->robustMedian; + params->data.F32[1] = sigma; + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(params); + PS_VECTOR_PRINT_F32(y); + } + psFree(xTmp); + rcBool = psMinimizeLMChi2(min, NULL, params, NULL, x, y, NULL, psMinimizeLMChi2Gauss1D); + + if (rcBool != true) { + psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(tmpStatsMinMax); + psFree(robustHistogram); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(newHistogram); + psFree(x); + psFree(y); + psFree(min); + psFree(params); + psFree(tmpMaskVec); + psTrace(__func__, 3, "---- %s(1) end ----\n", __func__); + return(1); + } + if (psTraceGetLevel(__func__) >= 8) { + PS_VECTOR_PRINT_F32(params); + } + + // + // The robust mean mean_r is derived directly from the fitted + // Gaussian mean. + // + stats->robustMean = params->data.F32[0]; + psTrace(__func__, 6, "The robust mean is %f.\n", params->data.F32[0]); + + // + // The robust standard deviation, SIGMA_r is determined by + // subtracting the smoothing scale in quadrature: + // SIGMA_r^2 = SIGMA^2 - sigma_s^2 + // + stats->robustStdev = sqrt(PS_SQR(params->data.F32[1]) - PS_SQR(sigma_s)); + psTrace(__func__, 6, "The robust stdev is %f.\n", stats->robustStdev); + + psFree(newHistogramSmoothed); + psFree(tmpStatsMinMax); + psFree(cumulativeRobustHistogram); + psFree(tmpScalar); + psFree(newHistogram); + psFree(x); + psFree(y); + psFree(min); + psFree(params); + psFree(tmpMaskVec); + psFree(robustHistogram); + + psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); + return(0); } @@ -2808,5 +2754,5 @@ // Since the various robust stats quantities share much computation, they // are grouped together in a single private function: - // p_psVectorRobustStats() + // p_psVectorRobustStatsNew() if ((stats->options & PS_STAT_ROBUST_MEAN) || (stats->options & PS_STAT_ROBUST_MEDIAN) || @@ -2814,8 +2760,12 @@ (stats->options & PS_STAT_ROBUST_STDEV) || (stats->options & PS_STAT_ROBUST_QUARTILE)) { - if (0 != p_psVectorRobustStats(inF32, errorsF32, mask, maskVal, stats)) { + if (0 != p_psVectorRobustStatsNew(inF32, errorsF32, mask, maskVal, stats)) { psError(PS_ERR_UNKNOWN, false, PS_ERRORTEXT_psStats_STATS_FAILED); // XXX: Set to NAN + // XXX: Is this the right thing to do? + // XXX: If so, do it for other funcs? + psFree(stats); + return(NULL); } }