Index: /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c
===================================================================
--- /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c	(revision 33070)
+++ /branches/eam_branches/ipp-20111122/psphot/src/psphotFitSourcesLinear.c	(revision 33071)
@@ -24,4 +24,11 @@
     assert (recipe);
 
+    pmSourceFitVarMode fitVarMode = psphotGetFitVarMode (recipe);
+    if (!fitVarMode) {
+	psError (PSPHOT_ERR_CONFIG, true, "failed to get LINEAR_FIT_VARIANCE_MODE");
+	return false;
+    }
+    pmSourceFitVarMode fitVarModePass1 = (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) ? PM_SOURCE_PHOTFIT_CONST : fitVarMode;
+
     int num = psphotFileruleCount(config, filerule);
 
@@ -50,8 +57,25 @@
         psAssert (psf, "missing psf?");
 
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarModePass1)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
             return false;
         }
+
+	// the MODEL_VAR weighting scheme requires knowledge of the model fluxes to generate the variance
+	// after we have determined the initial set of fits, then we can generate the variance image and 
+	// re-run the fit against that variance.
+	if (fitVarMode == PM_SOURCE_PHOTFIT_MODEL_VAR) {
+	    // generate the model variance image
+	    if (!psphotGenerateModelVariance (recipe, readout, sources, psf, final, fitVarMode)) {
+		psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
+		return false;
+	    }
+
+	    // rerun fit with correct fitVarMode
+	    if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final, fitVarMode)) {
+		psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
+		return false;
+	    }
+	}
 
 	psphotVisualShowResidualImage (readout, (num > 0)); 
@@ -62,5 +86,31 @@
 }
 
-bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) {
+pmSourceFitVarMode psphotGetFitVarMode (psMetadata *recipe) {
+
+    bool status = false;
+
+    char *fitVarModeString = psMetadataLookupStr(&status, recipe, "LINEAR_FIT_VARIANCE_MODE");
+    if (!status) {
+	bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+	if (!status) {
+	    psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS");
+	}
+	pmsourceFitVarMode fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
+	return fitVarMode;
+    } 
+    if (!strcasecmp(fitVarModeString, "CONSTANT") || !strcasecmp(fitVarModeString, "CONST")) {
+	return PM_SOURCE_PHOTFIT_CONST;
+    }
+    if (!strcasecmp(fitVarModeString, "IMAGE") || !strcasecmp(fitVarModeString, "IMAGE_VAR")) {
+	return PM_SOURCE_PHOTFIT_IMAGE_VAR;
+    }
+    if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
+	return PM_SOURCE_PHOTFIT_MODEL_VAR;
+    }
+    psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString);
+    return PM_SOURCE_PHOTFIT_MODEL_NONE;
+}
+
+bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final, pmSourceFitVarMode fitVarMode) {
 
     bool status;
@@ -99,30 +149,4 @@
     if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
-    pmSourceFitVarMode fitVarMode = PM_SOURCE_PHOTFIT_CONST;
-    char *fitVarModeString = psMetadataLookupStr(&status, recipe, "LINEAR_FIT_VARIANCE_MODE");
-    if (!status) {
-	bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
-	if (!status) {
-	    psAbort("You must provide a value for LINEAR_FIT_VARIANCE_MODE or CONSTANT_PHOTOMETRIC_WEIGHTS");
-	}
-	fitVarMode = CONSTANT_PHOTOMETRIC_WEIGHTS ? PM_SOURCE_PHOTFIT_CONST : PM_SOURCE_PHOTFIT_IMAGE_VAR;
-    } else {
-	if (!strcasecmp(fitVarModeString, "CONSTANT") || !strcasecmp(fitVarModeString, "CONST")) {
-	    fitVarMode = PM_SOURCE_PHOTFIT_CONST;
-	    goto gotit;
-	}
-	if (!strcasecmp(fitVarModeString, "IMAGE") || !strcasecmp(fitVarModeString, "IMAGE_VAR")) {
-	    fitVarMode = PM_SOURCE_PHOTFIT_IMAGE_VAR;
-	    goto gotit;
-	}
-	if (!strcasecmp(fitVarModeString, "MODEL") || !strcasecmp(fitVarModeString, "MODEL_VAR")) {
-	    fitVarMode = PM_SOURCE_PHOTFIT_MODEL_VAR;
-	    goto gotit;
-	}
-	psError (PSPHOT_ERR_CONFIG, false, "Invalid value for LINEAR_FIT_VARIANCE_MODE (%s)", fitVarModeString);
-	return false;
-    }
-
-gotit:
     int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
     if (!status) {
@@ -259,4 +283,12 @@
     psSparseBorder *border = psSparseBorderAlloc (sparse, nBorder);
 
+    // if fitVarMode is MODEL_VAR, then we need to generate the model image variance
+    // XXX we have two possibilities here: 
+
+    // 1) do 2 passes, where in the first case we use the CONST weighting, and in the second
+    // use the fitted model values to define the model
+
+    // 2) do a single pass, and use the model guess to define the model variance (but do I trust the Model Guess?)
+
     // fill out the sparse matrix elements and border elements (B)
     // SRCi is the current source of interest
@@ -266,10 +298,10 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
+        f = pmSourceModelDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
         psSparseMatrixElement (sparse, i, i, f);
 
-        // the formal error depends on the weighting scheme
-        if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
+        // if we have used CONSTANT errors, then we need to calculate the value of the parameter error
+        if (fitVarMode != PM_SOURCE_PHOTFIT_IMAGE_VAR) {
+            float var = pmSourceModelDotModel (SRCi, SRCi, PM_SOURCE_PHOTFIT_IMAGE_VAR, covarFactor, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -279,5 +311,5 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
+        f = pmSourceDataDotModel (SRCi, SRCi, fitVarMode, covarFactor, maskVal);
         psSparseVectorElement (sparse, i, f);
 
@@ -285,11 +317,11 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
+            f = pmSourceModelWeight (SRCi, 1, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
+            f = pmSourceModelWeight (SRCi, 2, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 2, f);
 
           case 0:
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
+            f = pmSourceModelWeight (SRCi, 0, fitVarMode, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -311,5 +343,5 @@
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
+            f = pmSourceModelDotModel (SRCi, SRCj, fitVarMode, covarFactor, maskVal);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -507,2 +539,83 @@
     return true;
 }
+
+bool psphotGenerateModelVariance (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    // create a model variance image
+    psImage *modelVar = psImageCopy (NULL, readout->variance, PS_TYPE_F32);
+
+    // find the pmFPAfile for the readout (background model is saved on PSPHOT.BACKMDL regardless of 'filename')
+    pmFPAfile *backModellFile = pmFPAfileSelectSingle(config->files, "PSPHOT.BACKMDL", index); // File of interest
+    assert (backModelFile);
+
+    // find the background model readout from the correct location
+    pmReadout *backModel = READOUT_OR_INTERNAL(view, backModelFile);
+    psAssert (backModel, "this must exist");
+
+    // find the binning information
+    psImageBinning *backBinning = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.BACKGROUND.BINNING");
+    assert (backBinning);
+
+    // linear interpolation to full-scale
+    if (!psImageUnbin (modelVar, backModel->image, backBinning)) {
+        psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for unbinning");
+        return false;
+    }
+
+    // XXX for a test:
+    psphotSaveImage (NULL, modelVar, "model.bck.fits");
+
+    // insert all of the source models
+    for (int i = 0; i < sources->n; i++) {
+
+	// source of interest
+	pmSource *source = sources->data[i];
+
+	// XXX which sources need to be injected?
+	// skip non-astronomical objects (very likely defects)
+	if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
+	if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+      
+	// do not include CRs in the full ensemble fit
+	if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+	
+	// do not include MOMENTS_FAILURES in the fit
+	if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+
+	// XXX need to identify the region appropriate for the source
+
+	// define the source->modelVar pixels (view on modelVar image)
+	psAssert (!source->modelVar, "programming error : modelVar should be NULL here");
+	source->modelVar = psImageSubset(readout->image, region);
+
+	// add the source model to the model variance image
+	pmSourceAdd (source, PM_MODEL_OP_XXX, maskVal);
+    }
+
+    // XXX add the readnoise to the image?  This is a bit problematic because I don't have the
+    // per-cell information.  Alternatively, I could use the median of the variance image or
+    // something equivalent to a background model to define the base level of the variance
+
+    // XXX for a test:
+    psphotSaveImage (NULL, modelVar, "model.var.fits");
+
+    return true;
+}
+
+bool psphotFreeModelVariance (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+
+    // find the binning information
+    psImage *modelVar = psMetadataLookupPtr(&status, model->analysis, "PSPHOT.MODEL.VAR");
+    assert (modelVar);
+
+    // clear modelVar pointers for all of the source models
+    for (int i = 0; i < sources->n; i++) {
+
+	// source of interest
+	pmSource *source = sources->data[i];
+	psFree (source->modelVar);
+    }
+    psFree (modelVar);
+
+    return true;
+}