Changeset 6193

Timestamp:

Jan 25, 2006, 2:31:19 PM (20 years ago)

Author:

gusciora

Message:

Several mods mostly to the psPolynomial fitting routines.

Location:

trunk/psLib

Files:

• src/imageops/psImageStats.c (modified) (2 diffs)
• src/math/psConstants.h (modified) (2 diffs)
• src/math/psMinimizePolyFit.c (modified) (12 diffs)
• src/math/psPolynomial.c (modified) (7 diffs)
• src/math/psStats.c (modified) (9 diffs)
• test/math/Makefile.am (modified) (3 diffs)
• test/math/tst_psPolyFit1D.c (modified) (3 diffs)

trunk/psLib/src/imageops/psImageStats.c

@@ -9 +9 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.87 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 22:25:31 $
+ *  @version $Revision: 1.88 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-26 00:31:19 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -432 +432 @@
                 for (y = 0; y < input->numCols; y++) {
                     double pixel;
-                    /*
-                                        if (input->type.type == PS_TYPE_S8) {
-                                            pixel = (double) input->data.S8[x][y];
-                                        } else if (input->type.type == PS_TYPE_U16) {
-                                            pixel = (double) input->data.U16[x][y];
-                                        } else if (input->type.type == PS_TYPE_F32) {
-                                            pixel = (double) input->data.F32[x][y];
-                                        } else if (input->type.type == PS_TYPE_F64) {
-                                            pixel = input->data.F64[x][y];
-                                        }
-                    */
+                    if (0) {
+                        if (input->type.type == PS_TYPE_S8) {
+                            pixel = (double) input->data.S8[x][y];
+                        } else if (input->type.type == PS_TYPE_U16) {
+                            pixel = (double) input->data.U16[x][y];
+                        } else if (input->type.type == PS_TYPE_F32) {
+                            pixel = (double) input->data.F32[x][y];
+                        } else if (input->type.type == PS_TYPE_F64) {
+                            pixel = input->data.F64[x][y];
+                        }
+                    }
                     pixel = nodes->data.F64[x][y];
                     sums[i][j] += pixel * psPolynomial1DEval(chebPolys[i], rScalingFactors[x]) *

trunk/psLib/src/math/psConstants.h

@@ -6 +6 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.82 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 22:25:31 $
+ *  @version $Revision: 1.83 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-26 00:31:19 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -74 +74 @@
 #define PS_ASSERT_INT_NONNEGATIVE(NAME, RVAL) \
 if ((NAME) < 0) { \
-    psError(PS_ERR_BAD_PARAMETER_VALUE, true, \
-            "Error: %s is less than 0.", #NAME); \
+    psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Error: %s is less than 0.", #NAME); \
     return(RVAL); \
 }

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 22:25:31 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-26 00:31:19 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
+ *
+ *  XXX: psMatrixLUSolve() does not return error codes when the results are NANs.
  *
  *  XXX: For clip-fit functions, what should we do if the mask is NULL?
@@ -279 +281 @@
  *****************************************************************************/
 
+static psPolynomial1D* VectorFitPolynomial1DOrd(
+    psPolynomial1D* myPoly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector *f,
+    const psVector *fErr,
+    const psVector *x);
+
+/******************************************************************************
+vectorFitPolynomial1DCheb():  This routine will fit a Chebyshev
+polynomial of degree myPoly->nX to the data points (x, y) and return the
+coefficients of that polynomial.
+ 
+    NOTE: We currently have implemented three algorithms.  This one is
+    non-standard.  It ignores the orthogonal properties of the Chebyshev
+    polys, and their known zero values.  Instead, we first fit a regular
+    ordinary polynomial to the data.  This produces the coefficients for the
+    various x^i terms.  We then "reverse-engineer" those coefficients to
+    determine the coefficients of the Chebyshev polys: basically for each
+    c_ix^i term in the ordinary polynomial, we sum all x^i terms in the
+    Chebshev polys: these must be equal.  This creates a set of nOrder+1
+    linear equations which can be easily solved.  The resulting vector is the
+    coefficients of the Chebyshev polys.
+    
+    This method is significantly slower than the standard NR algorithm.  It
+    was explicitly requested that we not use the NR algorithm.
+ 
+ Matrix A[[][]:
+     Element A[i][j] is the coefficient of the x^i term in the j-th cheby poly.
+ 
+    XXX: This can be improved significantly, performance-wise.  The second set
+    of linear equations which must be "solved" are already in upper-triangular
+    form.  One must only perform the reverse-substitution, LUD decomposition.
+ 
+    XXX: Also, we don't really need to generate the chebPolys data structure.
+    We simply need the matrix which corresponds to a transpose of each Cheby
+    polys coefficients.
+*****************************************************************************/
+static psPolynomial1D *vectorFitPolynomial1DCheb(
+    psPolynomial1D* myPoly,
+    const psVector *mask,
+    psMaskType maskValue,
+    const psVector* y,
+    const psVector* yErr,
+    const psVector* x)
+{
+    PS_ASSERT_POLY_NON_NULL(myPoly, NULL);
+    PS_ASSERT_INT_LARGER_THAN_OR_EQUAL(myPoly->nX, 0, NULL);
+    PS_ASSERT_VECTOR_NON_NULL(y, NULL);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL);
+    if (yErr != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, NULL);
+        PS_ASSERT_VECTOR_TYPE(yErr, PS_TYPE_F64, NULL);
+    }
+    if (x != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL);
+        PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL);
+    }
+    if (mask != NULL) {
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL);
+        PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL);
+    }
+
+    psS32 polyOrder = myPoly->nX;
+    psS32 numPoly = polyOrder + 1;
+
+    //
+    // We first fit an ordinary polynomial to the data.
+    //
+    psPolynomial1D *ordPoly = psPolynomial1DAlloc(polyOrder, PS_POLYNOMIAL_ORD);
+    psPolynomial1D *rc = VectorFitPolynomial1DOrd(ordPoly, mask, maskValue, y, yErr, x);
+    if (rc == NULL) {
+        psError(PS_ERR_UNKNOWN, false, "Could not fit a preliminary polynomial to the data vector.  Returning NULL.\n");
+        psFree(myPoly);
+        return(NULL);
+    }
+
+    //
+    // Create the A-matrix and B-vector which correspond to the linear equations
+    // which will be solved and will then yield the Cheby poly coefficients.
+    //
+    psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(numPoly);
+    psImage *A = psImageAlloc(numPoly, numPoly, PS_TYPE_F64);
+    psVector *B = psVectorAlloc(numPoly, PS_TYPE_F64);
+    for (psS32 i = 0 ; i < numPoly ; i++) {
+        for (psS32 j = 0 ; j < numPoly ; j++) {
+            A->data.F64[i][j] = 0.0;
+        }
+        B->data.F64[i] = ordPoly->coeff[i];
+    }
+
+    for (psS32 i = 0 ; i < numPoly ; i++) {
+        for (psS32 j = 0 ; j < numPoly ; j++) {
+            if (i <= chebPolys[j]->nX)
+                A->data.F64[i][j]+= chebPolys[j]->coeff[i];
+        }
+    }
+    // The following statement is essential.  It derives from (5.8.8) NR.
+    A->data.F64[0][0] = 0.5;
+    psFree(ordPoly);
+    if (psTraceGetLevel(__func__) >= 6) {
+        PS_IMAGE_PRINT_F64(A);
+        PS_VECTOR_PRINT_F64(B);
+    }
+
+    if (0) {
+        // GaussJordan version
+        if (false == psGaussJordan(A, B)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            // the first nTerm entries in B correspond directly to the desired
+            // polynomial coefficients.  this is only true for the 1D case
+            for (psS32 k = 0; k < numPoly; k++) {
+                myPoly->coeff[k] = B->data.F64[k];
+            }
+        }
+    } else {
+        // LUD version of the fit
+        psImage *ALUD = NULL;
+        psVector* outPerm = NULL;
+        psVector* coeffs = NULL;
+
+        ALUD = psImageAlloc(numPoly, numPoly, PS_TYPE_F64);
+        ALUD = psMatrixLUD(ALUD, &outPerm, A);
+        if (ALUD == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+            if (coeffs == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
+                psFree(myPoly);
+                myPoly = NULL;
+            } else {
+                for (psS32 k = 0; k < numPoly; k++) {
+                    myPoly->coeff[k] = coeffs->data.F64[k];
+                }
+            }
+        }
+
+
+        psFree(ALUD);
+        psFree(coeffs);
+        psFree(outPerm);
+    }
+
+    psFree(A);
+    psFree(B);
+    for (psS32 i=0;i<numPoly;i++) {
+        psFree(chebPolys[i]);
+    }
+    psFree(chebPolys);
+
+    return(myPoly);
+}
+
+
 /******************************************************************************
 vectorFitPolynomial1DChebSlow():  This routine will fit a Chebyshev polynomial
@@ -284 +446 @@
 polynomial.
  
-    NOTE: We currently have implemented two algorithms.  This one is
+    NOTE: We currently have implemented three algorithms.  This one is
     non-standard.  It ignores the orthogonal properties of the Chebyshev
     polys, and their known zero values.  Instead, we do build a system of
     linear equations based on minimizing the chi-squared for all data points
     and we then solve those equations.  This method is significantly slower
-    than the other algorithm.  It was explicitly requested that we implement
+    than the other algorithms.  It was explicitly requested that we implement
     this algorithm.
- 
 *****************************************************************************/
 static psPolynomial1D *vectorFitPolynomial1DChebySlow(
@@ -329 +490 @@
         NUM_DATA = x->n;
     }
