Changeset 6432

Timestamp:

Feb 15, 2006, 10:10:27 PM (20 years ago)

Author:

magnier

Message:

fits respect coeff masks, fixed coeff errors

File:

• branches/rel10_ifa/psLib/src/math/psMinimizePolyFit.c (modified) (13 diffs)

branches/rel10_ifa/psLib/src/math/psMinimizePolyFit.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-02 21:09:07 $
+ *  @version $Revision: 1.7.6.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-16 08:10:27 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -810 +810 @@
         }
         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;
         }
@@ -816 +818 @@
         // 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;
             }
@@ -822 +828 @@
     }
 
-    if (0) {
+    // 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;
+        }
+    }
+
+
+    if (1) {
         // GaussJordan version
         if (false == psGaussJordan(A, B)) {
@@ -833 +850 @@
             for (psS32 k = 0; k < nTerm; k++) {
                 myPoly->coeff[k] = B->data.F64[k];
+                myPoly->coeffErr[k] = sqrt(A->data.F64[k][k]);
             }
         }
@@ -857 +875 @@
                 for (psS32 k = 0; k < nTerm; k++) {
                     myPoly->coeff[k] = coeffs->data.F64[k];
+                    // XXX LUD does not give inverse of A
                 }
             }
@@ -1236 +1255 @@
         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;
             }
@@ -1242 +1263 @@
         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;
                 }
             }
@@ -1260 +1301 @@
             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]);
             }
         }
@@ -1712 +1754 @@
     }
 
-    if (0) {
+    if (1) {
         // does the solution in place
         // The matrices were overflowing, so I switched to LUD.
@@ -1764 +1806 @@
                             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
                         }
                     }
@@ -2265 +2308 @@
 
 
-    if (0) {
+    if (1) {
         // does the solution in place
         // The GaussJordan version was overflowing, so I'm using LUD.
@@ -2323 +2366 @@
                                 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)
                             }
                         }