Changeset 6193 for trunk/psLib/src/math/psMinimizePolyFit.c

Timestamp:

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

Author:

gusciora

Message:

Several mods mostly to the psPolynomial fitting routines.

File:

• trunk/psLib/src/math/psMinimizePolyFit.c (modified) (12 diffs)

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]);
+                            }
+                        }
                     }
                 }