Changeset 30784 for branches/eam_branches/ipp-20110213/psModules/src/objects

Timestamp:

Mar 3, 2011, 3:11:09 PM (15 years ago)

Author:

eugene

Message:

fix the calculation of chisq values for constant errors (re-calculated, do not attempt to correct raw value)

Location:

branches/eam_branches/ipp-20110213/psModules/src/objects

Files:

• pmModel.c (modified) (1 diff)
• pmModel.h (modified) (1 diff)
• pmPCMdata.c (modified) (2 diffs)
• pmPCMdata.h (modified) (1 diff)
• pmSourceFitModel.c (modified) (1 diff)
• pmSourceFitPCM.c (modified) (2 diffs)
• pmSourceFitSet.c (modified) (6 diffs)
• pmSourceOutputs.c (modified) (1 diff)
• pmSourcePhotometry.c (modified) (3 diffs)
• pmSourcePhotometry.h (modified) (1 diff)

branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.c

@@ -67 +67 @@
     tmp->chisqNorm = NAN;
     tmp->nDOF  = 0;
+    tmp->nPar  = 0;
     tmp->nPix  = 0;
     tmp->nIter = 0;

branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.h

@@ -39 +39 @@
     float magErr;                       ///< integrated model magnitude error
     int nPix;                           ///< number of pixels used for fit
-    int nDOF;                           ///< number of degrees of freedom
+    int nPar;                           ///< number of parameters in fit
+    int nDOF;                           ///< number of degrees of freedom (nDOF = nPix - nPar)
     int nIter;                          ///< number of iterations to reach min
     pmModelStatus flags;                ///< model status flags

branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c

@@ -254 +254 @@
 
     pcm->nPix = nPix;
-    pcm->nDOF = nPix - nParams - 1;
+    pcm->nPar = nParams;
+    pcm->nDOF = nPix - nParams;
 
     return pcm;
@@ -341 +342 @@
         return false;
     }
+    pcm->nPar = nParams;
     pcm->nDOF = pcm->nPix - nParams;
 

branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.h

@@ -33 +33 @@
     psMinConstraint *constraint;
     int nPix;
+    int nPar;
     int nDOF;
 } pmPCMdata;

branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c

@@ -236 +236 @@
         model->covar = psMemIncrRefCounter(covar);
     }
+    model->nIter = myMin->iter;
+    model->nPar = nParams;
+
     psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
 
     // save the resulting chisq, nDOF, nIter
+    // NOTE: if (!options->poissonErrors) chisq will be wrong : recalculate
     if (options->poissonErrors) {
         model->chisq = myMin->value;
         model->nPix  = y->n;
-        model->nDOF  = y->n - nParams;
+        model->nDOF  = y->n - model->nPar;
         model->chisqNorm = model->chisq / model->nDOF;
     } else {
-        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, nParams);
-    }
-    model->nIter = myMin->iter;
+        pmSourceChisqUnsubtracted (source, model, maskVal);
+    }
 
     // set the model success or failure status

branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c

@@ -84 +84 @@
         }
     }
+    pcm->modelConv->nIter = myMin->iter;
+    pcm->modelConv->nPar = pcm->nPar;
 
     // save the resulting chisq, nDOF, nIter
@@ -92 +94 @@
         pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
     } else {
-        pmSourceChisq (pcm->modelConv, source->pixels, source->maskObj, source->variance, maskVal, fitOptions->covarFactor, pcm->nPix - pcm->nDOF - 1);
+        // xxx this is wrong because it does not convolve with the psf
+        pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal);
     }
-    pcm->modelConv->nIter = myMin->iter;
 
     // set the model success or failure status

branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c

@@ -335 +335 @@
         psTrace ("psModules.objects", 4, " src %d", i);
 
+        model->nIter = myMin->iter;
+        // model->nPar is set by pmSourceFitSetMasks
+
         // save the resulting chisq, nDOF, nIter
         // these are not unique for any one source
@@ -340 +343 @@
             model->chisq = myMin->value;
             model->nPix  = nPix;
-            model->nDOF  = nPix - model->params->n;
+            model->nDOF  = nPix - model->nPar;
             model->chisqNorm = model->chisq / model->nDOF;
         } else {
-            pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, model->params->n);
+            pmSourceChisqUnsubtracted (source, model, maskVal);
         }
-        model->nIter = myMin->iter;
 
         // set the model success or failure status
@@ -399 +401 @@
     for (int i = 0; i < set->paramSet->n; i++) {
         psVector *paramOne = set->paramSet->data[i];
+        pmModel  *modelOne = set->modelSet->data[i];
 
         switch (mode) {
@@ -406 +409 @@
                 if (j == PM_PAR_I0) continue;
                 constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1;
+                modelOne->nPar = 1;
             }
             break;
@@ -415 +419 @@
                 if (j == PM_PAR_I0) continue;
                 constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1;
+                modelOne->nPar = 3;
             }
             break;
@@ -420 +425 @@
             // EXT model fits all params (except sky)
             constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + PM_PAR_SKY] = 1;
+            modelOne->nPar = paramOne->n - 1;
             break;
           default:

branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c

@@ -72 +72 @@
     }
 
-    *nImageOverlap = 1;
+    *nImageOverlap = psMetadataLookupS32 (&status2, header, "NINPUTS");
     return true;
 

branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c

@@ -741 +741 @@
 # endif
 
-// determine chisq, etc for linear normalization-only fit
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor, int nParams)
+// determine chisq, nPix, nDOF, chisqNorm : model->nPar must be set
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal)
 {
     PS_ASSERT_PTR_NON_NULL(model, false);
@@ -757 +757 @@
             if (variance->data.F32[j][i] <= 0)
                 continue;
-            // dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
             dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
             Npix ++;
         }
     }
-
     model->nPix = Npix;
-    model->nDOF = Npix - nParams - 1;
+    model->nDOF = Npix - model->nPar;
     model->chisq = dC;
     model->chisqNorm = dC / model->nDOF;
@@ -771 +769 @@
 }
 
+
+// return source aperture magnitude
+bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal)
+{
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(model, false);
+
+    float dC = 0.0;
+    int Npix = 0;
+
+    // the model function returns the source flux at a position
+    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+    psVector *params = model->params;
+    psImage  *image = source->pixels;
+    psImage  *mask = source->maskObj;
+    psImage  *variance = source->variance;
+
+    int dX = image->col0;
+    int dY = image->row0;
+
+    for (int iy = 0; iy < image->numRows; iy++) {
+        for (int ix = 0; ix < image->numCols; ix++) {
+
+            // skip pixels which are masked
+            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
+
+            if (variance->data.F32[iy][ix] <= 0) continue;
+
+            coord->data.F32[0] = (psF32) ix + dX + 0.5;
+            coord->data.F32[1] = (psF32) iy + dY + 0.5;
+
+            // for the full model, add all points
+            float value = model->modelFunc (NULL, params, coord);
+
+            // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n", 
+            // ix, iy, coord->data.F32[0], coord->data.F32[1], image->data.F32[iy][ix], value, dC);
+
+            dC += PS_SQR (image->data.F32[iy][ix] - value) / variance->data.F32[iy][ix];
+            Npix ++;
+        }
+    }
+    model->nPix = Npix;
+    model->nDOF = Npix - model->nPar;
+    model->chisq = dC;
+    model->chisqNorm = dC / model->nDOF;
+
+    psFree (coord);
+    return (true);
+}
 
 double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)

branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h

@@ -69 +69 @@
 bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius);
 
-bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor, int nParams);
+bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal);
+bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal);
 
 bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);