-    // psTrace    printf("vectorFitPolynomial1DChebySlow(): NUM_DATA is %d\n", NUM_DATA);
 
     psPolynomial1D **chebPolys = p_psCreateChebyshevPolys(NUM_POLY);
@@ -408 +568 @@
         psVector* coeffs = NULL;
 
-        // XXX: Check return codes.
         ALUD = psImageAlloc(NUM_POLY, NUM_POLY, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
-        coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-        for (psS32 k = 0; k < NUM_POLY; k++) {
-            myPoly->coeff[k] = coeffs->data.F64[k];
+        if (ALUD == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+            if (coeffs == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
+                psFree(myPoly);
+                myPoly = NULL;
+            } else {
+                for (psS32 k = 0; k < NUM_POLY; k++) {
+                    myPoly->coeff[k] = coeffs->data.F64[k];
+                }
+            }
         }
 
@@ -437 +608 @@
 polynomial.
  
-    NOTE: We currently have two algorithms.  This is standard method which
+    NOTE: We currently have three algorithms.  This is standard method which
     uses the orthogonal properties of the Chebyshev polys, and their known
-    zero values.  This is significantly faster than the chi-squared approach.
+    zero values.  This is significantly faster than the chi-squared
+    approaches.
  
     NOTE: This function will not work properly if the x-vector does not fully span
@@ -621 +793 @@
     // Build the B and A data structs.
     // XXX EAM : use temp pointers eg vB = B->data.F64 to save redirects
-    // XXX EAM : this function is only valid for data vectors of F64
     for (psS32 k = 0; k < f->n; k++) {
         if ((mask != NULL) && (mask->data.U8[k] && maskValue)) {
@@ -671 +842 @@
         psVector* coeffs = NULL;
 
-        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
-        coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-        for (psS32 k = 0; k < nTerm; k++) {
-            myPoly->coeff[k] = coeffs->data.F64[k];
-        }
+        if (ALUD == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+            if (coeffs == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
+                psFree(myPoly);
+                myPoly = NULL;
+            } else {
+                for (psS32 k = 0; k < nTerm; k++) {
+                    myPoly->coeff[k] = coeffs->data.F64[k];
+                }
+            }
+        }
+
         psFree(ALUD);
         psFree(coeffs);
@@ -771 +954 @@
 
         if (1) {
-            poly = vectorFitPolynomial1DChebySlow(poly, mask, maskValue, f64, fErr64, x64);
+            poly = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64);
         } else {
-            poly = vectorFitPolynomial1DChebyFast(poly, mask, maskValue, f64, fErr64, x64);
+            if (0) {
+                poly = vectorFitPolynomial1DChebySlow(poly, mask, maskValue, f64, fErr64, x64);
+                poly = vectorFitPolynomial1DChebyFast(poly, mask, maskValue, f64, fErr64, x64);
+            }
         }
         if (x == NULL) {
@@ -989 +1175 @@
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL);
     PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL);
-
     PS_ASSERT_VECTOR_NON_NULL(f, NULL);
     PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL);
@@ -1560 +1745 @@
         psVector* coeffs = NULL;
 
-        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
-        coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-
-        // select the appropriate solution entries
-        for (psS32 ix = 0; ix < nXterm; ix++) {
-            for (psS32 iy = 0; iy < nYterm; iy++) {
-                for (psS32 iz = 0; iz < nZterm; iz++) {
-                    psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
-                    myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
-                    myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+        if (ALUD == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+            if (coeffs == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
+                psFree(myPoly);
+                myPoly = NULL;
+            } else {
+                // select the appropriate solution entries
+                for (psS32 ix = 0; ix < nXterm; ix++) {
+                    for (psS32 iy = 0; iy < nYterm; iy++) {
+                        for (psS32 iz = 0; iz < nZterm; iz++) {
+                            psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm;
+                            myPoly->coeff[ix][iy][iz] = coeffs->data.F64[nx];
+                            myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+                        }
+                    }
                 }
             }
@@ -2108 +2303 @@
         psVector* coeffs = NULL;
 
-        // XXX: Check return codes.
         ALUD = psImageAlloc(nTerm, nTerm, PS_TYPE_F64);
         ALUD = psMatrixLUD(ALUD, &outPerm, A);
-        coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
-
-        // select the appropriate solution entries
-        for (psS32 ix = 0; ix < nXterm; ix++) {
-            for (psS32 iy = 0; iy < nYterm; iy++) {
-                for (psS32 iz = 0; iz < nZterm; iz++) {
-                    for (psS32 it = 0; it < nTterm; it++) {
-                        psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
-                        myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
-                        myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+        if (ALUD == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix.  Returning NULL.\n");
+            psFree(myPoly);
+            myPoly = NULL;
+        } else {
+            coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm);
+            if (coeffs == NULL) {
+                psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix.  Returning NULL.\n");
+                psFree(myPoly);
+                myPoly = NULL;
+            } else {
+                // select the appropriate solution entries
+                for (psS32 ix = 0; ix < nXterm; ix++) {
+                    for (psS32 iy = 0; iy < nYterm; iy++) {
+                        for (psS32 iz = 0; iz < nZterm; iz++) {
+                            for (psS32 it = 0; it < nTterm; it++) {
+                                psS32 nx = ix+iy*nXterm+iz*nXterm*nYterm+it*nXterm*nYterm*nZterm;
+                                myPoly->coeff[ix][iy][iz][it] = coeffs->data.F64[nx];
+                                myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+                            }
+                        }
                     }
                 }

trunk/psLib/src/math/psPolynomial.c

@@ -7 +7 @@
 *  polynomials.  It also contains a Gaussian functions.
 *
-*  @version $Revision: 1.136 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-01-23 22:25:31 $
+*  @version $Revision: 1.137 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-01-26 00:31:19 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -194 +194 @@
         }
 
+        if (psTraceGetLevel(__func__) >= 6) {
+            for (psS32 j = 0; j < numPolys; j++) {
+                PS_POLY_PRINT_1D(chebPolys[j]);
+            }
+        }
         return (chebPolys);
     }
@@ -606 +611 @@
 /*****************************************************************************
     This routine must allocate memory for the polynomial structures.
+ 
+    XXX: How do we check for an appropriate value for n?
  *****************************************************************************/
 psPolynomial1D* psPolynomial1DAlloc(
@@ -611 +618 @@
     psPolynomialType type)
 {
-    PS_ASSERT_INT_NONNEGATIVE(n, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(n, NULL);
     psU32 nOrder = n;
     psPolynomial1D *newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D));
@@ -636 +643 @@
                                      psPolynomialType type)
 {
-    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
 
     unsigned int x = 0;
@@ -674 +681 @@
                                      psPolynomialType type)
 {
-    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
 
     unsigned int x = 0;
@@ -723 +730 @@
                                      psPolynomialType type)
 {
-    PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
-    PS_ASSERT_INT_NONNEGATIVE(nT, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nX, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nY, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nZ, NULL);
+    //PS_ASSERT_INT_NONNEGATIVE(nT, NULL);
 
     unsigned int x = 0;

trunk/psLib/src/math/psStats.c

@@ -14 +14 @@
  *      stats->binsize
  *
- *
- *
- *
- *  @version $Revision: 1.160 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-23 22:25:31 $
+ *  @version $Revision: 1.161 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-01-26 00:31:19 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -574 +571 @@
     nValues = p_psVectorNValues(myVector, maskVector, maskVal, stats);
 
-    // XXX: Is a warning message appropriate?  What value should we set the
-    // sampleMedian to?  Should we generate an error?
     if (nValues <= 0) {
         psLogMsg(__func__, PS_LOG_WARN, "WARNING: p_psVectorSampleMedian(): no valid elements in input vector.\n");
@@ -784 +779 @@
 Returns
     NULL
- 
-XXX: Use static vectors.
  *****************************************************************************/
 bool p_psVectorSampleQuartiles(const psVector* myVector,
@@ -1133 +1126 @@
         p_psMemSetPersistent(statsTmp, true);
     } else {
-        // XXX EAM : initialize structure if already allocated
+        // EAM : initialize structure if already allocated
         statsTmp->sampleMean = NAN;
         statsTmp->sampleMedian = NAN;
@@ -1226 +1219 @@
         // Otherwise, use the new results and continue.
         if (isnan(statsTmp->sampleMean) || isnan(statsTmp->sampleStdev)) {
-            // Exit loop.  XXX: Should we throw an error/warning here?
+            // Exit loop.  XXX: throw an error/warning here?
             iter = stats->clipIter;
             rc = -1;
@@ -1487 +1480 @@
         //
         // XXX: This routine should probably be rewritten in a more general fashion
-        // so that the folloiwng checks are not necessary.
+        // so that the following checks are not necessary.
         //
         if (y->data.F64[0] < y->data.F64[1]) {
@@ -1591 +1584 @@
     psTrace(__func__, 5, "---- %s(%f) end ----\n", __func__, tmpFloat);
     return(tmpFloat);
-}
-
-/*****************************************************************************
-XXX: Is there a psLib function for this?
- *****************************************************************************/
-psVector *PsVectorDup(psVector *in)
-{
-    psTrace(__func__, 4, "---- %s() begin ----\n", __func__);
-    PS_ASSERT_VECTOR_NON_NULL(in, NULL);
-    psVector *out = psVectorAlloc(in->n, in->type.type);
-
-    if (in->type.type == PS_TYPE_F32) {
-        for (psS32 i = 0 ; i < in->n ; i++) {
-            out->data.F32[i] = in->data.F32[i];
-        }
-    } else if (in->type.type == PS_TYPE_F64) {
-        for (psS32 i = 0 ; i < in->n ; i++) {
-            out->data.F64[i] = in->data.F64[i];
-        }
-    } else {
-        psError(PS_ERR_UNKNOWN, true, "Unallowable vector type.\n");
-        psTrace(__func__, 4, "---- %s(NULL) end ----\n", __func__);
-        return(NULL);
-    }
-    psTrace(__func__, 4, "---- %s(psVector) end ----\n", __func__);
-    return(out);
 }
 
@@ -1832 +1799 @@
         // ADD: Step 3.
         // Interpolate to the exact 50% position: this is the robust histogram median.
-        // XXX: Check for errors here!
+        // XXX: Check return codes.
         //
         stats->robustMedian = fitQuadraticSearchForYThenReturnX(
@@ -2640 +2607 @@
  
 XXX: Should the default data type be F64?  Since we are buying Opterons...
+ 
+XXX: Use psLib functions instead.
  *****************************************************************************/
 psVector* p_psConvertToF32(psVector* in)

trunk/psLib/test/math/Makefile.am

@@ -5 +5 @@
 
 TESTS = \
-    tst_psFunc00 \
     tst_psFunc01 \
     tst_psFunc08 \
@@ -41 +40 @@
     tst_psSpline1D \
     tst_psRandom \
+    tst_psPolynomial \
     tst_psPolyFit1D \
     tst_psPolyFit2D \
     tst_psPolyFit3D \
-    tst_psPolyFit4D \
-    tst_psPolyFit2DOrd \
-    tst_psPolyFit3DOrd \
-    tst_psPolyFit4DOrd
+    tst_psPolyFit4D
 
 
@@ -85 +82 @@
 tst_psStats09_SOURCES =  tst_psStats09.c
 tst_psRandom_SOURCES =  tst_psRandom.c
+tst_psPolynomial_SOURCES = tst_psPolynomial.c
 tst_psPolyFit1D_SOURCES = tst_psPolyFit1D.c
 tst_psPolyFit2D_SOURCES = tst_psPolyFit2D.c
 tst_psPolyFit3D_SOURCES = tst_psPolyFit3D.c
 tst_psPolyFit4D_SOURCES = tst_psPolyFit4D.c
-tst_psPolyFit2DOrd_SOURCES = tst_psPolyFit2DOrd.c
-tst_psPolyFit3DOrd_SOURCES = tst_psPolyFit3DOrd.c
-tst_psPolyFit4DOrd_SOURCES = tst_psPolyFit4DOrd.c
 
 check_DATA =

trunk/psLib/test/math/tst_psPolyFit1D.c

@@ -14 +14 @@
 #include "psTest.h"
 #define NUM_DATA 35
-#define POLY_ORDER 4
+#define POLY_ORDER 5
 #define A 2.0
 #define B 3.0
@@ -350 +350 @@
     psTraceSetLevel("VectorFitPolynomial1DOrd", 0);
     psTraceSetLevel("VectorFitPolynomial1DCheb", 0);
+    psTraceSetLevel("vectorFitPolynomial1DCheb", 10);
     psTraceSetLevel("vectorFitPolynomial1DChebSlow", 0);
     psTraceSetLevel("vectorFitPolynomial1DChebFast", 0);
     psTraceSetLevel("psVectorFitPolynomial1D", 0);
+    psTraceSetLevel("p_psCreateChebyshevPolys", 0);
 
     printPositiveTestHeader(stdout, "psMinimize functions: 1D Polynomial Fitting Functions", "");
@@ -503 +505 @@
     testStatus &= genericTest(TS00_MASK_U8 | TS00_FERR_F64 | TS00_X_F32 | TS00_F_F64 | TS00_POLY_CHEB | TS00_CLIP_FIT, POLY_ORDER, NUM_DATA, false);
 
+
+    //    testStatus &= genericTest(TS00_MASK_U8 | TS00_FERR_F32 | TS00_X_F32 | TS00_F_F32 | TS00_POLY_CHEB, POLY_ORDER, NUM_DATA, true);
+    //    testStatus &= genericTest(TS00_MASK_U8 | TS00_FERR_F32 | TS00_X_F32 | TS00_F_F32 | TS00_POLY_ORD, POLY_ORDER, NUM_DATA, true);
+
     printFooter(stdout, "psMinimize functions: 1D Polynomial Fitting Functions", "", testStatus);
 }