Changeset 27839 for branches/simmosaic_branches/psLib/src/math/psStats.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• psLib/src/math/psStats.c (modified) (31 diffs)

branches/simmosaic_branches/psLib/src/math/psStats.c

@@ -129 +129 @@
 
 # define COUNT_WARNING(LIMIT, INTERVAL, ...) { \
-        static int nCalls = 1; \
-        if (nCalls < LIMIT) { \
-            psWarning(__VA_ARGS__); \
-        } \
-        if (!(nCalls % INTERVAL)) { \
-            psWarning(__VA_ARGS__); \
+        static int nCalls = 1; \
+        if (nCalls < LIMIT) { \
+            psWarning(__VA_ARGS__); \
+        } \
+        if (!(nCalls % INTERVAL)) { \
+            psWarning(__VA_ARGS__); \
             psWarning("(warning raised %d times)", nCalls); \
-        } \
-        nCalls ++; \
+        } \
+        nCalls ++; \
 }
- 
+
 
 /*****************************************************************************/
@@ -217 +217 @@
     for (long i = 0; i < numData; i++) {
         // Check if the data is with the specified range
-        if (!isfinite(data[i])) 
-            continue;
+        if (!isfinite(data[i]))
+            continue;
         if (useRange && (data[i] < stats->min))
             continue;
@@ -242 +242 @@
 
     if (!isnan(mean)) {
-        stats->results |= PS_STAT_SAMPLE_MEAN;
+        stats->results |= PS_STAT_SAMPLE_MEAN;
     }
     return true;
@@ -278 +278 @@
     for (long i = 0; i < num; i++) {
         // Check if the data is with the specified range
-        if (!isfinite(vector[i])) 
-            continue;
+        if (!isfinite(vector[i]))
+            continue;
         if (useRange && (vector[i] < stats->min))
             continue;
@@ -329 +329 @@
     // into the temporary vectors.
     for (long i = 0; i < inVector->n; i++) {
-        if (!isfinite(input[i])) 
-            continue;
+        if (!isfinite(input[i]))
+            continue;
         if (useRange && (input[i] < stats->min))
             continue;
@@ -353 +353 @@
     // Sort the temporary vector.
     if (!psVectorSort(outVector, outVector)) { // Sort in-place (since it's a copy, it's OK)
-        // an error in psVectorSort is a serious error
+        // an error in psVectorSort is a serious error
         psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to sort input data."));
         stats->sampleUQ = NAN;
@@ -407 +407 @@
     // If the mean is NAN, then generate a warning and set the stdev to NAN.
     if (isnan(stats->sampleMean)) {
-        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
+        COUNT_WARNING(10, 100, "WARNING: vectorSampleStdev(): sample mean is NAN. Setting stats->sampleStdev = NAN.");
         stats->sampleStdev = NAN;
         return true;
@@ -425 +425 @@
     for (long i = 0; i < myVector->n; i++) {
         // Check if the data is with the specified range
-        if (!isfinite(data[i])) 
-            continue;
+        if (!isfinite(data[i]))
+            continue;
         if (useRange && (data[i] < stats->min)) {
             continue;
@@ -502 +502 @@
     for (long i = 0; i < myVector->n; i++) {
         // Check if the data is with the specified range
-        if (!isfinite(data[i])) 
-            continue;
+        if (!isfinite(data[i]))
+            continue;
         if (useRange && (data[i] < stats->min)) {
             continue;
@@ -749 +749 @@
     // 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,
-                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
-                iterate);
+        psTrace(TRACE, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
+
+        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.
+          COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
+          goto escape;
+        }
 
         // Get the minimum and maximum values
@@ -771 +777 @@
             stats->robustLQ = min;
             stats->robustN50 = numValid;
-            // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits:
+            // XXX this is sort of an invalid / out-of-bounds result: to set or not to set these bits:
             stats->results |= PS_STAT_ROBUST_MEDIAN;
             stats->results |= PS_STAT_ROBUST_STDEV;
@@ -791 +797 @@
         psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
 
-        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can not specify the bin
-        // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize.  If
-        // we get here, we know that binSize != 0.0.
-        long numBins = (max - min) / binSize; // Number of bins
+        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can
+        // not specify the bin size precisely since the argument to psHistogramAlloc() is
+        // the number of bins, not the binSize.  If we get here, we know that binSize !=
+        // 0.0.  We can also have a floating-point round-off error such that the last bin
+        // of the histogram does not correspond exactly with the value of 'max'.  Let's be
+        // a bit generous and extend the histogram by two bins in either direction
+        long numBins = 4 + (max - min) / binSize; // Number of bins
         psTrace(TRACE, 6, "Numbins is %ld\n", numBins);
         psTrace(TRACE, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
         // Generate the histogram
-        histogram = psHistogramAlloc(min, max, numBins);
+        histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins);
         // XXXXX we need to consider this step if errors -> variance
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
-            // if psVectorHistogram returns false, we have a programming error
+            // if psVectorHistogram returns false, we have a programming error
             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for robust statistics.\n");
-            psFree(histogram);
-            psFree(cumulative);
-            psFree(statsMinMax);
-            psFree(mask);
-
-            return false;
+            psFree(histogram);
+            psFree(cumulative);
+            psFree(statsMinMax);
+            psFree(mask);
+            return false;
         }
         if (psTraceGetLevel("psLib.math") >= 8) {
             PS_VECTOR_PRINT_F32(histogram->bounds);
             PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+
+        // perversity check: if most of the values land in a single bin, then we probably
+        // have a perverse case (eg, small number of points at extremely large / small
+        // values; nearly bi-modal distribution).  if so, keep only points within 5? 10?
+        // bins of that excess bin:
+        int nMaxBin = 0;
+        int iMaxBin = 0;
+        for (long i = 1; i < histogram->nums->n; i++) {
+            if (histogram->nums->data.F32[i] > nMaxBin) {
+                nMaxBin = histogram->nums->data.F32[i];
+                iMaxBin = i;
+            }
+        }
+        if (nMaxBin > numValid / 2) {
+            float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize;
+            float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize;
+            int nInvalid = 0;
+            for (long i = 0; i < myVector->n; i++) {
+                // skip the already-masked values
+                if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue;
+                bool invalid = false;
+                invalid |= (myVector->data.F32[i] <= minKeep);
+                invalid |= (myVector->data.F32[i] >= maxKeep);
+                invalid |= (!isfinite(myVector->data.F32[i]));
+                if (!invalid) continue;
+                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
+                nInvalid ++;
+            }
+
+            if (nInvalid) {
+              psTrace(TRACE, 6, "data is concentrated in a single bin, masking %d extreme outliers and retrying\n", nInvalid);
+              psFree(histogram);
+              psFree(cumulative);
+              histogram = NULL;
+              cumulative = NULL;
+              continue;
+            }
+            // if we did not mask anything, give up.
         }
 
@@ -838 +885 @@
 
         // 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);
-        // fprintf (stderr, "robustBin : %f vs %f\n", robustBin, stats->robustMedian);
-
-        // convert bin to bin value: this is the robust histogram median.
+        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);
+
+        // convert bin to bin value: this is the robust histogram median.
         if (isnan(stats->robustMedian)) {
             COUNT_WARNING(10, 100, "Failed to fit a quadratic and calculate the 50-percent position.\n");
@@ -1007 +1054 @@
     stats->robustN50 = N50;
     psTrace(TRACE, 6, "The robustN50 is %ld.\n", N50);
+    psTrace(TRACE, 6, "The robust median and stdev are %f, %f\n", stats->robustMedian, stats->robustStdev);
 
     // Clean up
@@ -1054 +1102 @@
         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
-            return false;
-        }
+            return false;
+        }
     }
 
@@ -1107 +1155 @@
         psHistogram *histogram = psHistogramAlloc(min, max, numBins); // A new histogram (without outliers)
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
-            // if psVectorHistogram returns false, we have a programming error
+            // if psVectorHistogram returns false, we have a programming error
             psError(PS_ERR_UNKNOWN, false, "Unable to generate histogram for fitted statistics.\n");
             psFree(histogram);
@@ -1172 +1220 @@
             PS_VECTOR_PRINT_F32(y);
         }
-        
-        // psMinimizeLMChi2 can return false for bad data as well as for serious failures
+
+        // psMinimizeLMChi2 can return false for bad data as well as for serious failures
         if (!psMinimizeLMChi2(minimizer, NULL, params, NULL, x, y, NULL, minimizeLMChi2Gauss1D)) {
             psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
@@ -1235 +1283 @@
         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
-            return false;
-        }
+            return false;
+        }
     }
 
@@ -1355 +1403 @@
         // psVectorInit (fitMask, 0);
 
-        // XXX not sure if these should result in errors or not...
+        // XXX not sure if these should result in errors or not...
         if (!psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x)) {
             psError(PS_ERR_UNKNOWN, false, "Failed to fit a gaussian to the robust histogram.\n");
@@ -1432 +1480 @@
         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
-            return false;
-        }
+            return false;
+        }
     }
 
