Changeset 6437

Timestamp:

Feb 16, 2006, 2:56:48 PM (20 years ago)

Author:

gusciora

Message:

* empty log message *

Location:

trunk/psLib/src/math

Files:

• psConstants.h (modified) (5 diffs)
• psMathUtils.h (modified) (2 diffs)
• psMinimizeLMM.h (modified) (2 diffs)
• psPolynomial.c (modified) (6 diffs)
• psSpline.c (modified) (6 diffs)
• psStats.c (modified) (17 diffs)

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.84 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-26 21:10:22 $
+ *  @version $Revision: 1.85 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 00:56:48 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -169 +169 @@
 }
 
-// Return an error if the arg is lies outside the supplied range.
+// Return an error if the arg lies outside the supplied range.
 #define PS_ASSERT_FLOAT_WITHIN_RANGE(NAME, LOWER, UPPER, RVAL) \
 if ((NAME) < (LOWER) || (NAME) > (UPPER)) { \
@@ -186 +186 @@
 }
 
-// Return an error if the arg lies outside the supplied range
 #define PS_ASSERT_LONG_WITHIN_RANGE(NAME, LOWER, UPPER, RVAL) \
 if ((NAME) < (LOWER) || (NAME) > (UPPER)) { \
@@ -519 +518 @@
 } \
 
-/* XXX: This is not correct
-    #define PS_POLY_2D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
-    static psPolynomial2D *(NAME) = NULL; \
-    if ((NAME) == NULL) { \
-        (NAME) = psPolynomial2DAlloc(ORDER, TYPE); \
-        p_psMemSetPersistent((NAME), true); \
-    } \
- 
-    #define PS_POLY_3D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
-    static psPolynomial3D *(NAME) = NULL; \
-    if ((NAME) == NULL) { \
-        (NAME) = psPolynomial3DAlloc(ORDER, TYPE); \
-        p_psMemSetPersistent((NAME), true); \
-    } \
- 
-    #define PS_POLY_4D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
-    static psPolynomial4D *(NAME) = NULL; \
-    if ((NAME) == NULL) { \
-        (NAME) = psPolynomial4DAlloc(ORDER, TYPE); \
-        p_psMemSetPersistent((NAME), true); \
-    } \
-*/
-
 #define PS_POLY_1D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
 static psPolynomial1D *(NAME) = NULL; \
@@ -548 +524 @@
     p_psMemSetPersistent((NAME), true); \
 } \
-
-/* XXX: This is not correct
-    #define PS_POLY_2D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
-    static psPolynomial2D *(NAME) = NULL; \
-    if ((NAME) == NULL) { \
-        (NAME) = psPolynomial2DAlloc(ORDER, TYPE); \
-        p_psMemSetPersistent((NAME), true); \
-    } \
- 
-    #define PS_POLY_3D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
-    static psPolynomial3D *(NAME) = NULL; \
-    if ((NAME) == NULL) { \
-        (NAME) = psPolynomial3DAlloc(ORDER, TYPE); \
-        p_psMemSetPersistent((NAME), true); \
-    } \
- 
-    #define PS_POLY_4D_D_DECLARE_ALLOC_STATIC(NAME, ORDER, TYPE) \
-    static psPolynomial4D *(NAME) = NULL; \
-    if ((NAME) == NULL) { \
-        (NAME) = psPolynomial4DAlloc(ORDER, TYPE); \
-        p_psMemSetPersistent((NAME), true); \
-    } \
-*/
 
 #define PS_POLY_PRINT_1D(NAME) \

trunk/psLib/src/math/psMathUtils.h

@@ -7 +7 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-02 21:09:07 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 00:56:48 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -51 +51 @@
 );
 
-// XXX: Create a single, generic, version of the vector normalize function.
-// XXX: Ask IfA for a public psLib function.
-
 psS32 p_psNormalizeVectorRange(psVector* myData,
                                psF64 outLow,

trunk/psLib/src/math/psMinimizeLMM.h

@@ -8 +8 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-27 20:08:58 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 00:56:48 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -115 +115 @@
     const psArray *x,                  ///< Measurement ordinates of multiple vectors
     const psVector *y,                 ///< Measurement coordinates
-    const psVector *yErr,              ///< Errors in the measurement coordinates
+    const psVector *yWt,               ///< Errors in the measurement coordinates
     psMinimizeLMChi2Func func          ///< Specified function
 );

trunk/psLib/src/math/psPolynomial.c

@@ -4 +4 @@
 *         routines.
 *
-*  This file will hold the functions for allocated, freeing, and evaluating
+*  This file will hold the routiness for allocating, freeing, and evaluating
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.142 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-07 23:36:15 $
+*  @version $Revision: 1.143 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-17 00:56:48 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -159 +159 @@
 used frequently and the data structure created here does not contain the
 outer coefficients of the Chebyshev polynomials.
