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

Timestamp:

Mar 6, 2006, 10:57:42 AM (20 years ago)

Author:

gusciora

Message:

This merges in rel10_ifa. However, I changed the matrix solving
routines from GaussJordan to LUD.

File:

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

trunk/psLib/src/math/psMinimizePolyFit.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.8 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-24 23:43:15 $
+ *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-03-06 20:57:42 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -41 +41 @@
 #define PS_VECTOR_GEN_CHEBY_INDEX(VEC, SIZE, TYPE) \
 VEC = psVectorAlloc(SIZE, TYPE); \
+VEC->n = VEC->nalloc; \
 if (TYPE == PS_TYPE_F64) { \
     for (psS32 i = 0 ; i < SIZE ; i++) { \
         VEC->data.F64[i] = ((2.0 / ((psF64) (SIZE - 1))) * ((psF64) i)) - 1.0; \
-        VEC->n++; \
     }\
 } else if (TYPE == PS_TYPE_F32){ \
     for (psS32 i = 0 ; i < SIZE ; i++) { \
         VEC->data.F32[i] = ((2.0 / ((psF32) (SIZE - 1))) * ((psF32) i)) - 1.0; \
-        VEC->n++; \
     }\
 }\
@@ -309 +308 @@
     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.
@@ -370 +369 @@
     psImage *A = psImageAlloc(numPoly, numPoly, PS_TYPE_F64);
     psVector *B = psVectorAlloc(numPoly, PS_TYPE_F64);
+    B->n = B->nalloc;
     for (psS32 i = 0 ; i < numPoly ; i++) {
         for (psS32 j = 0 ; j < numPoly ; j++) {
@@ -375 +375 @@
         }
         B->data.F64[i] = ordPoly->coeff[i];
-        B->n++;
     }
 
@@ -817 +816 @@
         }
         for (psS32 i = 0; i < nTerm; i++) {
+            if (myPoly->mask[i])
+                continue;
             B->data.F64[i] += f->data.F64[k] * xSums->data.F64[i] * wt;
         }
@@ -823 +824 @@
         // we must handle masked orders
         for (psS32 i = 0; i < nTerm; i++) {
+            if (myPoly->mask[i])
+                continue;
             for (psS32 j = 0; j < nTerm; j++) {
+                if (myPoly->mask[j])
+                    continue;
                 A->data.F64[i][j] += xSums->data.F64[i + j] * wt;
             }
@@ -829 +834 @@
     }
 
+    // set masked elements in A,B appropriately
+    for (int i = 0; i < nTerm; i++) {
+        if (!myPoly->mask[i])
+            continue;
+        B->data.F64[i] = 0;
+        for (int j = 0; j < nTerm; j++) {
+            A->data.F64[i][j] = (i == j) ? 1 : 0;
+        }
+    }
+
+
+    // XXX: rel10_ifa used psGaussJordan().  However, this failed tests.  So, I'm using psMatrixLUD().
     if (0) {
         // GaussJordan version
@@ -840 +857 @@
             for (psS32 k = 0; k < nTerm; k++) {
                 myPoly->coeff[k] = B->data.F64[k];
+                myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
             }
         }
@@ -864 +882 @@
                 for (psS32 k = 0; k < nTerm; k++) {
                     myPoly->coeff[k] = coeffs->data.F64[k];
+                    // XXX LUD does not give inverse of A
                 }
             }
@@ -1245 +1264 @@
         for (psS32 n = 0; n < nXterm; n++) {
             for (psS32 m = 0; m < nYterm; m++) {
+                if (myPoly->mask[n][m])
+                    continue;
                 B->data.F64[n+m*nXterm] += f->data.F64[k] * Sums->data.F64[n][m] * wt;
             }
@@ -1251 +1272 @@
         for (psS32 i = 0; i < nXterm; i++) {
             for (psS32 j = 0; j < nYterm; j++) {
+                if (myPoly->mask[i][j])
+                    continue;
                 for (psS32 n = 0; n < nXterm; n++) {
                     for (psS32 m = 0; m < nYterm; m++) {
+                        if (myPoly->mask[n][m])
+                            continue;
                         A->data.F64[i+j*nXterm][n+m*nXterm] += Sums->data.F64[i+n][j+m] * wt;
                     }
+                }
+            }
+        }
+    }
+
+    // set masked elements appropriately
+    for (int i = 0; i < nXterm; i++) {
+        for (int j = 0; j < nYterm; j++) {
+            if (!myPoly->mask[i][j])
+                continue;
+            int nx = i+j*nXterm;
+            B->data.F64[nx] = 0;
+            for (int n = 0; n < nXterm; n++) {
+                for (int m = 0; m < nYterm; m++) {
+                    int ny = n+m*nXterm;
+                    A->data.F64[nx][ny] = (nx == ny) ? 1 : 0;
                 }
             }
@@ -1269 +1310 @@
             for (psS32 m = 0; m < nYterm; m++) {
                 myPoly->coeff[n][m] = B->data.F64[n+m*nXterm];
+                myPoly->coeffErr[n][m] = sqrt(A->data.F64[n+m*nXterm][n+m*nXterm]);
             }
         }
@@ -1723 +1765 @@
     }
 
+    // XXX: rel10_ifa used psGaussJordan().  However, this failed tests.  So, I'm using psMatrixLUD().
     if (0) {
         // does the solution in place
@@ -1775 +1818 @@
                             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]);
+                            // XXX myPoly->coeffErr[ix][iy][iz] = sqrt(A->data.F64[nx][nx]);
+                            // XXX this is wrong: A is not replaced with inverse(A), as for GaussJordan
                         }
                     }
@@ -2278 +2322 @@
 
 
+    // XXX: rel10_ifa used psGaussJordan().  However, this failed tests.  So, I'm using psMatrixLUD().
     if (0) {
         // does the solution in place
@@ -2336 +2381 @@
                                 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]);
+                                // myPoly->coeffErr[ix][iy][iz][it] = sqrt(A->data.F64[nx][nx]);
+                                // XXX this is wrong: LUD does not supply inverse(A)
                             }
                         }