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

Timestamp:

Sep 11, 2013, 6:12:49 PM (13 years ago)

Author:

eugene

Message:

use linear instead of quadratic fit for UQ and LQ as well (to avoid ill-conditioned matrices)

File:

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

trunk/psLib/src/math/psStats.c

@@ -174 +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);
 
 /******************************************************************************
@@ -918 +919 @@
         // 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.
+# if (0)
         { // Quick 5-point linear fit
           double Sx = (cumulative->nums->data.F32[binMedian - 2] + cumulative->nums->data.F32[binMedian - 1] +
@@ -937 +939 @@
           stats->robustMedian = linearMedian;
         }
+# endif
+        stats->robustMedian = fitLinearSearchForYThenReturnBin(cumulative->bounds, cumulative->nums, binMedian, totalDataPoints/2.0);
 
         // convert bin to bin value: this is the robust histogram median.
@@ -1108 +1112 @@
     // 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;
     }
@@ -2036 +2040 @@
 // other private functions used above
 
+# if (0)
 static psF32 QuadraticInverse(psF32 a,
                               psF32 b,
@@ -2057 +2062 @@
     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)
@@ -2313 +2334 @@
     return tmpFloat;
 }
-# endif
 
 /******************************************************************************
@@ -2484 +2504 @@
     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)
+    double Sx = (yVec->data.F32[binNum - 2] + yVec->data.F32[binNum - 1] +
+                 yVec->data.F32[binNum - 0] +
+                 yVec->data.F32[binNum + 1] + yVec->data.F32[binNum + 2]);
+    double Sy = (xVec->data.F32[binNum - 2] + xVec->data.F32[binNum - 1] +
+                 xVec->data.F32[binNum - 0] +
+                 xVec->data.F32[binNum + 1] + xVec->data.F32[binNum + 2]);
+    double Sxx = (PS_SQR(yVec->data.F32[binNum - 2]) + PS_SQR(yVec->data.F32[binNum - 1]) +
+                  PS_SQR(yVec->data.F32[binNum - 0]) +
+                  PS_SQR(yVec->data.F32[binNum + 1]) + PS_SQR(yVec->data.F32[binNum + 2]));
+    double Sxy = (xVec->data.F32[binNum - 2] * yVec->data.F32[binNum - 2] +
+                  xVec->data.F32[binNum - 1] * yVec->data.F32[binNum - 1] +
+                  xVec->data.F32[binNum - 0] * yVec->data.F32[binNum - 0] +
+                  xVec->data.F32[binNum + 1] * yVec->data.F32[binNum + 1] +
+                  xVec->data.F32[binNum + 2] * yVec->data.F32[binNum + 2]);
+    double value = ((Sy * Sxx - Sx * Sxy) + (5 * Sxy - Sx * Sy) * yVal)/(5 * Sxx - Sx * Sx);
+
+    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;
+                
+#if (CZW)
+        printf("   internal median: %f %f\n",tmpFloat,binValue);
+#endif
+        
+    } 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
+}