Changeset 25027 for branches/pap/psLib/src/math/psMinimizeLMM.c

Timestamp:

Aug 7, 2009, 4:08:25 PM (17 years ago)

Author:

Paul Price

Message:

Merging trunk (r25026) to get up-to-date on old branch.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psLib/src/math/psMinimizeLMM.c (modified) (4 diffs)

branches/pap/psLib/src/math/psMinimizeLMM.c

@@ -80 +80 @@
     assert (alpha->numCols == alpha->numRows);
 
+# define TESTGJ 0
+# if (TESTGJ)
+    lambda = 0.0;
+#endif
+
     // set new guess values, applying (1+lambda) scaling to pivots
     Beta = psVectorCopy(Beta, beta, PS_TYPE_F32);
@@ -88 +93 @@
 
     // error and clear above if kept?
-    if (false == psMatrixGJSolveF32(Alpha, Beta)) {
+    if (!psMatrixGJSolve(Alpha, Beta)) {
         psTrace ("psLib.math", 4, "singular matrix in Guess ABP\n");
         return(false);
     }
+
+// XXX check that the GJ solver works:
+# if (TESTGJ) 
+    psImage *out = psImageAlloc (alpha->numRows, alpha->numCols, PS_TYPE_F32);
+    for (int oy = 0; oy < out->numRows; oy++) {
+        for (int ox = 0; ox < out->numCols; ox++) {
+            float value = 0;
+            for (int i = 0; i < alpha->numCols; i++) {
+                value += alpha->data.F32[i][ox]*Alpha->data.F32[oy][i];
+            }
+            out->data.F32[oy][ox] = value;
+        }
+    }
+
+    psVector *vect = psVectorAlloc (beta->n, PS_TYPE_F32);
+    for (int oy = 0; oy < vect->n; oy++) {
+        float value = 0;
+        for (int i = 0; i < alpha->numCols; i++) {
+            value += alpha->data.F32[oy][i]*Beta->data.F32[i];
+        }
+        vect->data.F32[oy] = value;
+    }
+    
+    psFree (out);
+    psFree (vect);
+
+# endif
 
     // check for non-finite entries in result
@@ -191 +223 @@
     if (isnan(chisq)) {
         psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN chisq.\n");
+        psVectorInit (delta, NAN);
         retValue = false;
     }
@@ -205 +238 @@
     if (!status) {
         psTrace ("psLib.math", 5, "psMinLM_GuessABP() returned FALSE.\n");
+        psVectorInit (delta, NAN);
         retValue = false;
     }