Changeset 17354

Timestamp:

Apr 6, 2008, 2:44:17 PM (18 years ago)

Author:

eugene

Message:

make consistent with LMM code: Alpha,Beta do not include masked parameters

File:

• branches/eam_branch_20080324/psphot/src/psphotModelWithPSF.c (modified) (16 diffs)

branches/eam_branch_20080324/psphot/src/psphotModelWithPSF.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+# define SAVE_IMAGES 0
 
 bool psphotModelWithPSF_LMM (
@@ -36 +37 @@
 
     // allocate internal arrays (current vs Guess)
-    psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-    psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
-    psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F32);
+    psImage *Alpha = NULL;
+    psVector *Beta = NULL;
+
+    // Alpha & Beta only contain elements to represent the unmasked parameters
+    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
+        psAbort ("programming error: no unmasked parameters to be fit\n");
+    }
+    
+    // allocate internal arrays (current vs Guess)
+    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
 
@@ -79 +87 @@
         // dump some useful info if trace is defined
         if (psTraceGetLevel("psphot") >= 6) {
-            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
-            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
+            p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)");
         }
         if (psTraceGetLevel("psphot") >= 5) {
@@ -116 +125 @@
             beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
             params = psVectorCopy(params, Params, PS_TYPE_F32);
-            lambda *= 0.1;
+            lambda *= 0.25;
 
             // save the new convolved model image
@@ -130 +139 @@
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
-        if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
             psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
         }
+        // set covar values which are not masked
+        psImageInit (covar, 0.0);
+        for (int j = 0, J = 0; j < params->n; j++) {
+            if (paramMask && (paramMask->data.U8[j])) {
+                covar->data.F32[j][j] = 1.0;
+                continue;
+            }
+            for (int k = 0, K = 0; k < params->n; k++) {
+                if (paramMask && (paramMask->data.U8[k])) continue;
+                covar->data.F32[j][k] = Alpha->data.F32[J][K];
+                K++;
+            }
+            J++;
+        }
     }
 
@@ -177 +200 @@
     }
 
-    // generate the model and derivative images for this parameter set
+    // 1 *** generate the model and derivative images for this parameter set
 
     // storage for model derivatives
@@ -223 +246 @@
 
             for (int n = 0; n < params->n; n++) {
-              if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-              psImage *dmodel = pcm->dmodels->data[n];
-              dmodel->data.F32[i][j] = deriv->data.F32[n];
+                if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
+                psImage *dmodel = pcm->dmodels->data[n];
+                dmodel->data.F32[i][j] = deriv->data.F32[n];
             }
         }
@@ -231 +254 @@
     psFree(coord);
     psFree(deriv);
-
 
     // convolve model and dmodel arrays with PSF
     psImageConvolveDirect (pcm->modelConv, pcm->model, psf);
     for (int n = 0; n < pcm->dmodels->n; n++) {
-      if (pcm->dmodels->data[n] == NULL) continue;
-      psImage *dmodel = pcm->dmodels->data[n];
-      psImage *dmodelConv = pcm->dmodelsConv->data[n];
-      psImageConvolveDirect (dmodelConv, dmodel, psf);
+        if (pcm->dmodels->data[n] == NULL) continue;
+        psImage *dmodel = pcm->dmodels->data[n];
+        psImage *dmodelConv = pcm->dmodelsConv->data[n];
+        psImageConvolveDirect (dmodelConv, dmodel, psf);
     }
 
     // XXX TEST : SAVE IMAGES
-# if (1) 
+# if (SAVE_IMAGES) 
     psphotSaveImage (NULL, psf->image, "psf.fits");
     psphotSaveImage (NULL, pcm->model, "model.fits");
@@ -252 +274 @@
 # endif
 
+    // 2 *** accumulate alpha & beta 
+
     // zero alpha and beta for summing below
     psImageInit (alpha, 0.0);
