Changeset 36345

Timestamp:

Dec 4, 2013, 2:11:09 PM (13 years ago)

Author:

bills

Message:

Integrate updated psStats code from the trunk

Location:

branches/eam_branches/ipp-20130904/psLib/src/math

Files:

• psMinimizePolyFit.c (modified) (4 diffs, 1 prop)
• psStats.c (modified) (37 diffs, 1 prop)

branches/eam_branches/ipp-20130904/psLib/src/math/psMinimizePolyFit.c

@@ -715 +715 @@
         sigma  = (sorted->data.F64[3 * sorted->n / 4] - sorted->data.F64[sorted->n / 4]); 
         psFree(sorted);
+
+        if ((!isfinite(median))||
+            (!isfinite(sigma))) {
+          median = 0.0;
+          sigma = 1.0;
+        }
+        if ((median == 0.0)&&(sigma == 1.0)) {
+          scale = false;
+        }
         // 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);
@@ -726 +735 @@
         
         result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, z64);
+        psFree(z64); // Done with this.
       }
       else {
@@ -736 +746 @@
       
       if (scale) {
-        psFree(z64); // Done with this.
-
         // Undo scaling in the polynomial values.
         psF64 *Zcoeff = psAlloc((1 + poly->nX) * sizeof(psF64));
@@ -757 +765 @@
 #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),
+                   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)
                    );

branches/eam_branches/ipp-20130904/psLib/src/math/psStats.c

@@ -140 +140 @@
 }
 
+// Debug information
+#define CZW 0
 
 /*****************************************************************************/
@@ -172 +174 @@
 *****************************************************************************/
 
-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 +232 @@
         }
         count++;
+
     }
     if (errors) {
@@ -793 +797 @@
         } 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 +880 @@
         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 +913 @@
 
         // 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 +933 @@
         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 +948 @@
             goto escape;
         }
-
+    
         // ADD step 4b: Interpolate Sigma (linearly) to find these two positions exactly: these are the 1sigma
         // positions.
@@ -938 +960 @@
         // (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 +970 @@
                             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 +989 @@
                 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 +1030 @@
 
         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 Robust (%d): median %f sigma %f delta: %f \n\t %f %f %f %f %f %f %f \n\t %f %f %f %f %f %f %f\n",
+               iterate,
+               stats->robustMedian,stats->robustStdev,
+               fabs(cumulative->bounds->data.F32[binMedian] - cumulative->bounds->data.F32[binMedian + 1]),
+               
+               cumulative->bounds->data.F32[binMedian-3],cumulative->bounds->data.F32[binMedian-2],
+               cumulative->bounds->data.F32[binMedian-1],
+               cumulative->bounds->data.F32[binMedian],
+               cumulative->bounds->data.F32[binMedian+1],
+               cumulative->bounds->data.F32[binMedian+2],cumulative->bounds->data.F32[binMedian+3],
+               
+               cumulative->nums->data.F32[binMedian-3],cumulative->nums->data.F32[binMedian-2],
+               cumulative->nums->data.F32[binMedian-1],
+               cumulative->nums->data.F32[binMedian],
+               cumulative->nums->data.F32[binMedian+1],
+               cumulative->nums->data.F32[binMedian+2],cumulative->nums->data.F32[binMedian+3]);
+        //      PS_VECTOR_PRINT_F32(histogram->bounds);
+        //      PS_VECTOR_PRINT_F32(histogram->nums);
+        //      PS_VECTOR_PRINT_F32(cumulative->bounds);
+        //      PS_VECTOR_PRINT_F32(cumulative->nums);
+#endif
 
         // ADD step 6: If the measured SIGMA is less than 2 times the bin size, exclude points which are more
@@ -983 +1079 @@
         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 +1137 @@
         }
     }
-
+    
     // 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 +1147 @@
     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 +1207 @@
  * "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 +1242 @@
         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 +1262 @@
             // 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));
+
+          // CZW 2013-11-20: We know that the histogram is going to be basically Gaussian.
+          // Furthermore, we only use the inner +/- 2 sigma parts.  Therefore, define the
+          // binsize such that the bin at 2 sigma contains ~50 points (S/N ~ 7).  robustN50
+          // contains half the total points, so 2 * robustN50 / 50 is the fraction of all
+          // points in the 2 sigma bin.  Dance the erf() relations around, and it looks like
+          // there's a factor of about 1/20 to include.  Keep the PS_MAX to ensure we never bin
+          // wider than 1 sigma when the number of points is small.
+          psF32 dN = PS_MAX(1, (stats->robustN50 / 500.0));
+          binSize = guessStdev / dN;
         }
 
@@ -1182 +1298 @@
         // XXX can we calculate the binMin, binMax **before** building this histogram?
         // if the range is too absurd, adjust numBins & binSize