- 
  *****************************************************************************/
 psPolynomial1D **p_psCreateChebyshevPolys(psS32 numPolys)
@@ -165 +164 @@
     PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(numPolys, 1, NULL);
 
-    if (0) {
-        psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
-        for (psS32 i = 0; i < numPolys; i++) {
-            chebPolys[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, i);
-            chebPolys[i]->coeff[i] = 1;
-        }
-        return (chebPolys);
-    } else {
-        psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
-        for (psS32 i = 0; i < numPolys; i++) {
-            chebPolys[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, i);
-        }
-
-        // Create the Chebyshev polynomials.
-        // Polynomial i has i-th order.
-        chebPolys[0]->coeff[0] = 1.0;
-        if (numPolys >= 2) {
-            chebPolys[1]->coeff[1] = 1.0;
-
-            for (psS32 i = 2; i < numPolys; i++) {
-                for (psS32 j = 0; j < chebPolys[i - 1]->nX+1; j++) {
-                    chebPolys[i]->coeff[j + 1] = 2.0 * chebPolys[i - 1]->coeff[j];
-                }
-                for (psS32 j = 0; j < chebPolys[i - 2]->nX+1; j++) {
-                    chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
-                }
-            }
-        }
-
-        if (psTraceGetLevel(__func__) >= 6) {
-            for (psS32 j = 0; j < numPolys; j++) {
-                PS_POLY_PRINT_1D(chebPolys[j]);
-            }
-        }
-        return (chebPolys);
-    }
+    psPolynomial1D **chebPolys = (psPolynomial1D **) psAlloc(numPolys * sizeof(psPolynomial1D *));
+    for (psS32 i = 0; i < numPolys; i++) {
+        chebPolys[i] = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, i);
+    }
+
+    // Create the Chebyshev polynomials.
+    // Polynomial i has i-th order.
+    chebPolys[0]->coeff[0] = 1.0;
+    if (numPolys >= 2) {
+        chebPolys[1]->coeff[1] = 1.0;
+
+        for (psS32 i = 2; i < numPolys; i++) {
+            for (psS32 j = 0; j < chebPolys[i - 1]->nX+1; j++) {
+                chebPolys[i]->coeff[j + 1] = 2.0 * chebPolys[i - 1]->coeff[j];
+            }
+            for (psS32 j = 0; j < chebPolys[i - 2]->nX+1; j++) {
+                chebPolys[i]->coeff[j] -= chebPolys[i - 2]->coeff[j];
+            }
+        }
+    }
+
+    if (psTraceGetLevel(__func__) >= 6) {
+        for (psS32 j = 0; j < numPolys; j++) {
+            PS_POLY_PRINT_1D(chebPolys[j]);
+        }
+    }
+    return (chebPolys);
 }
 