@@ -259 +283 @@
     for (psS32 i = 0; i < source->pixels->numRows; i++) {
         for (psS32 j = 0; j < source->pixels->numCols; j++) {
-          // XXX are we doing the right thing with the mask?
+            // XXX are we doing the right thing with the mask?
             // skip masked points
             if (source->maskObj->data.U8[i][j]) {
@@ -279 +303 @@
             chisq += PS_SQR(delta) * yweight;
 
-            if (isnan(delta))
-              psAbort("nan in delta");
-            if (isnan(chisq))
-              psAbort("nan in chisq");
-
-            for (psS32 n1 = 0; n1 < params->n; n1++) {
-              if ((paramMask != NULL) && (paramMask->data.U8[n1])) {
-                continue;
-              }
-              psImage *dmodel = pcm->dmodelsConv->data[n1];
-              float weight = dmodel->data.F32[i][j] * yweight;
-              for (psS32 n2 = 0; n2 <= n1; n2++) {
-                if ((paramMask != NULL) && (paramMask->data.U8[n2])) {
-                  continue;
-                }
-                dmodel = pcm->dmodelsConv->data[n2];
-                alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j];
-              }
-              beta->data.F32[n1] += weight * delta;
+            if (isnan(delta)) psAbort("nan in delta");
+            if (isnan(chisq)) psAbort("nan in chisq");
+
+            // alpha & beta only contain unmasked elements 
+            for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
+                if ((paramMask != NULL) && (paramMask->data.U8[n1])) continue;
+                psImage *dmodel = pcm->dmodelsConv->data[n1];
+                float weight = dmodel->data.F32[i][j] * yweight;
+                for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
+                    if ((paramMask != NULL) && (paramMask->data.U8[n2])) continue;
+                    dmodel = pcm->dmodelsConv->data[n2];
+                    alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
+                    N2++;
+                }
+                beta->data.F32[N1] += weight * delta;
+                N1++;
             }
         }
@@ -303 +324 @@
 
     // calculate lower-left half of alpha
-    for (psS32 j = 1; j < params->n; j++) {
+    for (psS32 j = 1; j < alpha->numCols; j++) {
         for (psS32 k = 0; k < j; k++) {
             alpha->data.F32[k][j] = alpha->data.F32[j][k];
-        }
-    }
-
-    // fill in pivots if we apply a mask
-    if (paramMask != NULL) {
-        for (psS32 j = 0; j < params->n; j++) {
-            if (paramMask->data.U8[j]) {
-                alpha->data.F32[j][j] = 1;
-                beta->data.F32[j] = 1;
-            }
         }
     }
@@ -345 +356 @@
     pcm->dmodels = psArrayAlloc (params->n);
     for (psS32 n = 0; n < params->n; n++) {
-      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-      pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+        pcm->dmodels->data[n] = NULL;
+        if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
+        pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     }
 
@@ -353 +365 @@
     pcm->dmodelsConv = psArrayAlloc (params->n);
     for (psS32 n = 0; n < params->n; n++) {
-      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-      pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+        pcm->dmodelsConv->data[n] = NULL;
+        if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
+        pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     }
 
@@ -378 +391 @@
  
  while () {
-   GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)
-   dLinear = dLinear(Beta, beta, lambda);
-   chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func)
-   convergence tests...
+ GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)
+ dLinear = dLinear(Beta, beta, lambda);
+ chisq = SetABX (alpha, beta, params, paramMask, x, y, dy, func)
+ convergence tests...
  }
  
@@ -388 +401 @@
  ** GuessABP:
 
-      f_c = sum_i (kern_i * func (x_i; p_o))
-
-df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))]
-
-df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)])
-
-df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)])
-
-- generate image arrays for func, dfunc/dp_j (not masked)
-- convolve each with psf
-- measure delta = f_conv - data
-- etc
+ f_c = sum_i (kern_i * func (x_i; p_o))
+
+ df_c/dp_o = d/dp_o [sum_i (kern_i * func (x_i; p_o))]
+
+ df_c/dp_o = sum_i (d/dp_o [kern_i * func (x_i; p_o)])
+
+ df_c/dp_o = sum_i (kern_i * d/dp_o [func (x_i; p_o)])
+
+ - generate image arrays for func, dfunc/dp_j (not masked)
+ - convolve each with psf
+ - measure delta = f_conv - data
+ - etc
 */