@@ -1729 +1777 @@
         if (!vectorRobustStats(myVector, errors, mask, maskVal, stats)) {
             psError(PS_ERR_UNKNOWN, false, "failure to measure robust stats\n");
-            return false;
-        }
+            return false;
+        }
     }
 
     // If the mean is NAN, then generate a warning and set the stdev to NAN.
-    if (isnan(stats->robustMedian)) {
-        stats->fittedStdev = NAN;
-        stats->fittedStdev = NAN;
-        return true;
-    }
+    if (isnan(stats->robustMedian)) goto escape;
 
     float guessStdev = stats->robustStdev;  // pass the guess sigma
@@ -1770 +1814 @@
             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
             psFree(statsMinMax);
-            stats->fittedStdev = NAN;
-            stats->fittedStdev = NAN;
+            goto escape;
+        }
+
+        // If all data points have the same value, then we set the appropriate members of stats and return.
+        if (fabs(max - min) <= FLT_EPSILON) {
+            COUNT_WARNING(10, 100, "All data points have the same value: %f.\n", min);
+            stats->fittedMean = min;
+            stats->fittedStdev = 0.0;
+            stats->results |= PS_STAT_FITTED_MEAN_V4;
+            stats->results |= PS_STAT_FITTED_STDEV_V4;
             return true;
         }
