Index: /branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c (revision 36261) +++ /branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c (revision 36262) @@ -140,6 +140,4 @@ } -// Debug information -#define CZW 0 /*****************************************************************************/ @@ -174,6 +172,5 @@ *****************************************************************************/ -// static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); -static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); +static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); /****************************************************************************** @@ -232,5 +229,4 @@ } count++; - } if (errors) { @@ -797,5 +793,5 @@ } else { // Determine the bin size of the robust histogram, using the pre-defined number of bins - binSize = (max - min) / INITIAL_NUM_BINS; + binSize = (max - min) / INITIAL_NUM_BINS; } psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize); @@ -880,22 +876,8 @@ cumulative = psHistogramAlloc(min, max, numBins); cumulative->nums->data.F32[0] = histogram->nums->data.F32[0]; - cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1]; - - // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins). - long Nc = 1; // track the current bin of cumulative - // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the - // upper end of the current histogram bin - for (long i = 1; i < histogram->nums->n - 1; i++) { - if (histogram->nums->data.F32[i] == 0.0) continue; - cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i]; - cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1]; - Nc ++; - } - long Nlast = Nc - 1; // last valid cumulative bin - for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled. - cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast]; - cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0; - } - + for (long i = 1; i < histogram->nums->n; i++) { + cumulative->nums->data.F32[i] = cumulative->nums->data.F32[i-1] + histogram->nums->data.F32[i]; + cumulative->bounds->data.F32[i-1] = histogram->bounds->data.F32[i]; + } if (psTraceGetLevel("psLib.math") >= 8) { PS_VECTOR_PRINT_F32(cumulative->bounds); @@ -913,10 +895,7 @@ // ADD step 3: Interpolate to the exact 50% position in bin units - // stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); - // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + stats->robustMedian = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); + // float robustBin = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian); - // There's no reason to do a quadratic fit near the 50% bin, as it's approximately linear there. - // Instead, do a 5-point linear fit. - stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0); // convert bin to bin value: this is the robust histogram median. @@ -933,8 +912,7 @@ PS_BIN_FOR_VALUE(binL2, cumulative->nums, totalDataPoints * 0.308538f, 0); PS_BIN_FOR_VALUE(binH2, cumulative->nums, totalDataPoints * 0.691462f, 0); - PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0); - PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0); - - + PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022481f, 0); + PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977519f, 0); + psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n", binLo, binHi); @@ -948,5 +926,5 @@ goto escape; } - + // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma // positions. @@ -960,5 +938,4 @@ // (extrapolation should not be needed and will result in errors) float binLoF32, binHiF32, binL2F32, binH2F32, binL4F32, binH4F32; -#if (0) PS_BIN_INTERPOLATE (binLoF32, cumulative->nums, cumulative->bounds, binLo, totalDataPoints * 0.158655f); @@ -970,15 +947,8 @@ totalDataPoints * 0.691462f); PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4, - totalDataPoints * 0.022750f); + totalDataPoints * 0.022481f); PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, - totalDataPoints * 0.977250f); -#else - binLoF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo, totalDataPoints * 0.158655); - binHiF32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi, totalDataPoints * 0.841345); - binL2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL2, totalDataPoints * 0.308538); - binH2F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH2, totalDataPoints * 0.691462); - binL4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binL4, totalDataPoints * 0.022750); - binH4F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binH4, totalDataPoints * 0.977250); -#endif + totalDataPoints * 0.977519f); + // report +/- 1 sigma points psTrace(TRACE, 5, @@ -989,39 +959,16 @@ binL2F32, binH2F32); psTrace(TRACE, 5, - "The exact 02.2275 and 97.7250 percent data point positions are: (%f, %f)\n", + "The exact 02.22481 and 97.7519 percent data point positions are: (%f, %f)\n", binL4F32, binH4F32); - // If some of the fits failed, attempt to fix this - if (!isfinite(binLoF32) && isfinite(binHiF32)) { binLoF32 = -1.0 * binHiF32; } - if (!isfinite(binHiF32) && isfinite(binLoF32)) { binHiF32 = -1.0 * binLoF32; } - if (!isfinite(binL2F32) && isfinite(binH2F32)) { binL2F32 = -1.0 * binH2F32; } - if (!isfinite(binH2F32) && isfinite(binL2F32)) { binH2F32 = -1.0 * binL2F32; } - if (!isfinite(binL4F32) && isfinite(binH4F32)) { binL4F32 = -1.0 * binH4F32; } - if (!isfinite(binH4F32) && isfinite(binL4F32)) { binH4F32 = -1.0 * binL4F32; } - // ADD step 5: Determine SIGMA as the distance between binL2 and binH2 (+/- 0.5 sigma) - - float sigma1 = (binH2F32 - binL2F32); float sigma2 = (binHiF32 - binLoF32) / 2.0; float sigma4 = (binH4F32 - binL4F32) / 4.0; - // Fix again? - if (!isfinite(sigma1) && isfinite(sigma2) && isfinite(sigma4)) { sigma1 = (sigma2 + sigma4) / 2.0; } - if (!isfinite(sigma2) && isfinite(sigma1) && isfinite(sigma4)) { sigma2 = (sigma1 + sigma4) / 2.0; } - if (!isfinite(sigma4) && isfinite(sigma2) && isfinite(sigma1)) { sigma4 = (sigma2 + sigma1) / 2.0; } - // take the smallest of the three: if we have a clump with wide outliers, sigma2 and // sigma4 will be biased high; if we have a bi-modal distribution, sigma1 and sigma2 // will be biased high. - // sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4)); - // CZW: Instead, take the median. Taking the MIN forces a bias on unbiased data. - // It seems like occasionally getting the wrong answer on a complex distribution - // is more acceptable than always getting the wrong answer for simple ones. - - - sigma = PS_MAX( PS_MIN(sigma1,sigma2), - PS_MIN( PS_MAX(sigma1,sigma2), - sigma4)); + sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4)); psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1); @@ -1030,48 +977,5 @@ psTrace(TRACE, 6, "The current sigma is %f.\n", sigma); - // stats->robustStdev = sigma; - stats->robustStdev = sigma; - -#if (CZW && 0) - // Skewness check: Find least biased sample for each pair. - sigma1 = 2.0 * PS_MIN(binH2F32 - stats->robustMedian, - stats->robustMedian - binL2F32); - sigma2 = 1.0 * PS_MIN(binHiF32 - stats->robustMedian, - stats->robustMedian - binLoF32); - sigma4 = 0.5 * PS_MIN(binH4F32 - stats->robustMedian, - stats->robustMedian - binL4F32); - // Kurtosis check: Take median sample as the solution. - stats->robustStdev = PS_MAX( PS_MIN(sigma1,sigma2), - PS_MIN( PS_MAX(sigma1,sigma2), - sigma4)); -#endif - - -#if (CZW) - printf("CZW: bad sigma?: %f %f %f %f %f %f %f %f %f %f\n", - binH2F32,binL2F32,binHiF32,binLoF32,binH4F32,binL4F32, - sigma1,sigma2,sigma4,sigma); - - printf("CZW (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n", - iterate, - stats->robustMedian,stats->robustStdev, - fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]), - - cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2], - cumulative->bounds->data.F32[binMedian-1], - cumulative->bounds->data.F32[binMedian], - cumulative->bounds->data.F32[binMedian+1], - cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3], - - cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2], - cumulative->nums->data.F32[binMedian-1], - cumulative->nums->data.F32[binMedian], - cumulative->nums->data.F32[binMedian+1], - cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]); - // PS_VECTOR_PRINT_F32(histogram->bounds); - // PS_VECTOR_PRINT_F32(histogram->nums); - // PS_VECTOR_PRINT_F32(cumulative->bounds); - // PS_VECTOR_PRINT_F32(cumulative->nums); -#endif + stats->robustStdev = sigma; // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more @@ -1079,31 +983,20 @@ if (sigma < (3.0 * binSize)) { psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); - - // these limits are supposed to be 25 x the raw bin size, NOT 25 of the cumulative histogram bins - psF32 medianLo = stats->robustMedian - 25*binSize; - psF32 medianHi = stats->robustMedian + 25*binSize; - - // long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region - // long maskHi = PS_MIN(cumulative->bounds->n - 1, (binMedian + 25)); // High index for masking - // psF32 medianLo = cumulative->bounds->data.F32[maskLo]; // Value at low index - // psF32 medianHi = cumulative->bounds->data.F32[maskHi]; // Value at high index + long maskLo = PS_MAX(0, (binMedian - 25)); // Low index for masking region + long maskHi = PS_MIN(histogram->bounds->n - 1, (binMedian + 25)); // High index for masking + psF32 medianLo = histogram->bounds->data.F32[maskLo]; // Value at low index + psF32 medianHi = histogram->bounds->data.F32[maskHi]; // Value at high index psTrace(TRACE, 6, "Masking data more than 25 bins from the median\n"); - // psTrace(TRACE, 6, "The median is at bin number %ld. We mask bins outside the bin range (%ld:%ld)\n", binMedian, maskLo, maskHi); + psTrace(TRACE, 6, + "The median is at bin number %ld. We mask bins outside the bin range (%ld:%ld)\n", + binMedian, maskLo, maskHi); psTrace(TRACE, 6, "Masking data outside (%f %f)\n", medianLo, medianHi); - int Nmasked = 0; for (long i = 0 ; i < myVector->n ; i++) { if ((myVector->data.F32[i] < medianLo) || (myVector->data.F32[i] > medianHi)) { - if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue; + // XXXX is this correct? is MASK_MARK safe? mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK; psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]); - Nmasked ++; } } - - if (Nmasked == 0) { - // no significant change to the sigma & binsize -- we are done here - iterate = -1; - continue; - } // Free the histograms; they will be recreated on the next iteration, with new bounds @@ -1137,5 +1030,5 @@ } } - + // XXX test lines while studying algorithm errors // fprintf (stderr, "robust stats test %7.1f +/- %7.1f : %4ld %4ld %4ld %4ld %4ld : %f %f %f\n", @@ -1147,12 +1040,12 @@ PS_BIN_FOR_VALUE (binLo25, cumulative->nums, totalDataPoints * 0.25f, 0); PS_BIN_FOR_VALUE (binHi25, cumulative->nums, totalDataPoints * 0.75f, 0); - psTrace(TRACE, 6, "The 25-percent and 75-percent data point bins are (%ld, %ld).\n", binLo25, binHi25); + psTrace(TRACE, 6, "The 25-percent and 75-precent data point bins are (%ld, %ld).\n", binLo25, binHi25); // ADD step 8: Interpolate to find these two positions exactly: these are the upper and lower quartile // positions. - psF32 binLo25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f); - psF32 binHi25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f); + psF32 binLo25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f); + psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f); if (isnan(binLo25F32) || isnan(binHi25F32)) { - COUNT_WARNING(10, 100, "could not determine the robustUQ or LQ: fitLinearSearchForYThenReturnBin() returned a NAN.\n"); + COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n"); goto escape; } @@ -1207,5 +1100,5 @@ * "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 1D polynomial to ln(y) this + * 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 ********************/ @@ -1242,5 +1135,5 @@ return true; } - if (myVector->n < 1) { printf("There are no elements in this vector.\n"); abort(); } + float guessStdev = stats->robustStdev; // pass the guess sigma float guessMean = stats->robustMedian; // pass the guess mean @@ -1262,6 +1155,6 @@ // 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; + psF32 dN = PS_MAX (1, PS_MIN (4, stats->robustN50 / 8)); + binSize = guessStdev / dN; } @@ -1289,21 +1182,10 @@ // XXX can we calculate the binMin, binMax **before** building this histogram? // if the range is too absurd, adjust numBins & binSize - // We no longer want to reset the binSize here, as it can cause odd things. Better to select - // a number of bins, and then set the min/max values to put those bins sanely around the mean. long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize)); - // binSize = (max - min) / (float) numBins; + binSize = (max - min) / (float) numBins; psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max); psTrace(TRACE, 6, "The new bin size is %f.\n", binSize); psTrace(TRACE, 6, "The numBins is %ld\n", numBins); - -#define FITTED_CLIPPING_NUM 5.0 - if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) { - min = guessMean - FITTED_CLIPPING_NUM * guessStdev; - } - if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) { - max = guessMean + FITTED_CLIPPING_NUM * guessStdev; - } - psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers) if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) { @@ -1340,5 +1222,4 @@ 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"); @@ -1367,5 +1248,4 @@ psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); - float lowfitMean = NAN; float lowfitStdev = NAN; @@ -1405,5 +1285,5 @@ } y->n = x->n = j; - + // fit 2nd order polynomial to ln(y) = -(x-xo)^2/2sigma^2 // XXX this fit may fail with an error for an ill-conditioned matrix (bad data) @@ -1417,5 +1297,4 @@ psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); - psFree(poly); psFree(histogram); @@ -1499,12 +1378,4 @@ psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x); -#if (CZW && 0) - for (long i = 0; i < x->n; i++) { - printf("CZW: Dcheck: %ld %f %f %f\n", - i,x->data.F32[i],y->data.F32[i], - poly->coeff[0] + poly->coeff[1] * x->data.F32[i] + - poly->coeff[2] * pow(x->data.F32[i],2)); - } -#endif psFree(x); psFree(y); @@ -1522,5 +1393,4 @@ 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]); @@ -1545,5 +1415,4 @@ } - psFree (poly); } @@ -1568,18 +1437,4 @@ done = true; } - - -#if (CZW && 1) - printf("CZW IN FITTED: iter %d %f \n" - " low %f %f \n" - " full %f %f \n" - " robust %f %f \n" - " final %f %f\n", - iteration,minValueSym, - lowfitMean,lowfitStdev, - fullfitMean,fullfitStdev, - stats->robustMedian,stats->robustStdev, - guessMean,guessStdev); -#endif // Clean up after fitting @@ -2110,5 +1965,4 @@ // other private functions used above -# if (0) static psF32 QuadraticInverse(psF32 a, psF32 b, @@ -2132,20 +1986,4 @@ return 0.5 * (xLo + xHi); } - -static psF32 LinearInverse(psF32 a, - psF32 b, - psF32 y, - psF32 xLo, - psF32 xHi - ) -{ - psF64 x = (y - b) / a; - - if (xLo <= x && x <= xHi) { - return x; - } - return 0.5 * (xLo + xHi); -} -# endif # if (0) @@ -2404,4 +2242,5 @@ return tmpFloat; } +# endif /****************************************************************************** @@ -2437,34 +2276,23 @@ PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN); - // psVector *x = psVectorAlloc(3, PS_TYPE_F64); - // psVector *y = psVectorAlloc(3, PS_TYPE_F64); - psVector *x = psVectorAlloc(5, PS_TYPE_F64); - psVector *y = psVectorAlloc(5, PS_TYPE_F64); + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); psF32 tmpFloat = 0.0f; - // if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) { - if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) { + if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) { // The general case. We have all three points. - // x->data.F64[0] = binNum - 1; - // x->data.F64[1] = binNum; - // x->data.F64[2] = binNum + 1; - x->data.F64[0] = xVec->data.F32[binNum - 2]; - x->data.F64[1] = xVec->data.F32[binNum - 1]; - x->data.F64[2] = xVec->data.F32[binNum + 0]; - x->data.F64[3] = xVec->data.F32[binNum + 1]; - x->data.F64[4] = xVec->data.F32[binNum + 2]; - y->data.F64[0] = yVec->data.F32[binNum - 2]; - y->data.F64[1] = yVec->data.F32[binNum - 1]; - y->data.F64[2] = yVec->data.F32[binNum + 0]; - y->data.F64[3] = yVec->data.F32[binNum + 1]; - y->data.F64[4] = yVec->data.F32[binNum + 2]; + x->data.F64[0] = binNum - 1; + x->data.F64[1] = binNum; + x->data.F64[2] = binNum + 1; + y->data.F64[0] = yVec->data.F32[binNum - 1]; + y->data.F64[1] = yVec->data.F32[binNum]; + y->data.F64[2] = yVec->data.F32[binNum + 1]; psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]); psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); - // Ensure that the y value lies within range of the y values. - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || - ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) { + if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) || + ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2503,5 +2331,5 @@ psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); - float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]); + float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[2]); psFree(myPoly); @@ -2513,15 +2341,13 @@ return(NAN); } - + // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1 - // assert (binValue >= binNum - 1); - // assert (binValue <= binNum + 1); - - // int fitBin = binValue; - // float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; - // float dY = binValue - fitBin; - // tmpFloat = xVec->data.F32[fitBin] + dY * dX; - tmpFloat = binValue; - + assert (binValue >= binNum - 1); + assert (binValue <= binNum + 1); + + int fitBin = binValue; + float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; + float dY = binValue - fitBin; + tmpFloat = xVec->data.F32[fitBin] + dY * dX; } else { // These are special cases where the bin is at the beginning or end of the vector. @@ -2555,187 +2381,2 @@ return tmpFloat; } -# endif - - -/****************************************************************************** -fitQuadraticSearchForYThenReturnXusingValues(*xVec, *yVec, binNum, yVal): A general routine -which fits a quadratic to three points and returns the x bin value corresponding to the input -y-value. This routine takes psVectors of x/y pairs as input, and fits a quadratic to the 3 -points surrounding element binNum in the vectors. This version uses the values of x[i] for the -x coordinates (not the midpoints). This is appropriate for a cumulative histogram. It then -determines for what value x does that quadratic f(x) = yVal (the input parameter). - -XXX this function is used a fair amount in an inner loop: the polynomial fitting and evaluation -could easily be done with statically allocated doubles, skipping the psLib versions of -polynomial fitting, etc. - -*****************************************************************************/ -static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec, - psVector *yVec, - psS32 binNum, - psF32 yVal - ) -{ - -# if (1) -# define HALF_SIZE 2 - double Sx = 0.0; - - double Sy = 0.0; - double Sxx = 0.0; - double Sxy = 0.0; - double deltaY = 0.0; - int N = 0; - - for (int u = binNum - HALF_SIZE; u <= binNum + HALF_SIZE; u++) { - if ((u >= 0)&&(u < yVec->n)) { - if (u+1 < xVec->n) { - Sx += yVec->data.F32[u]; - Sxx += PS_SQR(yVec->data.F32[u]); - - deltaY = xVec->data.F32[u]; - //deltaY = 0.5 * (xVec->data.F32[u] + xVec->data.F32[u+1]); - Sy += deltaY; - Sxy += yVec->data.F32[u] * deltaY; - N += 1; - } - } - } - double Det = N * Sxx - Sx * Sx; - if (Det == 0.0) return NAN; - if (N == 0) return NAN; - - double C0 = (Sy*Sxx - Sx*Sxy) / Det; - double C1 = (Sxy*N - Sx*Sy) / Det; - - double value = C0 + yVal*C1; - return value; - - -# else - psTrace(TRACE, 5, "binNum, yVal is (%d, %f)\n", binNum, yVal); - if (psTraceGetLevel("psLib.math") >= 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); - PS_ASSERT_VECTOR_TYPE(xVec, PS_TYPE_F32, NAN); - PS_ASSERT_VECTOR_TYPE(yVec, PS_TYPE_F32, NAN); - PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(xVec->n - 1), NAN); - PS_ASSERT_INT_WITHIN_RANGE(binNum, 0, (int)(yVec->n - 1), NAN); - - // psVector *x = psVectorAlloc(3, PS_TYPE_F64); - // psVector *y = psVectorAlloc(3, PS_TYPE_F64); - psVector *x = psVectorAlloc(5, PS_TYPE_F64); - psVector *y = psVectorAlloc(5, PS_TYPE_F64); - psF32 tmpFloat = 0.0f; - - if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) { - x->data.F64[0] = xVec->data.F32[binNum - 2]; - x->data.F64[1] = xVec->data.F32[binNum - 1]; - x->data.F64[2] = xVec->data.F32[binNum + 0]; - x->data.F64[3] = xVec->data.F32[binNum + 1]; - x->data.F64[4] = xVec->data.F32[binNum + 2]; - - y->data.F64[0] = yVec->data.F32[binNum - 2]; - y->data.F64[1] = yVec->data.F32[binNum - 1]; - y->data.F64[2] = yVec->data.F32[binNum + 0]; - y->data.F64[3] = yVec->data.F32[binNum + 1]; - y->data.F64[4] = yVec->data.F32[binNum + 2]; - psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]); - psTrace(TRACE, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]); - psTrace(TRACE, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]); - - // Ensure that the y value lies within range of the y values. - if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || - ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) { - psError(PS_ERR_BAD_PARAMETER_VALUE, true, - _("Specified yVal, %g, is not within y-range, %g to %g."), - (psF64)yVal, y->data.F64[0], y->data.F64[2]); - return NAN; - } - - // Ensure that the y values are monotonic. - if (((y->data.F64[0] < y->data.F64[1]) && !(y->data.F64[1] <= y->data.F64[2])) || - ((y->data.F64[0] > y->data.F64[1]) && !(y->data.F64[1] >= y->data.F64[2]))) { - psError(PS_ERR_UNKNOWN, true, - "This routine must be called with monotonically increasing or decreasing data points.\n"); - psFree(x); - psFree(y); - return NAN; - } - - // Determine the coefficients of the polynomial. - psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 1); - if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { - psError(PS_ERR_UNEXPECTED_NULL, false, - _("Failed to fit a 1-dimensional polynomial to the three specified data points. " - "Returning NAN.")); - psFree(x); - psFree(y); - return NAN; - } - - psTrace(TRACE, 6, "myPoly->coeff[0] is %f\n", myPoly->coeff[0]); - psTrace(TRACE, 6, "myPoly->coeff[1] is %f\n", myPoly->coeff[1]); - psTrace(TRACE, 6, "Fitted y vec is (%f %f)\n", - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[0]), - (psF32) psPolynomial1DEval(myPoly, (psF64) x->data.F64[1])); - - psTrace(TRACE, 6, "We fit the polynomial, now find x such that f(x) equals %f\n", yVal); - float binValue = LinearInverse(myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]); - psFree(myPoly); - - if (isnan(binValue)) { - psError(PS_ERR_UNEXPECTED_NULL, - false, _("Failed to determine the median of the fitted polynomial. Returning NAN.")); - psFree(x); - psFree(y); - return(NAN); - } - - // I believe that mathematically the fitted bin position must be between binNum - 1 and binNum + 1 - // assert (binValue >= binNum - 1); - // assert (binValue <= binNum + 1); - - // int fitBin = binValue; - // float dX = xVec->data.F32[fitBin+1] - xVec->data.F32[fitBin]; - // float dY = binValue - fitBin; - // tmpFloat = xVec->data.F32[fitBin] + dY * dX; - tmpFloat = binValue; - - - } else { - // These are special cases where the bin is at the beginning or end of the vector. - if (binNum == 0) { - // We have two points only at the beginning of the vectors x and y. - // X = (dX/dY)(Y - Yo) + Xo - float dX = xVec->data.F32[1] - xVec->data.F32[0]; - float dY = yVec->data.F32[1] - yVec->data.F32[0]; - if (dY == 0.0) { - tmpFloat = xVec->data.F32[0]; - } else { - tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0]; - } - } else if (binNum == (xVec->n - 1)) { - // We have two points only at the end of the vectors x and y. - // X = (dX/dY)(Y - Yo) + Xo - float dX = xVec->data.F32[binNum] - xVec->data.F32[binNum-1]; - float dY = yVec->data.F32[binNum] - yVec->data.F32[binNum-1]; - if (dY == 0.0) { - tmpFloat = xVec->data.F32[binNum-1]; - } else { - tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1]; - } - } - } - - psTrace(TRACE, 6, "FIT: return %f\n", tmpFloat); - psFree(x); - psFree(y); - - return tmpFloat; -# endif -}