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

Timestamp:

Feb 3, 2006, 12:05:22 PM (20 years ago)

Author:

gusciora

Message:

Misc code cleaning

File:

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

trunk/psLib/src/math/psStats.c

@@ -16 +16 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.166 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-03 01:30:14 $
+ *  @version $Revision: 1.167 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-03 22:05:22 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -303 +303 @@
 max of the input vector.  If there was a problem with the max calculation,
 this routine sets stats->max to NAN.
- 
-XXX: Do we need to factor errors into it?
  *****************************************************************************/
 psS32 p_psVectorMax(const psVector* myVector,
@@ -549 +547 @@
 median of the input vector.  Returns true on success (including if there were
 no valid input vector elements).
- 
-XXX: Use static vectors for sort arrays.
  *****************************************************************************/
 bool p_psVectorSampleMedian(const psVector* myVector,
@@ -783 +779 @@
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    psVector* unsortedVector = NULL;    // Temporary vector
-    psVector* sortedVector = NULL;      // Temporary vector
-    psS32 i = 0;                  // Loop index variable
-    psS32 count = 0;              // # of points in this mean.
-    psS32 nValues = 0;            // # data points
+    psVector* unsortedVector = NULL;
+    psVector* sortedVector = NULL;
+    psS32 i = 0;
+    psS32 count = 0;
+    psS32 nValues = 0;
 
     // Determine how many data points fit inside this min/max range
-    // and are not maxed, IF the maskVector is not NULL.
+    // and are not masked, if the maskVector is not NULL.
     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
 
@@ -796 +792 @@
     unsortedVector = psVectorAlloc(nValues, PS_TYPE_F32);
 
-    // Determine if we must only use data points within a min/max range.
     if (stats->options & PS_STAT_USE_RANGE) {
         // Store all non-masked data points within the min/max range
@@ -1001 +996 @@
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
     psTrace(__func__, 4, "Trace level is %d\n", psTraceGetLevel(__func__));
-    psF32 clippedMean = 0.0;    // self-explanatory
-    psF32 clippedStdev = 0.0;   // self-explanatory
-    psVector* tmpMask = NULL;   // Temporary vector for masks during iterations.
-    static psStats *statsTmp = NULL;   // Temporary psStats struct.
-    psS32 rc = 0;               // Return code.
+    psF32 clippedMean = 0.0;
+    psF32 clippedStdev = 0.0;
+    psVector* tmpMask = NULL;
+    static psStats *statsTmp = NULL;
+    psS32 rc = 0;
 
     // Ensure that stats->clipIter is within the proper range.
@@ -1189 +1184 @@
 parameter).
  
-XXX: After you fit the polynomial, solve for X analytically.
- 
-XXX: the vectors do not have to be the same length.  Must insert the proper
-tests to ensure that binNum is within acceptable ranges for both vectors.
- 
-XXX: This currently assumes that the three points are monotonically increasing
-or decreasing: so, it works for the cumulative histogram vectors, but not for
-arbitrary vectors.  We should probably test that condition.
 *****************************************************************************/
 psF32 fitQuadraticSearchForYThenReturnX(
@@ -1221 +1208 @@
     psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     psVector *y = psVectorAlloc(3, PS_TYPE_F64);
-    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
-
     psF32 tmpFloat = 0.0f;
 
-    if ((binNum > 0) && (binNum < (yVec->n - 2))) {
+    if ((binNum >= 1) && (binNum < (yVec->n - 2)) && (binNum < (xVec->n - 2))) {
         // The general case.  We have all three points.
         x->data.F64[0] = (psF64) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum]));
@@ -1234 +1219 @@
         y->data.F64[2] = yVec->data.F32[binNum + 1];
         psTrace(__func__, 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(__func__, 6, "x vec is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
-        psTrace(__func__, 6, "y vec is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
+                xVec->data.F32[binNum], xVec->data.F32[binNum+1], xVec->data.F32[binNum+2]);
+        psTrace(__func__, 6, "x data is (%f %f %f)\n", x->data.F64[0], x->data.F64[1], x->data.F64[2]);
+        psTrace(__func__, 6, "y data is (%f %f %f)\n", y->data.F64[0], y->data.F64[1], y->data.F64[2]);
 
         //
@@ -1246 +1229 @@
         // so that the following checks are not necessary.
         //
-        if (y->data.F64[0] < y->data.F64[1]) {
-            if (!(y->data.F64[1] <= y->data.F64[2])) {
-                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
-                psFree(x);
-                psFree(y);
-                psFree(yErr);
-                psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
-                return(NAN);
-            }
-            // Ensure that yVal is within the range of the bins we are using.
-            if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                        PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
-                        (psF64)yVal, y->data.F64[2], y->data.F64[0]);
-            }
-        } else if (y->data.F64[0] > y->data.F64[1]) {
-            if (!(y->data.F64[1] >= y->data.F64[2])) {
-                psError(PS_ERR_UNKNOWN, true, "This routine must be called with montically increasing or decreasing data points.\n");
-                psFree(x);
-                psFree(y);
-                psFree(yErr);
-                psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
-                return(NAN);
-            }
-            // Ensure that yVal is within the range of the bins we are using.
-            if (!((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) {
-                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
-                        PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
-                        (psF64)yVal, y->data.F64[2], y->data.F64[0]);
-            }
-        }
-
-        yErr->data.F64[0] = 1.0;
-        yErr->data.F64[1] = 1.0;
-        yErr->data.F64[2] = 1.0;
+        if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
+                ((y->data.F64[2] <= yVal) && (yVal <= y->data.F64[0]))) ) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    PS_ERRORTEXT_psStats_YVAL_OUT_OF_RANGE,
+                    (psF64)yVal, y->data.F64[0], y->data.F64[2]);
+        }
+
+        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);
+            psTrace(__func__, 5, "---- %s() end ----\n", __func__);
+            return(NAN);
+        }
 
         // Determine the coefficients of the polynomial.
         psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2);
