Index: /branches/eam_branches/ipp-20110213/psLib/src/math/psMinimizePolyFit.c =================================================================== --- /branches/eam_branches/ipp-20110213/psLib/src/math/psMinimizePolyFit.c (revision 30862) +++ /branches/eam_branches/ipp-20110213/psLib/src/math/psMinimizePolyFit.c (revision 30863) @@ -767,6 +767,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); @@ -1211,6 +1211,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); @@ -1621,6 +1621,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); @@ -2055,6 +2055,6 @@ // XXX enforce consistency? // XXX psStatsGetValue() probably has inverted precedence - psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); - psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV); if (!meanOption) { psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); Index: /branches/eam_branches/ipp-20110213/psLib/src/math/psStats.c =================================================================== --- /branches/eam_branches/ipp-20110213/psLib/src/math/psStats.c (revision 30862) +++ /branches/eam_branches/ipp-20110213/psLib/src/math/psStats.c (revision 30863) @@ -172,8 +172,4 @@ *****************************************************************************/ -// static prototypes: -static psF32 minimizeLMChi2Gauss1D(psVector *deriv, const psVector *params, const psVector *coords); -// static psF32 fitQuadraticSearchForYThenReturnX(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); -// static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); @@ -1087,13 +1083,16 @@ } -/* - * vectorFittedStats requires guess for fittedMean and fittedStdev - * robustN50 should also be set - */ +/******************** + * perform an asymmetric fit to the population. In development, this was called + * "vectorFittedStats_v4" all versions of fitted stats now resolve to this function (only v4 + * has really been used) vectorFittedStats requires guess for fittedMean and fittedStdev + * robustN50 should also be set gaussian fit is performed using 2D polynomial to ln(y) this + * 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) + const psVector* errors, + psVector* mask, + psVectorMaskType maskVal, + psStats* stats) { @@ -1121,401 +1120,4 @@ stats->results |= PS_STAT_FITTED_MEAN; stats->results |= PS_STAT_FITTED_STDEV; - return true; - } - - float guessStdev = stats->robustStdev; // pass the guess sigma - float guessMean = stats->robustMedian; // pass the guess mean - - psTrace(TRACE, 6, "The guess mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The guess stdev is %f.\n", guessStdev); - - bool done = false; - 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) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { - // construct a histogram with (sigma/10 < binsize < sigma) - // set roughly so that the lowest bins have about 2 cnts - // Nsmallest ~ N50 / (4*dN)) - psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); - binSize = guessStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->max; - if (numValid == 0 || isnan(min) || isnan(max)) { - psTrace(TRACE, 5, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - return true; - } - - // Calculate the number of bins. - // XXX can we calculate the binMin, binMax **before** building this histogram? - long numBins = (max - min) / binSize; - 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); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - // if psVectorHistogram returns false, we have a programming error - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // Fit a Gaussian to the bins in the requested range about the guess mean - // min,max overrides clipSigma - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - // 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); - - // Generate the variables that will be used in the Gaussian fitting - // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins - psVector *y = psVectorAlloc((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates - psArray *x = psArrayAlloc((1 + (binMax - binMin))); // Array of ordinates - for (long i = binMin, j = 0; i <= binMax ; i++, j++) { - y->data.F32[j] = histogram->nums->data.F32[i]; - psVector *ordinate = psVectorAlloc(1, PS_TYPE_F32); // The ordinate value - ordinate->data.F32[0] = PS_BIN_MIDPOINT(histogram, i); - x->data[j] = ordinate; - } - if (psTraceGetLevel("psLib.math") >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - psTrace(TRACE, 6, "The clipped numBins is %ld\n", y->n); - psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); - - // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0]) - // XXX EAM : why not just fit the normalization as well?? - p_psNormalizeVectorRange(y, 0.0, 1.0); - - // Fit a Gaussian to the data. - psMinimization *minimizer = psMinimizationAlloc(100, 0.01, 1.0); // The minimizer information - psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian - // Initial guess for the mean (index 0) and var (index 1). - params->data.F32[0] = guessMean; - params->data.F32[1] = PS_SQR(guessStdev); - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(params); - PS_VECTOR_PRINT_F32(y); - } - - // psMinimizeLMChi2 can return false for bad data as well as for serious failures - if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(params); - psFree(minimizer); - psFree(x); - psFree(y); - psFree(histogram); - psFree(statsMinMax); - return true; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(params); - } - - guessMean = params->data.F32[0]; - guessStdev = sqrt(params->data.F32[1]); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - - // Clean up after fitting - psFree (params); - psFree (minimizer); - psFree (x); - psFree (y); - psFree (histogram); - psFree (statsMinMax); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = guessMean; - psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = guessStdev; - psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - stats->results |= PS_STAT_FITTED_MEAN; - stats->results |= PS_STAT_FITTED_STDEV; - - return true; -} - - -/******************** - * vectorFittedStats_v2 requires guess for fittedMean and fittedStdev - * robustN50 should also be set - * gaussian fit is performed using 2D polynomial to ln(y) - ********************/ -static bool vectorFittedStats_v2 (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)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } - } - - // 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_V2; - stats->results |= PS_STAT_FITTED_STDEV_V2; - return true; - } - - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN_V2; - stats->results |= PS_STAT_FITTED_STDEV_V2; - return true; - } - - 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); - psTrace(TRACE, 6, "The ** starting ** guess stdev is %f.\n", guessStdev); - - bool done = false; - 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) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } else { - // construct a histogram with (sigma/10 < binsize < sigma) - // set roughly so that the lowest bins have about 2 cnts - // Nsmallest ~ N50 / (4*dN)) - psF32 dN = PS_MAX (1, PS_MIN (10, stats->robustN50 / 8)); - binSize = guessStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->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); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return true; - } - - // Calculate the number of bins. - // XXX can we calculate the binMin, binMax **before** building this histogram? - long numBins = (max - min) / binSize; - 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); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistcs v2.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // Fit a Gaussian to the bins in the requested range about the guess mean - // min,max overrides clipSigma - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - // 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); - - // Generate the variables that will be used in the Gaussian fitting - // XXX EAM : we should test / guarantee that there are at least 3 (right?) bins - psVector *y = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of coordinates - psVector *x = psVectorAllocEmpty((1 + (binMax - binMin)), PS_TYPE_F32); // Vector of ordinates - long j = 0; - for (long i = binMin; i <= binMax ; i++) { - if (histogram->nums->data.F32[i] <= 0.0) - continue; - x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); - y->data.F32[j] = log(histogram->nums->data.F32[i]); - // XXX note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) - j++; - } - y->n = x->n = j; - if (psTraceGetLevel("psLib.math") >= 8) { - // XXX: Print the x array somehow. - PS_VECTOR_PRINT_F32(y); - } - psTrace(TRACE, 6, "The clipped numBins is %ld\n", y->n); - psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax), binMax); - - // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 - psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - - // XXX how can we protect against some extreme outliers? the robust histogram - // should have already dealt with those, but we are again sensitive to them... - // psStats *fitStats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - // fitStats->clipIter = 3.0; - // fitStats->clipSigma = 3.0; - // psVector *fitMask = psVectorAlloc(y->n, PS_TYPE_VECTOR_MASK); - // psVectorInit (fitMask, 0); - - // XXX not sure if these should result in errors or not... - if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(x); - psFree(y); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - if (poly->coeff[2] >= 0.0) { - psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psError(PS_ERR_UNKNOWN, false, "fit did not converge\n"); - psFree(x); - psFree(y); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } else { - psTrace(TRACE, 6, "Parabolic fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The new guess mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The new guess stdev is %f.\n", guessStdev); - } - - // Clean up after fitting - psFree (x); - psFree (y); - psFree (poly); - // psFree (fitStats); - // psFree (fitMask); - psFree (histogram); - psFree (statsMinMax); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = guessMean; - psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = guessStdev; - psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - stats->results |= PS_STAT_FITTED_MEAN_V2; - stats->results |= PS_STAT_FITTED_STDEV_V2; - - return true; -} - -/******************** - * perform an asymmetric fit to the population - * vectorFittedStats_v3 requires guess for fittedMean and fittedStdev - * robustN50 should also be set - * gaussian fit is performed using 2D polynomial to ln(y) - ********************/ -static bool vectorFittedStats_v3 (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)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } - } - - // 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_V3; - stats->results |= PS_STAT_FITTED_STDEV_V3; - return true; - } - - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN_V3; - stats->results |= PS_STAT_FITTED_STDEV_V3; return true; } @@ -1551,310 +1153,4 @@ COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return true; - } - - // Calculate the number of bins. - // XXX can we calculate the binMin, binMax **before** building this histogram? - long numBins = (max - min) / binSize; - 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); - - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) - if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics v3.\n"); - psFree(histogram); - psFree(statsMinMax); - return false; - } - if (psTraceGetLevel("psLib.math") >= 8) { - PS_VECTOR_PRINT_F32(histogram->nums); - } - - // now fit a Gaussian to the upper and lower halves about the peak independently - - // set the full-range upper and lower limits - psF32 maxFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - maxFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->max)) { - maxFitSigma = fabs(stats->max); - } - - psF32 minFitSigma = 2.0; - if (isfinite(stats->clipSigma)) { - minFitSigma = fabs(stats->clipSigma); - } - if (isfinite(stats->min)) { - minFitSigma = fabs(stats->min); - } - - // 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) { - COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n"); - psFree(statsMinMax); - return true; - } - - // search for mode (peak of histogram within range mean-2sigma - mean+2sigma - long binPeak = binMin; - float valPeak = histogram->nums->data.F32[binPeak]; - 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]); - } - psTrace (TRACE, 6, "\n"); - - // 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); - psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); - psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binPeak), binPeak); - psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); - - { - // fit the lower half of the distribution - // run down until we drop below 0.25*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) { - binS = i; - break; - } - } - psTrace(TRACE, 6, "Lower bound for lower half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binS), binS); - psTrace(TRACE, 6, "Upper bound for lower half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binE), binE); - - psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates - psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates - long j = 0; - for (long i = binS; i < binE; i++) { - if (histogram->nums->data.F32[i] <= 0.0) - continue; - x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); - // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) - y->data.F32[j] = log(histogram->nums->data.F32[i]); - j++; - } - y->n = x->n = j; - - // 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); - psFree(x); - psFree(y); - - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - if (poly->coeff[2] >= 0.0) { - psTrace(TRACE, 6, "Failed parabolic fit: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - - // sometimes, the guessStdev is much too large. in this case, the entire real population - // 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; - psTrace(TRACE, 6, "*** retry, new stdev is %f.\n", guessStdev); - continue; - } - - psError(PS_ERR_UNKNOWN, false, "fit did not converge\n"); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - // calculate lower mean & stdev from parabolic fit -- use this as the result - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); - - psFree(poly); - } - - // for test, measure the same result for the upper section - { - // fit the upper half of the distribution - // run up until we drop below 0.25*valPeak - long binS = PS_MAX (binPeak - 3, 0); - long binE = binMax; - for (int i = binPeak+3; i < binMax; i++) { - if (histogram->nums->data.F32[i] < 0.25*valPeak) { - binE = i; - break; - } - } - psTrace(TRACE, 6, "Lower bound for upper half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binS), binS); - psTrace(TRACE, 6, "Upper bound for upper half: %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binE), binE); - - psVector *y = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of coordinates - psVector *x = psVectorAllocEmpty(binE - binS, PS_TYPE_F32); // Vector of ordinates - long j = 0; - for (long i = binS; i < binE; i++) { - if (histogram->nums->data.F32[i] <= 0.0) - continue; - x->data.F32[j] = PS_BIN_MIDPOINT(histogram, i); - // note this is the natural log: expected distribution is A exp(-(x-xo)^2/2sigma^2) - y->data.F32[j] = log(histogram->nums->data.F32[i]); - j++; - } - y->n = x->n = j; - - // 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); - psFree(x); - psFree(y); - - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(poly); - psFree(histogram); - psFree(statsMinMax); - psTrace(TRACE, 4, "---- %s(false) end ----\n", __func__); - return false; - } - - // calculate upper mean & stdev from parabolic fit -- ignore this value - float upperStdev = sqrt(-0.5/poly->coeff[2]); - float upperMean = poly->coeff[1]*PS_SQR(upperStdev); -#ifndef PS_NO_TRACE - psTrace(TRACE, 6, "Parabolic Upper fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The upper mean is %f.\n", upperMean); - psTrace(TRACE, 6, "The upper stdev is %f.\n", upperStdev); -#endif - - // if the resulting value is outside of the range binMin - binMax, use the upper value - if (done && (guessMean > PS_BIN_MIDPOINT(histogram, binMax - 1))) { - guessMean = upperMean; - guessStdev = upperStdev; - } - - psFree (poly); - } - - // Clean up after fitting - psFree (histogram); - psFree (statsMinMax); - } - - // The fitted mean is the Gaussian mean. - stats->fittedMean = guessMean; - psTrace(TRACE, 6, "The fitted mean is %f.\n", stats->fittedMean); - - // The fitted standard deviation - stats->fittedStdev = guessStdev; - psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - - stats->results |= PS_STAT_FITTED_MEAN_V3; - stats->results |= PS_STAT_FITTED_STDEV_V3; - - return true; -} - -/******************** - * perform an asymmetric fit to the population - * vectorFittedStats_v4 requires guess for fittedMean and fittedStdev - * robustN50 should also be set - * gaussian fit is performed using 2D polynomial to ln(y) - * this version follows the upper portion of the distribution until it passes 0.5*peak - ********************/ -static bool vectorFittedStats_v4 (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)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n"); - return false; - } - } - - // 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_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; - return true; - } - - if (stats->robustStdev <= FLT_EPSILON) { - stats->fittedMean = stats->robustMedian; - stats->fittedStdev = stats->robustStdev; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; - return true; - } - - 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); - psTrace(TRACE, 6, "The ** starting ** guess stdev is %f.\n", guessStdev); - - bool done = false; - 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) { - // Set initial bin size to the specified value. - binSize = stats->binsize; - psTrace(TRACE, 6, "Setting initial robust bin size to %.2f\n", binSize); - } 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)); - binSize = guessStdev / dN; - } - - // Determine the min/max of the vector (which prior outliers masked out) - int numValid = vectorMinMax(myVector, mask, maskVal, statsMinMax); // Number of values - float min = statsMinMax->min; - float max = statsMinMax->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); goto escape; } @@ -1865,6 +1161,6 @@ stats->fittedMean = min; stats->fittedStdev = 0.0; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; } @@ -1934,11 +1230,14 @@ // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak + float clippedMean = PS_BIN_MIDPOINT(histogram, binPeak); + psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin); psTrace(TRACE, 6, "The clipped min is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMin), binMin); - psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", - PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); - psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binPeak), binPeak); + psTrace(TRACE, 6, "The clipped max is %f (%ld)\n", PS_BIN_MIDPOINT(histogram, binMax - 1), binMax - 1); + psTrace(TRACE, 6, "The clipped peak is %f (%ld)\n", clippedMean, binPeak); psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); + float lowfitMean = NAN; + float lowfitStdev = NAN; { // fit the lower half of the distribution @@ -2014,24 +1313,21 @@ // calculate lower mean & stdev from parabolic fit -- use this as the result - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); - if (guessStdev > 0.75*stats->robustStdev) { - done = true; - } - psTrace(TRACE, 6, "Parabolic Lower fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The lower mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The lower stdev is %f.\n", guessStdev); + 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]); + psTrace(TRACE, 6, "The lower mean is %f.\n", lowfitMean); + psTrace(TRACE, 6, "The lower stdev is %f.\n", lowfitStdev); psFree(poly); } - // if we converge on a solution outside the range binMin - binMax, use a more conservative range - float minValue = PS_BIN_MIDPOINT(histogram, binMin); - float maxValue = PS_BIN_MIDPOINT(histogram, binMax - 1); - - if (done && ((guessMean < minValue) || (guessMean > maxValue))) { - psTrace(TRACE, 6, "Inconsistent result, re-trying the fit\n"); - + 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 @@ -2085,26 +1381,25 @@ // calculate upper mean & stdev from parabolic fit -- ignore this value - guessStdev = sqrt(-0.5/poly->coeff[2]); - guessMean = poly->coeff[1]*PS_SQR(guessStdev); + fullfitStdev = sqrt(-0.5/poly->coeff[2]); + fullfitMean = poly->coeff[1]*PS_SQR(fullfitStdev); #ifndef PS_NO_TRACE - psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", - poly->coeff[0], poly->coeff[1], poly->coeff[2]); - psTrace(TRACE, 6, "The symmetric mean is %f.\n", guessMean); - psTrace(TRACE, 6, "The symmetric stdev is %f.\n", guessStdev); + psTrace(TRACE, 6, "Parabolic Symmetric fit results: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]); + psTrace(TRACE, 6, "The symmetric mean is %f.\n", fullfitMean); + psTrace(TRACE, 6, "The symmetric stdev is %f.\n", fullfitStdev); #endif // if we converge on a solution outside the range binMin - binMax, use a more conservative range - float minValueSym = PS_BIN_MIDPOINT(histogram, binS); - float maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1); + minValueSym = PS_BIN_MIDPOINT(histogram, binS); + maxValueSym = PS_BIN_MIDPOINT(histogram, binE - 1); // saturate on min or max value - if (guessMean < minValueSym) { - guessMean = minValueSym; + if (fullfitMean < minValueSym) { + fullfitMean = minValueSym; psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); } // saturate on min or max value - if (guessMean > maxValueSym) { - guessMean = maxValueSym; + if (fullfitMean > maxValueSym) { + fullfitMean = maxValueSym; psTrace(TRACE, 6, "The symmetric mean is out of bounds, saturating to %f.\n", guessMean); } @@ -2112,4 +1407,24 @@ 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; + } // Clean up after fitting @@ -2126,6 +1441,6 @@ psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev); - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; @@ -2134,6 +1449,6 @@ stats->fittedMean = NAN; stats->fittedStdev = NAN; - stats->results |= PS_STAT_FITTED_MEAN_V4; - stats->results |= PS_STAT_FITTED_STDEV_V4; + stats->results |= PS_STAT_FITTED_MEAN; + stats->results |= PS_STAT_FITTED_STDEV; return true; @@ -2442,37 +1757,4 @@ // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2)) { - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V2"); - } - if (!vectorFittedStats_v2(inF32, errorsF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); - status &= false; - } - } - - // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_STDEV_V3)) { - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V3"); - } - if (!vectorFittedStats_v3(inF32, errorsF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); - status &= false; - } - } - - // ************************************************************************ - if (stats->options & (PS_STAT_FITTED_MEAN_V4 | PS_STAT_FITTED_STDEV_V4)) { - if (stats->options & (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV)) { - psAbort("you may not specify both FITTED_MEAN and FITTED_MEAN_V4"); - } - if (!vectorFittedStats_v4(inF32, errorsF32, maskVector, maskVal, stats)) { - psError(PS_ERR_UNKNOWN, false, _("Failed to calculate fitted statistics")); - status &= false; - } - } - - // ************************************************************************ if ((stats->options & PS_STAT_CLIPPED_MEAN) || (stats->options & PS_STAT_CLIPPED_STDEV)) { if (!vectorClippedStats(inF32, errorsF32, maskVector, maskVal, stats)) { @@ -2513,16 +1795,16 @@ 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_V2); - READ_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN_V2); - READ_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2); - READ_STAT("FITTED_V3", PS_STAT_FITTED_MEAN_V3); - READ_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN_V3); - READ_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3); - READ_STAT("FITTED_V4", PS_STAT_FITTED_MEAN_V4); - READ_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN_V4); - READ_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4); + 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_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_STDEV_V4", PS_STAT_FITTED_STDEV); READ_STAT("CLIPPED", PS_STAT_CLIPPED_MEAN); READ_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); @@ -2554,10 +1836,4 @@ WRITE_STAT("FITTED_MEAN", PS_STAT_FITTED_MEAN); WRITE_STAT("FITTED_STDEV", PS_STAT_FITTED_STDEV); - WRITE_STAT("FITTED_MEAN_V2", PS_STAT_FITTED_MEAN_V2); - WRITE_STAT("FITTED_STDEV_V2", PS_STAT_FITTED_STDEV_V2); - WRITE_STAT("FITTED_MEAN_V3", PS_STAT_FITTED_MEAN_V3); - WRITE_STAT("FITTED_STDEV_V3", PS_STAT_FITTED_STDEV_V3); - WRITE_STAT("FITTED_MEAN_V4", PS_STAT_FITTED_MEAN_V4); - WRITE_STAT("FITTED_STDEV_V4", PS_STAT_FITTED_STDEV_V4); WRITE_STAT("CLIPPED_MEAN", PS_STAT_CLIPPED_MEAN); WRITE_STAT("CLIPPED_STDEV", PS_STAT_CLIPPED_STDEV); @@ -2610,10 +1886,4 @@ case PS_STAT_FITTED_MEAN: case PS_STAT_FITTED_STDEV: - case PS_STAT_FITTED_MEAN_V2: - case PS_STAT_FITTED_STDEV_V2: - case PS_STAT_FITTED_MEAN_V3: - case PS_STAT_FITTED_STDEV_V3: - case PS_STAT_FITTED_MEAN_V4: - case PS_STAT_FITTED_STDEV_V4: case PS_STAT_CLIPPED_MEAN: case PS_STAT_CLIPPED_STDEV: @@ -2631,6 +1901,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_FITTED_MEAN_V2 | - PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_MEAN_V4); + PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN); } @@ -2638,6 +1907,5 @@ { return options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | - PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2 | PS_STAT_FITTED_STDEV_V3 | - PS_STAT_FITTED_STDEV_V4); + PS_STAT_FITTED_STDEV); } @@ -2665,16 +1933,4 @@ case PS_STAT_FITTED_STDEV: return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V2: - return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V2: - return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V3: - return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V3: - return stats->fittedStdev; - case PS_STAT_FITTED_MEAN_V4: - return stats->fittedMean; - case PS_STAT_FITTED_STDEV_V4: - return stats->fittedStdev; case PS_STAT_CLIPPED_MEAN: return stats->clippedMean; @@ -3115,37 +2371,2 @@ return tmpFloat; } - -/****************************************************************************** -NOTE: We assume unnormalized gaussians. -*****************************************************************************/ -static psF32 minimizeLMChi2Gauss1D(psVector *deriv, - const psVector *params, - const psVector *coords - ) -{ - psTrace(TRACE, 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_VECTOR_NON_NULL(params, NAN); - PS_ASSERT_VECTOR_SIZE(params, (long)2, NAN); - PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_NON_NULL(coords, NAN); - PS_ASSERT_VECTOR_SIZE(coords, (long)1, NAN); - PS_ASSERT_VECTOR_TYPE(coords, PS_TYPE_F32, NAN); - - psF32 x = coords->data.F32[0]; - psF32 mean = params->data.F32[0]; - psF32 var = params->data.F32[1]; - psF32 dx = (x - mean); - - psF32 gauss = exp (-0.5*PS_SQR(dx)/var); - if (deriv) { - PS_ASSERT_VECTOR_SIZE(deriv, (long)2, NAN); - PS_ASSERT_VECTOR_TYPE(deriv, PS_TYPE_F32, NAN); - psF32 tmp = dx * gauss; - deriv->data.F32[0] = tmp / var; - deriv->data.F32[1] = tmp * dx / (var * var); - } - - - psTrace(TRACE, 4, "---- %s() end ----\n", __func__); - return gauss; -} Index: /branches/eam_branches/ipp-20110213/psLib/src/math/psStats.h =================================================================== --- /branches/eam_branches/ipp-20110213/psLib/src/math/psStats.h (revision 30862) +++ /branches/eam_branches/ipp-20110213/psLib/src/math/psStats.h (revision 30863) @@ -43,14 +43,8 @@ PS_STAT_FITTED_MEAN = 0x001000, ///< Fitted Mean PS_STAT_FITTED_STDEV = 0x002000, ///< Fitted Standard Deviation - PS_STAT_FITTED_MEAN_V2 = 0x004000, ///< Fitted Mean - PS_STAT_FITTED_STDEV_V2 = 0x008000, ///< Fitted Standard Deviation - PS_STAT_FITTED_MEAN_V3 = 0x010000, ///< Fitted Mean - PS_STAT_FITTED_STDEV_V3 = 0x020000, ///< 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 - PS_STAT_FITTED_MEAN_V4 = 0x400000, ///< Fitted Mean - PS_STAT_FITTED_STDEV_V4 = 0x800000, ///< Fitted Standard Deviation } psStatsOptions; Index: /branches/eam_branches/ipp-20110213/psLib/test/math/tap_psStats_Sample_01.c =================================================================== --- /branches/eam_branches/ipp-20110213/psLib/test/math/tap_psStats_Sample_01.c (revision 30862) +++ /branches/eam_branches/ipp-20110213/psLib/test/math/tap_psStats_Sample_01.c (revision 30863) @@ -586,5 +586,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); @@ -622,5 +622,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); @@ -657,5 +657,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); @@ -694,5 +694,5 @@ psFree (stats); - stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2 | PS_STAT_USE_BINSIZE); + stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV | PS_STAT_USE_BINSIZE); stats->binsize = 1.0; psVectorStats (stats, y, NULL, NULL, 1); Index: /branches/eam_branches/ipp-20110213/psLib/test/optime/tap_psStatsTiming.c =================================================================== --- /branches/eam_branches/ipp-20110213/psLib/test/optime/tap_psStatsTiming.c (revision 30862) +++ /branches/eam_branches/ipp-20110213/psLib/test/optime/tap_psStatsTiming.c (revision 30863) @@ -678,5 +678,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *rnd2 = psVectorAlloc (1000, PS_TYPE_F32); for (int i = 0; i < rnd2->n; i++) @@ -702,5 +702,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *rnd2 = psVectorAlloc (3000, PS_TYPE_F32); for (int i = 0; i < rnd2->n; i++) @@ -725,5 +725,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *rnd2 = psVectorAlloc (10000, PS_TYPE_F32); for (int i = 0; i < rnd2->n; i++) @@ -790,5 +790,5 @@ psMemId id = psMemGetId(); - psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_STDEV_V2); + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); psVector *sample = psVectorAlloc (1000, PS_TYPE_F32); psVector *robust = psVectorAlloc (1000, PS_TYPE_F32);