@@ -1789 +1841 @@
             psFree(histogram);
             psFree(statsMinMax);
-            stats->fittedStdev = NAN;
-            stats->fittedStdev = NAN;
-            return true;
+            goto escape;
         }
         if (psTraceGetLevel("psLib.math") >= 8) {
@@ -1823 +1873 @@
             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
             psFree(statsMinMax);
-            stats->fittedStdev = NAN;
-            stats->fittedStdev = NAN;
-            return true;
+            psFree(histogram);
+            goto escape;
         }
 
@@ -1886 +1935 @@
 
             // 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)
-            // we probably should be able to handle the data errors gracefully
+            // XXX this fit may fail with an error for an ill-conditioned matrix (bad data)
+            // we probably should be able to handle the data errors gracefully
             psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
             bool status = psVectorFitPolynomial1D (poly, NULL, 0, y, NULL, x);
@@ -1894 +1943 @@
 
             if (!status) {
+                psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
                 psFree(poly);
                 psFree(histogram);
                 psFree(statsMinMax);
-                stats->fittedStdev = NAN;
-                stats->fittedStdev = NAN;
-                return true;
+                goto escape;
             }
 
@@ -1919 +1967 @@
 
                 COUNT_WARNING(10, 100, "fit did not converge\n");
-                stats->fittedStdev = NAN;
-                stats->fittedStdev = NAN;
-                return true;
+                goto escape;
             }
 
@@ -1987 +2033 @@
 
             if (!status) {
+                psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
                 psFree(poly);
                 psFree(histogram);
                 psFree(statsMinMax);
-                stats->fittedStdev = NAN;
-                stats->fittedStdev = NAN;
-                return true;
+                goto escape;
             }
 
@@ -2038 +2083 @@
     psTrace(TRACE, 6, "The fitted stdev is %f.\n", stats->fittedStdev);
 
+    stats->results |= PS_STAT_FITTED_MEAN_V4;
+    stats->results |= PS_STAT_FITTED_STDEV_V4;
+
+    return true;
+
+escape:
+    stats->fittedMean = NAN;
+    stats->fittedStdev = NAN;
     stats->results |= PS_STAT_FITTED_MEAN_V4;
     stats->results |= PS_STAT_FITTED_STDEV_V4;
@@ -2896 +2949 @@
 *****************************************************************************/
 static psF32 fitQuadraticSearchForYThenReturnBin(const psVector *xVec,
-                                                 psVector *yVec,
-                                                 psS32 binNum,
-                                                 psF32 yVal
+                                                 psVector *yVec,
+                                                 psS32 binNum,
+                                                 psF32 yVal
     )
 {
@@ -2980 +3033 @@
         }
 
-        // 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;
+        // 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;
     } 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];       
-            tmpFloat = (yVal - yVec->data.F32[0]) * (dX / dY) + xVec->data.F32[0];
+            // 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];           
-            tmpFloat = (yVal - yVec->data.F32[binNum-1]) * (dX / dY) + xVec->data.F32[binNum-1];
-        } 
+            // 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];
+            }
+        }
     }