Changeset 18155 for trunk/psLib/src/math/psStats.c

Timestamp:

Jun 16, 2008, 4:26:34 PM (18 years ago)

Author:

Paul Price

Message:

Robust statistics should fall over gracefully when the iteration limit is reached.

File:

• trunk/psLib/src/math/psStats.c (modified) (11 diffs)

trunk/psLib/src/math/psStats.c

@@ -13 +13 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.223 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-05-09 21:34:42 $
+ *  @version $Revision: 1.224 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-06-17 02:26:34 $
  *
  *  Copyright 2006 IfA, University of Hawaii
@@ -68 +68 @@
 #define PS_POLY_MEDIAN_MAX_ITERATIONS 30
 #define MASK_MARK 0x80   // bit to use internally to mark data as bad
+#define PS_ROBUST_MAX_ITERATIONS 20     // Maximum number of iterations for robust statistics
 
 #define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \
@@ -752 +753 @@
     long binMedian;
 
-    // Iterate to get the best bin size
+    // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop.
     for (int iterate = 1; iterate > 0; iterate++) {
         psTrace(TRACE, 6,
@@ -768 +769 @@
         }
         psTrace(TRACE, 6, "Data min/max is (%.2f, %.2f)\n", min, max);
-        assert (iterate < 20);
-        // fprintf (stderr, "%d @ %d : %f - %f\n", iterate, numValid, min, max);
 
         // If all data points have the same value, then we set the appropriate members of stats and return.
@@ -817 +816 @@
 
         // ADD step 1: Convert the specific histogram to a cumulative histogram
-        // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
+        // The cumulative histogram data points correspond to the UPPER bound value (N < Bin[i+1])
         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[i-1] = histogram->bounds->data.F32[i];
         }
         if (psTraceGetLevel("psLib.math") >= 8) {
@@ -833 +832 @@
         psTrace(TRACE, 6, "Total data points is %ld\n", totalDataPoints);
 
-        // find bin which is the lower bound of median (value[binMedian] < median < value[binMedian+1]
+        // find bin which is the lower bound of median (value[binMedian] < median < value[binMedian+1]
         PS_BIN_FOR_VALUE(binMedian, cumulative->nums, totalDataPoints/2.0, 0);
 
@@ -840 +839 @@
 
         // ADD step 3: Interpolate to the exact 50% position: this is the robust histogram median.
-        stats->robustMedian = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
+        stats->robustMedian = fitQuadraticSearchForYThenReturnXusingValues(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
         if (isnan(stats->robustMedian)) {
             psError(PS_ERR_UNKNOWN, false,
@@ -940 +939 @@
                 }
             }
+
             // Free the histograms; they will be recreated on the next iteration, with new bounds
             psFree(histogram);
@@ -946 +946 @@
             psFree(cumulative);
             cumulative = NULL;
+
+            if (iterate >= PS_ROBUST_MAX_ITERATIONS) {
+                // This occurs when a large number of the values are identical --- a bin size cannot be found
+                // that will spread out the distribution.  Therefore, set what we can, and fall over
+                // gracefully.
+                // stats->robustMedian has already been set
+                stats->robustStdev = sigma;
+                stats->robustUQ = stats->robustMedian;
+                stats->robustLQ = stats->robustMedian;
+                stats->robustN50 = numValid;
+                stats->results |= PS_STAT_ROBUST_MEDIAN;
+                stats->results |= PS_STAT_ROBUST_STDEV;
+                stats->results |= PS_STAT_ROBUST_QUARTILE;
+                psTrace(TRACE, 5, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
+                psTrace(TRACE, 4, "---- %s(0) end  ----\n", __func__);
+                psFree(mask);
+                psFree(statsMinMax);
+                return false;
+            }
         } else {
             // We've got the bin size correct now
@@ -2729 +2748 @@
 *****************************************************************************/
 static psF32 fitQuadraticSearchForYThenReturnXusingValues(const psVector *xVec,
-                                                          psVector *yVec,
-                                                          psS32 binNum,
-                                                          psF32 yVal
+                                                          psVector *yVec,
+                                                          psS32 binNum,
+                                                          psF32 yVal
     )
 {
@@ -2824 +2843 @@
         if (binNum == 0) {
             // We have two points only at the beginning of the vectors x and y.
-            // XXX this does not seem to be doing the linear interpolation / extrapolation
+            // XXX this does not seem to be doing the linear interpolation / extrapolation
             tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                               xVec->data.F32[binNum + 1]);