Index: /branches/eam_branch_20080324/psphot/src/psphotModelWithPSF.c
===================================================================
--- /branches/eam_branch_20080324/psphot/src/psphotModelWithPSF.c	(revision 17353)
+++ /branches/eam_branch_20080324/psphot/src/psphotModelWithPSF.c	(revision 17354)
@@ -1,3 +1,4 @@
 # include "psphotInternal.h"
+# define SAVE_IMAGES 0
 
 bool psphotModelWithPSF_LMM (
@@ -36,8 +37,15 @@
 
     // 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,6 +87,7 @@
         // 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,5 +125,5 @@
             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,7 +139,21 @@
     // 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,5 +200,5 @@
     }
 
-    // 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,7 +246,7 @@
 
 	    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,17 +254,16 @@
     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,4 +274,6 @@
 # endif
 
+    // 2 *** accumulate alpha & beta 
+
     // zero alpha and beta for summing below
     psImageInit (alpha, 0.0);
@@ -259,5 +283,5 @@
     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,23 +303,20 @@
 	    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,17 +324,7 @@
 
     // 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,6 +356,7 @@
     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,6 +365,7 @@
     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,8 +391,8 @@
  
  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,16 +401,16 @@
  ** 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
 */