@@ -223 +213 @@
     for (loop_x = 0; loop_x < poly->nX+1; loop_x++) {
         if (poly->mask[loop_x] == 0) {
-            // XXX: If you set the tracelevel to 10 here, and later set the tracelevel to
-            // 2 or higher in the test code, you get seg faults.
             psTrace(__func__, 8,
                     "polysum+= sum*coeff [%lf+= (%lf * %lf)\n", polySum, xSum, poly->coeff[loop_x]);
@@ -238 +226 @@
 // XXX: You can do this without having to psAlloc() vector d.
 // XXX: How does the mask vector effect Crenshaw's formula?
-// XXX: We assume that x is scaled between -1.0 and 1.0;
+// NOTE: We assume that x is scaled between -1.0 and 1.0;
 // XXX: Create a faster version for low-order Chebyshevs.
 static psF64 chebPolynomial1DEval(
@@ -298 +286 @@
         psFree(d);
     } else {
-        // This is old code that does not use Clenshaw's formula.  Get rid of it.
+        // XXX: This is old code that does not use Clenshaw's formula.  Get rid of it.
         psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(1 + poly->nX);
 

trunk/psLib/src/math/psSpline.c

@@ -6 +6 @@
 *  This file contains the routines that allocate, free, and evaluate splines.
 *
-*  @version $Revision: 1.135 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-02 21:09:08 $
+*  @version $Revision: 1.136 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-17 00:56:48 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -101 +101 @@
  
 The first and second derivatives at the endpoints were previously undefined in
-the SDR.  From bugzilla #???, they are required to be 0.0, impementing natural
+the SDR.  From bugzilla #???, they are required to be 0.0, implementing natural
 splines.
  
@@ -228 +228 @@
     // The following code ensures that xPtr and yPtr points to a psF32 psVector.
     //
-    // XXX: When you debug, use the vector copy and create routines here:
+    // XXX: Use the vector copy and create routines here:
     //
 
@@ -345 +345 @@
 /*****************************************************************************
 psSpline1DEval(): this routine takes an existing spline of arbitrary order
-and an independent x value.  Each determines which spline that x corresponds
+and an independent x value.  It determines which spline that x corresponds
 to by doing a bracket disection on the knots of the spline data structure
 (vectorBinDisectF32()).  Then it evaluates the spline at that x location
@@ -360 +360 @@
     float x)
 {
-    psTrace(__func__, 3, "---- %s(%f) begin ----\n", __func__, x);
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(spline, NAN);
     PS_ASSERT_INT_NONNEGATIVE(spline->n, NAN);
@@ -366 +366 @@
 
     psS32 n = spline->n;
+    if ((x < spline->knots->data.F32[0]) || (x > spline->knots->data.F32[n-1])) {
+        // If x is outside the range of spline->knots, generate a warning
+        // message, then return the left, or right, endpoint.
+        psLogMsg(__func__, PS_LOG_WARN,
+                 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).",
+                 x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n);
+
+        psS32 binNum = (x < spline->knots->data.F32[0]) ? 0 : n-1;
+        psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+        return(psPolynomial1DEval(spline->spline[binNum], x));
+    }
+
     psScalar tmpScalar;
     tmpScalar.type.type = PS_TYPE_F32;
     tmpScalar.data.F32 = x;
     psS32 binNum = p_psVectorBinDisect(spline->knots, &tmpScalar);
-
     if (binNum < 0) {
-        //
-        // If x is outside the range of spline->knots, generate a warning
-        // message, then return the left, or right, endpoint.
-        //
-        psLogMsg(__func__, PS_LOG_WARN,
-                 "psSpline1DEval(): x ordinate (%f) is outside the spline range (%f - %f) (%d).",
-                 x, spline->knots->data.F32[0], spline->knots->data.F32[n-1], n);
-
-        if (x < spline->knots->data.F32[0]) {
-            psF32 rcF32 = psPolynomial1DEval(spline->spline[0], x);
-            psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, rcF32);
-            return(rcF32);
-        } else if (x > spline->knots->data.F32[n-1]) {
-            psF32 rcF32 = psPolynomial1DEval(spline->spline[n-1], x);
-            psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, rcF32);
-            return(rcF32);
-        }
-    }
-
-    psF32 rcF32 = psPolynomial1DEval(spline->spline[binNum], x);
-    psTrace(__func__, 3, "---- %s(%f) end ----\n", __func__, rcF32);
-    return(rcF32);
+        psError(PS_ERR_UNKNOWN, false, "Could not perform bin dissection on spline->knots.\n");
+        return(NAN);
+    }
+
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return(psPolynomial1DEval(spline->spline[binNum], x));
 }
 

trunk/psLib/src/math/psStats.c

@@ -16 +16 @@
  * use ->min and ->max (PS_STAT_USE_RANGE)
  *
- *  @version $Revision: 1.167 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-03 22:05:22 $
+ *  @version $Revision: 1.168 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 00:56:48 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -440 +440 @@
 unmasked element within the specified min/max range).  Otherwise, return
 "false".
- 
-XXX: Can you use psVectorCountPixelMask here?
  *****************************************************************************/
-bool p_psVectorCheckNonEmpty(const psVector* myVector,
-                             const psVector* maskVector,
-                             psU32 maskVal,
-                             psStats* stats)
+bool p_psVectorCheckNonEmpty(
+    const psVector* myVector,
+    const psVector* maskVector,
+    psU32 maskVal,
+    psStats* stats)
 {
     psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
@@ -496 +495 @@
 number of non-masked pixels in the vector that fall within the min/max
 range, if specified.
- 
-XXX: Can you use psVectorCountPixelMask here?
  *****************************************************************************/
 psS32 p_psVectorNValues(const psVector* myVector,
@@ -644 +641 @@
 robustHistogram with a Gaussian of width sigma.  It returns a psVector of the
 smoothed data.
- 
-XXX: Only PS_TYPE_F32 is supported.
- 
-XXX: Write a general routine which smoothes a psVector.  This routine should
-call that.  Is that possible?
  *****************************************************************************/
 psVector *p_psVectorSmoothHistGaussian(
@@ -985 +977 @@
     -1: error
     -2: warning
- 
-XXX: Use static vectors for tmpMask.
  *****************************************************************************/
 psS32 p_psVectorClippedStats(const psVector* myVector,
@@ -1224 +1214 @@
 
         //
-        // Ensure that the y values are monotonic.
-        //
-        // XXX: This routine should probably be rewritten in a more general fashion
-        // so that the following checks are not necessary.
+        // Ensure that the y value lies within range of the y values.
         //
         if (! (((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2])) ||
@@ -1236 +1223 @@
         }
 
+        //
+        // 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]))) {
@@ -1289 +1279 @@
         } else if (binNum == (xVec->n - 2)) {
             // XXX: Is this right?
-            tmpFloat = 0.5 * (xVec->data.F32[binNum] +
-                              xVec->data.F32[binNum + 1]);
+            tmpFloat = 0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum + 1]);
         }
     }
@@ -1677 +1666 @@
     tmpScalar.data.F32 = totalDataPoints * 0.25f;
     if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
-        // XXX: Special case where its in first bin.  Must code last bin.
+        // Special case where its in first bin.  The last bin should be handled automatically.
         binLo25 = 0;
     } else {
@@ -1684 +1673 @@
     tmpScalar.data.F32 = totalDataPoints * 0.75f;
     if (tmpScalar.data.F32 <= cumulativeRobustHistogram->nums->data.F32[0]) {
-        // XXX: Special case where its in first bin.  Must code last bin.
+        // Special case where its in first bin.  The last bin should be handled automatically.
         binHi25 = 0;
     } else {
@@ -1704 +1693 @@
     // Interpolate to find these two positions exactly: these are the upper
     // and lower quartile positions.
-    // XXX: Check for errors.
     //
     psF32 binLo25F32 = fitQuadraticSearchForYThenReturnX(
@@ -1711 +1699 @@
                            binLo25,
                            totalDataPoints * 0.25f);
+
     psF32 binHi25F32 = fitQuadraticSearchForYThenReturnX(
                            *(psVector* *)&cumulativeRobustHistogram->bounds,
@@ -1716 +1705 @@
                            binHi25,
                            totalDataPoints * 0.75f);
+
+    if (isnan(binLo25F32) || isnan(binHi25F32)) {
+        psError(PS_ERR_UNKNOWN, false, "could not determine the robustUQ: fitQuadraticSearchForYThenReturnX() returned a NAN.\n");
+        psFree(tmpStatsMinMax);
+        psFree(robustHistogram);
+        psFree(cumulativeRobustHistogram);
+        psFree(tmpMaskVec);
+        psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+        return(1);
+    }
+
     stats->robustLQ = binLo25F32;
     stats->robustUQ = binHi25F32;
     psTrace(__func__, 6, "The 25 and 75 percent data point exact positions are (%f, %f).\n", binLo25F32, binHi25F32);
-
     psS32 N50 = 0;
     for (psS32 i = 0 ; i < myVector->n ; i++) {
@@ -1863 +1862 @@
         }
 
-        // XXX: Use the min/max routines for this
-        psF32 minY = FLT_MAX;
-        psF32 maxY = -FLT_MAX;
-        for (psS32 i = 0 ; i < y->n ; i++) {
-            if (y->data.F32[i] > maxY) {
-                maxY = y->data.F32[i];
-            }
-            if (y->data.F32[i] < minY) {
-                minY = y->data.F32[i];
-            }
-        }
+        rc = p_psVectorMin(y, NULL, 0, tmpStatsMinMax);
+        rc|= p_psVectorMax(y, NULL, 0, tmpStatsMinMax);
+        if ((rc != 0) || isnan(tmpStatsMinMax->min) || isnan(tmpStatsMinMax->max)) {
+            psError(PS_ERR_UNKNOWN, false, "Failed to calculate the min/max of the input vector.\n");
+            psFree(tmpMaskVec);
+            psFree(robustHistogram);
+            psFree(cumulativeRobustHistogram);
+            psFree(tmpStatsMinMax);
+            psTrace(__func__, 4, "---- %s(1) end  ----\n", __func__);
+            return(1);
+        }
+
         //
         // Normalize y to [0.0:1.0] (since the psMinimizeLMChi2Gauss1D() functions is [0.0:1.0])
@@ -1879 +1879 @@
         //
         for (psS32 i = 0 ; i < y->n ; i++) {
-            y->data.F32[i]= (y->data.F32[i] - minY) / (maxY - minY);
+            y->data.F32[i]= (y->data.F32[i] - tmpStatsMinMax->min) / (tmpStatsMinMax->max - tmpStatsMinMax->min);
         }
 
@@ -2295 +2295 @@
                         }
                     } else {
-                        // XXX: This if-statement really shouldn't be necessary.
+                        // This if-statement really shouldn't be necessary.
                         // However, due to numerical lack of precision, we
                         // occasionally produce a binNum outside the range.
@@ -2435 +2435 @@
 XXX: Should we free stats if the asserts fail?
  *****************************************************************************/
-psStats* psVectorStats(psStats* stats,
-                       const psVector* in,
-                       const psVector* errors,
-                       const psVector* mask,
-                       psMaskType maskVal)
+psStats* psVectorStats(
+    psStats* stats,
+    const psVector* in,
+    const psVector* errors,
+    const psVector* mask,
+    psMaskType maskVal)
 {
     psTrace(__func__, 3,"---- %s() begin  ----\n", __func__);