-        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
-        binSize = (max - min) / (float) numBins;
+        // We no longer want to reset the binSize here, as it can cause odd things.  Better to select
+        // a number of bins, and then set the min/max values to put those bins sanely around the mean.
+        //        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
+        //        binSize = (max - min) / (float) numBins;
         psTrace(TRACE, 6, "The new min/max values are (%f, %f).\n", min, max);
         psTrace(TRACE, 6, "The new bin size is %f.\n", binSize);
-        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
-
+        //        psTrace(TRACE, 6, "The numBins is %ld\n", numBins);
+
+
+#define FITTED_CLIPPING_NUM 5.0
+        if (min < guessMean - FITTED_CLIPPING_NUM * guessStdev) {
+          min = guessMean - FITTED_CLIPPING_NUM * guessStdev;
+        }
+        if (max > guessMean + FITTED_CLIPPING_NUM * guessStdev) {
+          max = guessMean + FITTED_CLIPPING_NUM * guessStdev;
+        }
+        long numBins = PS_MAX (50, PS_MIN (10000, (max - min) / binSize));
+        if (CZW) { printf("I've clipped: %f %f => %f %f ; %f %ld\n",guessMean,guessStdev,min,max,binSize,stats->robustN50); }
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
@@ -1222 +1350 @@
         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");
@@ -1238 +1367 @@
             }
             psTrace (TRACE, 6, "(%f = %.0f) ", histogram->bounds->data.F32[i], histogram->nums->data.F32[i]);
+            if (CZW) { printf("CENTERED_HISTOGRAM: %f %f\n",
+                              PS_BIN_MIDPOINT(histogram,i),
+                              histogram->nums->data.F32[i]); }
         }
         psTrace (TRACE, 6, "\n");
 
+        if (CZW) { printf("Bin selection done: %ld %f %f %ld %f %f %ld %f %f\n",
+                          binMin,PS_BIN_MIDPOINT(histogram,binMin),histogram->nums->data.F32[binMin],
+                          binMax,PS_BIN_MIDPOINT(histogram,binMax),histogram->nums->data.F32[binMax],
+                          binPeak,PS_BIN_MIDPOINT(histogram,binPeak),histogram->nums->data.F32[binPeak]);
+        }
+        
         // assume a reasonably well-defined gaussian-like population; run from peak out until val < 0.25*peak
         psTrace(TRACE, 6, "The clipped numBins is %ld\n", binMax - binMin);
@@ -1248 +1386 @@
         psTrace(TRACE, 6, "The clipped peak value is %f\n", histogram->nums->data.F32[binPeak]);
 
+        
         float lowfitMean = NAN;
         float lowfitStdev = NAN;
@@ -1285 +1424 @@
             }
             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)
@@ -1291 +1430 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+#if (CZW && 1)
+            printf("CZW: LowfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]);
+            for (long i = 0; i < x->n; i++) {
+              printf("CZW: Lowfit: %d %ld %f %f %f\n",
+                     status,i,x->data.F32[i],y->data.F32[i],
+                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+                     poly->coeff[2] * pow(x->data.F32[i],2));
+            }
+#endif
             psFree(x);
             psFree(y);
@@ -1305 +1453 @@
             if (poly->coeff[2] >= 0.0) {
                 COUNT_WARNING(10, 100, "Failed parabolic fit: %f + %f x + %f x^2\n", poly->coeff[0], poly->coeff[1], poly->coeff[2]);
+
                 psFree(poly);
                 psFree(histogram);
@@ -1357 +1506 @@
                 }
             }
+
             psTrace(TRACE, 6, "Lower bound for symmetric range: %f (%ld)\n",
                     PS_BIN_MIDPOINT(histogram, binS), binS);
@@ -1378 +1528 @@
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
+#if (CZW && 1)
+            printf("CZW: FullfitPoly: %f %f %f\n",poly->coeff[0],poly->coeff[1],poly->coeff[2]);
+            for (long i = 0; i < x->n; i++) {
+              printf("CZW: Fullfit: %d %ld %f %f %f\n",
+                     status,i,x->data.F32[i],y->data.F32[i],
+                     poly->coeff[0] + poly->coeff[1] * x->data.F32[i] +
+                     poly->coeff[2] * pow(x->data.F32[i],2));
+            }
+#endif
             psFree(x);
             psFree(y);
@@ -1393 +1552 @@
             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 +1575 @@
             }
 
+            
             psFree (poly);
         }
@@ -1437 +1598 @@
             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 +2140 @@
 // other private functions used above
 
+# if (0)
 static psF32 QuadraticInverse(psF32 a,
                               psF32 b,
@@ -1986 +2162 @@
     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 +2434 @@
     return tmpFloat;
 }
-# endif
 
 /******************************************************************************
@@ -2276 +2467 @@
     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 +2533 @@
 
         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 +2543 @@
             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 +2585 @@
     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
+}