-        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, yErr, x);
+        myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x);
         if (myPoly == NULL) {
-            psError(PS_ERR_UNEXPECTED_NULL,
-                    false,
+            psError(PS_ERR_UNEXPECTED_NULL, false,
                     PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT);
             psFree(x);
             psFree(y);
-            psFree(yErr);
             psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
             return(NAN);
@@ -1309 +1270 @@
         if (isnan(tmpFloat)) {
             psError(PS_ERR_UNEXPECTED_NULL,
-                    false,
-                    PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
+                    false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN);
             psFree(x);
             psFree(y);
-            psFree(yErr);
             psTrace(__func__, 5, "---- %s(NAN) end ----\n", __func__);
             return(NAN);
         }
-
     } else {
-        // The special case where we have two points only at the beginning of
-        // the vectors x and y.
+        // 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.
             tmpFloat = 0.5 * (xVec->data.F32[binNum] +
                               xVec->data.F32[binNum + 1]);
@@ -1339 +1297 @@
     psFree(x);
     psFree(y);
-    psFree(yErr);
 
     psTrace(__func__, 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
@@ -1542 +1499 @@
         // ADD: Step 3.
         // Interpolate to the exact 50% position: this is the robust histogram median.
-        // XXX: Check return codes.
         //
         stats->robustMedian = fitQuadraticSearchForYThenReturnX(
@@ -1549 +1505 @@
                                   binMedian,
                                   totalDataPoints/2.0);
+        if (isnan(stats->robustMedian)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to fit a quadratic and calculate the 50-percent position.\n");
+            psFree(tmpStatsMinMax);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpMaskVec);
+            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+            return(1);
+        }
         psTrace(__func__, 6, "Current robust median is %f\n", stats->robustMedian);
 
@@ -2324 +2289 @@
                     binNum = (psS32)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize);
                     if (errorsF32 != NULL) {
-                        // XXX: Check return codes.
-                        UpdateHistogramBins(binNum, out,
-                                            inF32->data.F32[i],
-                                            errorsF32->data.F32[i]);
+                        psS32 rc = UpdateHistogramBins(binNum, out, inF32->data.F32[i],
+                                                       errorsF32->data.F32[i]);
+                        if (rc < 0) {
+                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n");
+                        }
                     } else {
                         // XXX: This if-statement really shouldn't be necessary.
@@ -2348 +2314 @@
                     } else {
                         if (errorsF32 != NULL) {
-                            // XXX: Check return codes.
-                            UpdateHistogramBins(binNum, out,
-                                                inF32->data.F32[i],
-                                                errors->data.F32[i]);
+                            psS32 rc = UpdateHistogramBins(binNum, out,
+                                                           inF32->data.F32[i],
+                                                           errors->data.F32[i]);
+                            if (rc < 0) {
+                                psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram bins with the errors vector.\n");
+                            }
                         } else {
                             (out->nums->data.F32[binNum])+= 1.0;