Changeset 26552

Timestamp:

Jan 8, 2010, 4:45:56 PM (16 years ago)

Author:

Paul Price

Message:

Updating to use new API for psMatrixSolveSVD.

File:

• branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c (modified) (7 diffs)

branches/eam_branches/20091201/psModules/src/imcombine/pmSubtractionEquation.c

@@ -543 +543 @@
             for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) {
                 // Contribution to chi^2: a_i^2 P_i
-                if (!isfinite(penalties->data.F32[i])) {
-                    psAbort ("invalid penalty");
-                }
+                psAssert(isfinite(penalties->data.F32[i]), "Invalid penalty");
                 matrix->data.F64[index][index] -= norm * penalties->data.F32[i];
             }
@@ -1360 +1358 @@
 #if 0
         {
-            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
+            solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-5);
             if (!solution) {
                 psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
@@ -1367 +1365 @@
                 return NULL;
             }
+
+#ifdef TESTING
+            {
+                psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
+                psFits *fits = psFitsOpen("solnVector.fits", "w");
+                for (int i = 0; i < numParams + numParams2; i++) {
+                    fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);
+                    image->data.F64[0][i] = solution->data.F64[i];
+                }
+                psFitsWriteImage(fits, NULL, image, 0, NULL);
+                psFree(image);
+                psFitsClose(fits);
+            }
+#endif
 
             int numSpatial = PM_SUBTRACTION_POLYTERMS(kernels->spatialOrder); // Number of spatial variations
@@ -1405 +1417 @@
             }
 
-            #define TOL 1.0e-5
+            #define TOL 1.0e-3
             int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation
             double norm = fabs(solution->data.F64[normIndex]);  // Normalisation
@@ -1429 +1441 @@
 
                 if (!mask1 && !mask2) {
-                    if (fabs(param1) < fabs(param2)) {
-                        psTrace("psModules.imcombine", 7, "Parameter %d: 1 < 2\n", i);
+                    if (fabs(param1) > fabs(param2)) {
+                        psTrace("psModules.imcombine", 7, "Parameter %d: 1 > 2\n", i);
+                        MASK_BASIS_2(i);
+                    } else {
+                        psTrace("psModules.imcombine", 7, "Parameter %d: 2 > 1\n", i);
                         MASK_BASIS_1(i);
-                    } else {
-                        psTrace("psModules.imcombine", 7, "Parameter %d: 2 < 1\n", i);
-                        MASK_BASIS_2(i);
                     }
                 }
             }
+
+            calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
+                                  sumMatrixX, sumVector1, sumVector2);
 
 #ifdef TESTING
@@ -1456 +1471 @@
             }
 #endif
-
-            calculateEquationDual(&sumMatrix, &sumVector, sumMatrix1, sumMatrix2,
-                                  sumMatrixX, sumVector1, sumVector2);
-        }
-#endif
-
-
-        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector);
+        }
+#endif
+
+        solution = psMatrixSolveSVD(solution, sumMatrix, sumVector, 1e-7);
         if (!solution) {
             psError(PS_ERR_UNKNOWN, false, "SVD solution of least-squares equation failed.\n");
@@ -1484 +1495 @@
         {
             psImage *image = psImageAlloc(1, numParams + numParams2, PS_TYPE_F64);
-            psFits *fits = psFitsOpen("solnVector.fits", "w");
+            psFits *fits = psFitsOpen("solnVectorFix.fits", "w");
             for (int i = 0; i < numParams + numParams2; i++) {
                 fprintf(stderr, "Solution %d: %lf\n", i, solution->data.F64[i]);