Changeset 25990

Timestamp:

Nov 1, 2009, 4:04:00 PM (17 years ago)

Author:

eugene

Message:

covariance needed when calculating chisq

File:

• trunk/psphot/src/psphotFitSourcesLinear.c (modified) (8 diffs)

trunk/psphot/src/psphotFitSourcesLinear.c

@@ -59 +59 @@
     }
 
+    // XXX test: choose a larger-than expected radius:
+    float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
+
+    // XXX do not apply covarFactor for the moment...
+    // covarFactor = 1.0;
+
     // select the sources which will be used for the fitting analysis
     for (int i = 0; i < sources->n; i++) {
@@ -120 +127 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
         psSparseMatrixElement (sparse, i, i, f);
 
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -133 +140 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
         psSparseVectorElement (sparse, i, f);
 
@@ -139 +146 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 2, f);
 
           case 0:
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -165 +172 @@
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -219 +226 @@
 
     // measure chisq for each source
-    for (int i = 0; final && (i < fitSources->n); i++) {
+    // for (int i = 0; final && (i < fitSources->n); i++) {
+    for (int i = 0; i < fitSources->n; i++) {
         pmSource *source = fitSources->data[i];
         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
         pmModel *model = pmSourceGetModel (NULL, source);
-        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
+        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor);
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -238 +246 @@
 
     psphotVisualShowResidualImage (readout);
+    psphotVisualPlotChisq (sources);
     // psphotVisualShowFlags (sources);
 
@@ -244 +253 @@
 
 // XXX do we need this?
+// XXX disallow the simultaneous sky fit and remove this code...
 
 // Calculate the weight terms for the sky fit component of the matrix.  This function operates