Changeset 29598 for trunk/psModules/src/imcombine/pmSubtractionEquation.c

Timestamp:

Oct 28, 2010, 5:50:53 PM (16 years ago)

Author:

eugene

Message:

adjust the moments penalty scaling by a factor to make the recipe penalties in the vicinity of 1.0

File:

• trunk/psModules/src/imcombine/pmSubtractionEquation.c (modified) (4 diffs)

trunk/psModules/src/imcombine/pmSubtractionEquation.c

@@ -25 +25 @@
 # define PENALTY false
 # define MOMENTS (!PENALTY)
+# define MOMENTS_PENALTY_SCALE 2e-5 // XXX this value is not completely arbitrary, but I don't understand why it needs to be this value...
 
 //////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -325 +326 @@
                     // add in second moments
                     if (MOMENTS) {
-                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA;
-                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA;
-
-                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA;
-                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA;
-
-                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB;
-                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB;
-
-                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB;
-                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB;
-
-                        // XXX this uses the single moments, first try
-                        // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
-                        // matrix->data.F64[iIndex][jIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
-                        // 
-                        // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MxxI1->data.F32[i]*stamp->MxxI1->data.F32[j];
-                        // matrix->data.F64[jIndex][iIndex] += kernels->penalty * stamp->MyyI1->data.F32[i]*stamp->MyyI1->data.F32[j];
-                        // 
-                        // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
-                        // matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
-                        // 
-                        // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MxxI2->data.F32[i]*stamp->MxxI2->data.F32[j];
-                        // matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * stamp->MyyI2->data.F32[i]*stamp->MyyI2->data.F32[j];
+                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
+                        matrix->data.F64[iIndex][jIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
+
+                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MxxAA * MOMENTS_PENALTY_SCALE;
+                        matrix->data.F64[jIndex][iIndex] += kernels->penalty * MyyAA * MOMENTS_PENALTY_SCALE;
+
+                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
+                        matrix->data.F64[iIndex + numParams][jIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
+
+                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MxxBB * MOMENTS_PENALTY_SCALE;
+                        matrix->data.F64[jIndex + numParams][iIndex + numParams] += kernels->penalty * MyyBB * MOMENTS_PENALTY_SCALE;
                     }
                     matrix->data.F64[iIndex][jIndex + numParams] = ab;
@@ -465 +453 @@
             // add in second moments
             if (MOMENTS) {
-                vector->data.F64[iIndex]             -= kernels->penalty * MxxAI1;
-                vector->data.F64[iIndex]             -= kernels->penalty * MyyAI1;
-
-                vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2;
-                vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2;
+                vector->data.F64[iIndex]             -= kernels->penalty * MxxAI1 * MOMENTS_PENALTY_SCALE;
+                vector->data.F64[iIndex]             -= kernels->penalty * MyyAI1 * MOMENTS_PENALTY_SCALE;
+
+                vector->data.F64[iIndex + numParams] -= kernels->penalty * MxxBI2 * MOMENTS_PENALTY_SCALE;
+                vector->data.F64[iIndex + numParams] -= kernels->penalty * MyyBI2 * MOMENTS_PENALTY_SCALE;
             }
         }
@@ -793 +781 @@
         }
     }
+
     stamp->norm = normI2 / normI1;
     stamp->normI1 = normI1;