Index: /branches/eam_branches/ipp-20130904/psLib/src/math/Makefile.am =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/math/Makefile.am (revision 36255) +++ /branches/eam_branches/ipp-20130904/psLib/src/math/Makefile.am (revision 36256) @@ -7,4 +7,5 @@ psUnaryOp.c \ psBinaryOp.c \ + psBinomialCoeff.c \ psClip.c \ psCompare.c \ @@ -36,4 +37,5 @@ psUnaryOp.h \ psBinaryOp.h \ + psBinomialCoeff.h \ psClip.h \ psCompare.h \ Index: /branches/eam_branches/ipp-20130904/psLib/src/math/psBinomialCoeff.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/math/psBinomialCoeff.c (revision 36256) +++ /branches/eam_branches/ipp-20130904/psLib/src/math/psBinomialCoeff.c (revision 36256) @@ -0,0 +1,30 @@ +/** @file psBinomialCoeff.c + * + * @brief Calculate binomial coefficient + * + * @author Chris Waters, IfA + * + */ + +#ifdef HAVE_CONFIG_H +# include "config.h" +#endif + +#include "psError.h" + +#include "psBinomialCoeff.h" + +psS64 psBinomialCoeff( psS64 N, + psS64 k) { + psS64 i; + psS64 b = 1; + if (N == k) { return (b); } + if (k == 0) { return (b); } + + for (i = 1; i <= k; i++) { + b *= (N - (k - i))/i; + } + + return(b); +} + Index: /branches/eam_branches/ipp-20130904/psLib/src/math/psBinomialCoeff.h =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/math/psBinomialCoeff.h (revision 36256) +++ /branches/eam_branches/ipp-20130904/psLib/src/math/psBinomialCoeff.h (revision 36256) @@ -0,0 +1,8 @@ +#ifndef PS_BINCOEFF_H +#define PS_BINCOEFF_H + +// Calculate the binomial coefficient +psS64 psBinomialCoeff( psS64 N, + psS64 k); + +#endif Index: /branches/eam_branches/ipp-20130904/psLib/src/math/psMinimizePolyFit.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/math/psMinimizePolyFit.c (revision 36255) +++ /branches/eam_branches/ipp-20130904/psLib/src/math/psMinimizePolyFit.c (revision 36256) @@ -42,7 +42,9 @@ #include "psLogMsg.h" #include "psMathUtils.h" +#include "psBinomialCoeff.h" /*****************************************************************************/ /* DEFINE STATEMENTS */ /*****************************************************************************/ +#define CZW 0 # define USE_GAUSS_JORDAN 1 @@ -603,5 +605,26 @@ bool status = false; if (USE_GAUSS_JORDAN) { + +#if (CZW) + printf("CZW: about to do GJ: %d\n",status); + for (psS32 k = 0; k < nTerm; k++) { + printf("CZW: %d %f \t",k,B->data.F64[k]); + for (psS32 kk = 0; kk < nTerm; kk++) { + printf(" %f ",(A->data.F64[k][kk])); + } + printf("\n"); + } +#endif status = psMatrixGJSolve(A, B); +#if (CZW) + printf("CZW: just did GJ: %d\n",status); + for (psS32 k = 0; k < nTerm; k++) { + printf("CZW: %d %f \t",k,B->data.F64[k]); + for (psS32 kk = 0; kk < nTerm; kk++) { + printf(" %f ",(A->data.F64[k][kk])); + } + printf("\n"); + } +#endif } else { status = psMatrixLUSolve(A, B); @@ -676,10 +699,80 @@ bool result = true; + // Define values that may be used by PS_POLYNOMIAL_ORD form. + bool scale = false; + double median = 0.0; + double sigma = 1.0; + psVector *sorted = NULL; + psVector *z64 = NULL; + switch (poly->type) { case PS_POLYNOMIAL_ORD: - result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); - if (!result) { - psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); - } + if ((f64->n < 10000)&&(poly->nX > 1)) { + scale = true; + sorted = psVectorSort(NULL,x64); + median = sorted->data.F64[sorted->n / 2]; + // CZW: I'm not bothering to scale this because it doesn't really matter. + sigma = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]); + psFree(sorted); + // I can't see a way to not clobber x if it's already F64, so make a copy.x + z64 = psVectorCopy(NULL,x64,PS_TYPE_F64); + psBinaryOp(z64,z64,"-",psScalarAlloc(median,PS_TYPE_F64)); + psBinaryOp(z64,z64,"/",psScalarAlloc(sigma,PS_TYPE_F64)); + +#if (CZW) + printf("poly1d: Scale parameters: %f %f\n",median,sigma); +#endif + + + result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64); + } + else { + result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); + } + + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); + } + + if (scale) { + psFree(z64); // Done with this. + + // Undo scaling in the polynomial values. + psF64 *Zcoeff = psAlloc((1 + poly->nX) * sizeof(psF64 *)); + psF64 *ZcoeffErr = psAlloc((1 + poly->nX) * sizeof(psF64 *)); + + for (psS32 i = 0; i <= poly->nX; i++) { + Zcoeff[i] = poly->coeff[i]; + ZcoeffErr[i] = poly->coeffErr[i]; +#if (CZW) + printf("poly1d: fit parameters: %d %f %f\n", + i,poly->coeff[i],poly->coeffErr[i]); +#endif + } + for (psS32 i = 0; i <= poly->nX; i++) { + poly->coeff[i] = 0.0; + poly->coeffErr[i] = 0.0; + + for (psS32 j = 0; j <= poly->nX; j++) { +#if (CZW) + printf(" %d %d %f %f %f %f => %f\n", + i,j,Zcoeff[j],pow(1.0 / sigma,j) * pow(-1,j - i),pow(median,j - i),1.0 * psBinomialCoeff(j,i), + Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j -i) * pow(median,j - i) * 1.0 * psBinomialCoeff(j,i) + ); +#endif + poly->coeff[i] += Zcoeff[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - i) * psBinomialCoeff(j,i); + poly->coeffErr[i] += pow(ZcoeffErr[j] * pow(1.0 / sigma,j) * pow(-1,j - i) * pow(median,j - 1) * psBinomialCoeff(j,i),2); + } + poly->coeffErr[i] = sqrt(poly->coeffErr[i]); +#if (CZW) + printf("poly1d: unscaled parameters: %d %f %f\n", + i,poly->coeff[i], poly->coeffErr[i]); +#endif + } + psFree(Zcoeff); + psFree(ZcoeffErr); + + } + break; case PS_POLYNOMIAL_CHEB: Index: /branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c =================================================================== --- /branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c (revision 36255) +++ /branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c (revision 36256) @@ -140,4 +140,6 @@ } +// Debug information +#define CZW 0 /*****************************************************************************/ @@ -172,5 +174,6 @@ *****************************************************************************/ -static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); +// static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); +static psF32 fitLinearSearchForYThenReturnBin(const psVector *xVec, psVector *yVec, psS32 binNum, psF32 yVal); /****************************************************************************** @@ -229,4 +232,5 @@ } count++; + } if (errors) { @@ -793,5 +797,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); @@ -876,8 +880,22 @@ cumulative = psHistogramAlloc(min, max, numBins); cumulative->nums->data.F32[0] = histogram->nums->data.F32[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]; - } + cumulative->bounds->data.F32[0] = histogram->bounds->data.F32[1]; + + // Correctly fill the cumulative distribution with monotonically increasing values (skip zero valued bins). + long Nc = 1; // track the current bin of cumulative + // the boundaries for the current cumulative bin are from upper end of the last valid histogram bin to the + // upper end of the current histogram bin + for (long i = 1; i < histogram->nums->n - 1; i++) { + if (histogram->nums->data.F32[i] == 0.0) continue; + cumulative->nums->data.F32[Nc] = cumulative->nums->data.F32[Nc - 1] + histogram->nums->data.F32[i]; + cumulative->bounds->data.F32[Nc] = histogram->bounds->data.F32[i+1]; + Nc ++; + } + long Nlast = Nc - 1; // last valid cumulative bin + for (long i = Nc; i < histogram->nums->n; i++) { // Ensure the unused entries are filled. + cumulative->nums->data.F32[i] = cumulative->nums->data.F32[Nlast]; + cumulative->bounds->data.F32[i] = cumulative->bounds->data.F32[i-1] + 1.0; + } + if (psTraceGetLevel("psLib.math") >= 8) { PS_VECTOR_PRINT_F32(cumulative->bounds); @@ -895,7 +913,10 @@ // 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. @@ -912,7 +933,8 @@ 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.022481f, 0); - PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977519f, 0); - + PS_BIN_FOR_VALUE(binL4, cumulative->nums, totalDataPoints * 0.022750f, 0); + PS_BIN_FOR_VALUE(binH4, cumulative->nums, totalDataPoints * 0.977250f, 0); + + psTrace(TRACE, 6, "The 15.8655%% and 84.1345%% data point bins are (%ld, %ld).\n", binLo, binHi); @@ -926,5 +948,5 @@ goto escape; } - + // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma // positions. @@ -938,4 +960,5 @@ // (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); @@ -947,8 +970,15 @@ totalDataPoints * 0.691462f); PS_BIN_INTERPOLATE (binL4F32, cumulative->nums, cumulative->bounds, binL4, - totalDataPoints * 0.022481f); + totalDataPoints * 0.022750f); PS_BIN_INTERPOLATE (binH4F32, cumulative->nums, cumulative->bounds, binH4, - totalDataPoints * 0.977519f); - + 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 // report +/- 1 sigma points psTrace(TRACE, 5, @@ -959,16 +989,39 @@ binL2F32, binH2F32); psTrace(TRACE, 5, - "The exact 02.22481 and 97.7519 percent data point positions are: (%f, %f)\n", + "The exact 02.2275 and 97.7250 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)); + // sigma = PS_MIN (sigma1, PS_MIN (sigma2, sigma4)); + // CZW: Instead, take the median. Taking the MIN forces a bias on unbiased data. + // It seems like occasionally getting the wrong answer on a complex distribution + // is more acceptable than always getting the wrong answer for simple ones. + + + sigma = PS_MAX( PS_MIN(sigma1,sigma2), + PS_MIN( PS_MAX(sigma1,sigma2), + sigma4)); psTrace(TRACE, 6, "The 1x sigma is %f.\n", sigma1); @@ -977,5 +1030,48 @@ psTrace(TRACE, 6, "The current sigma is %f.\n", sigma); - stats->robustStdev = 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 // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more @@ -983,20 +1079,31 @@ if (sigma < (3.0 * binSize)) { psTrace(TRACE, 6, "*************: Do another iteration (%f %f).\n", sigma, binSize); - 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 + + // 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 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)) { - // XXXX is this correct? is MASK_MARK safe? + if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & MASK_MARK) continue; mask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= MASK_MARK; psTrace(TRACE, 6, "Masking element %ld is %f\n", i, myVector->data.F32[i]); + Nmasked ++; } } + + 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 @@ -1030,5 +1137,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", @@ -1040,12 +1147,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-precent data point bins are (%ld, %ld).\n", binLo25, binHi25); + psTrace(TRACE, 6, "The 25-percent and 75-percent 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 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f); - psF32 binHi25F32 = fitQuadraticSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f); + psF32 binLo25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binLo25, totalDataPoints * 0.25f); + psF32 binHi25F32 = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binHi25, totalDataPoints * 0.75f); if (isnan(binLo25F32) || isnan(binHi25F32)) { - COUNT_WARNING(10, 100, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnBin() returned a NAN.\n"); + COUNT_WARNING(10, 100, "could not determine the robustUQ or LQ: fitLinearSearchForYThenReturnBin() returned a NAN.\n"); goto escape; } @@ -1100,5 +1207,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 2D polynomial to ln(y) this + * robustN50 should also be set gaussian fit is performed using 1D polynomial to ln(y) this * version follows the upper portion of the distribution until it passes 0.5*peak ********************/ @@ -1135,5 +1242,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 @@ -1155,6 +1262,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; } @@ -1182,10 +1289,21 @@ // 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)) { @@ -1222,4 +1340,5 @@ 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"); @@ -1248,4 +1367,5 @@ psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]); + float lowfitMean = NAN; float lowfitStdev = NAN; @@ -1285,5 +1405,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) @@ -1297,4 +1417,5 @@ psErrorClear(); COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n"); + psFree(poly); psFree(histogram); @@ -1378,4 +1499,12 @@ 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); @@ -1393,4 +1522,5 @@ 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]); @@ -1415,4 +1545,5 @@ } + psFree (poly); } @@ -1437,4 +1568,18 @@ 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 @@ -1965,4 +2110,5 @@ // other private functions used above +# if (0) static psF32 QuadraticInverse(psF32 a, psF32 b, @@ -1986,4 +2132,20 @@ 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) @@ -2242,5 +2404,4 @@ return tmpFloat; } -# endif /****************************************************************************** @@ -2276,23 +2437,34 @@ 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(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 >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) { + // if ((binNum >= 1) && (binNum <= (yVec->n - 2)) && (binNum <= (xVec->n - 2))) { + if ((binNum >= 2) && (binNum <= (yVec->n - 3)) && (binNum <= (xVec->n - 3))) { // The general case. We have all three points. - x->data.F64[0] = binNum - 1; - x->data.F64[1] = binNum; - x->data.F64[2] = binNum + 1; - 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]; + // x->data.F64[0] = binNum - 1; + // x->data.F64[1] = binNum; + // x->data.F64[2] = binNum + 1; + x->data.F64[0] = xVec->data.F32[binNum - 2]; + x->data.F64[1] = xVec->data.F32[binNum - 1]; + x->data.F64[2] = xVec->data.F32[binNum + 0]; + x->data.F64[3] = xVec->data.F32[binNum + 1]; + x->data.F64[4] = xVec->data.F32[binNum + 2]; + y->data.F64[0] = yVec->data.F32[binNum - 2]; + y->data.F64[1] = yVec->data.F32[binNum - 1]; + y->data.F64[2] = yVec->data.F32[binNum + 0]; + y->data.F64[3] = yVec->data.F32[binNum + 1]; + y->data.F64[4] = yVec->data.F32[binNum + 2]; psTrace(TRACE, 6, "x vec (orig) is (%f %f %f %f)\n", xVec->data.F32[binNum - 1], xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]); 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[2])) || - ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) { + if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[4])) || + ((y->data.F64[4] <= yVal) && (yVal <= y->data.F64[0]))) ) { psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Specified yVal, %g, is not within y-range, %g to %g."), @@ -2331,5 +2503,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[2]); + float binValue = QuadraticInverse(myPoly->coeff[2], myPoly->coeff[1], myPoly->coeff[0], yVal, x->data.F64[0], x->data.F64[4]); psFree(myPoly); @@ -2341,13 +2513,15 @@ 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; + // 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. @@ -2381,2 +2555,187 @@ 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 +}