Changeset 28692 for branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c

Timestamp:

Jul 20, 2010, 7:21:53 PM (16 years ago)

Author:

eugene

Message:

update PCM fitting (not quite ready)

File:

• branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c (modified) (13 diffs)

branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -1 @@
-# include "psphotInternal.h"
-# define SAVE_IMAGES 0
-
-bool psphotModelWithPSF_LMM (
+/* @file  pmPCM_MinimizeChisq.c
+ * structures and functions to support PSF-convolved model fitting
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: 1.29 $
+ * @date $Date: 2009-02-16 22:30:50 $
+ * Copyright 2010 Institute for Astronomy, University of Hawaii
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <strings.h>
+#include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmFPAMaskWeight.h"
+
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmSpan.h"
+#include "pmFootprintSpans.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModelFuncs.h"
+#include "pmModel.h"
+#include "pmModelUtils.h"
+#include "pmModelClass.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
+#include "pmSource.h"
+#include "pmSourceFitModel.h"
+#include "pmPCMdata.h"
+
+# define FACILITY "psModules.objects"
+
+bool pmPCM_MinimizeChisq (
     psMinimization *min,
     psImage *covar,
     psVector *params,
-    psMinConstraint *constraint,
     pmSource *source,
-    const psKernel *psf,
-    psMinimizeLMChi2Func func)
+    pmPCMdata *pcm)
 {
-    psTrace("psphot", 3, "---- begin ----\n");
+    psTrace(FACILITY, 3, "---- begin ----\n");
     PS_ASSERT_PTR_NON_NULL(min, false);
     PS_ASSERT_VECTOR_NON_NULL(params, false);
     PS_ASSERT_VECTOR_NON_EMPTY(params, false);
     PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
-    psVector *paramMask = NULL;
-    if (constraint != NULL) {
-        paramMask = constraint->paramMask;
-        if (paramMask != NULL) {
-          PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
-            PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false);
-        }
-    }
-    PS_ASSERT_PTR_NON_NULL(func, false);
     PS_ASSERT_PTR_NON_NULL(source, false);
-
-    psMinimizeLMLimitFunc checkLimits = NULL;
-    if (constraint) {
-        checkLimits = constraint->checkLimits;
-    }
+    PS_ASSERT_PTR_NON_NULL(pcm, false);
+    PS_ASSERT_VECTOR_TYPE(pcm->constraint->paramMask, PS_TYPE_VECTOR_MASK, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(params, pcm->constraint->paramMask, false);
+
+    psVector *paramMask = pcm->constraint->paramMask;
+
+    psMinimizeLMLimitFunc checkLimits = pcm->constraint->checkLimits;
 
     // this function has test values and current values for several things
@@ -54 +86 @@
     psF32 dLinear = 0.0;
 
-    // generate PCM data storage structure
-    pmPCMData *pcm = pmPCMDataAlloc (params, paramMask, source);
-
     // calculate initial alpha and beta, set chisq (min->value)
-    min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, pcm, source, psf, func);
+    min->value = pmPCM_SetABX(alpha, beta, params, paramMask, pcm, source);
     if (isnan(min->value)) {
         min->iter = min->maxIter;
@@ -64 +93 @@
     }
     // dump some useful info if trace is defined
-    if (psTraceGetLevel("psphot") >= 6) {
+    if (psTraceGetLevel(FACILITY) >= 6) {
         p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
         p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
     }
-    if (psTraceGetLevel("psphot") >= 5) {
+    if (psTraceGetLevel(FACILITY) >= 5) {
         p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
     }
 
     // iterate until the tolerance is reached, or give up
-    while ((min->iter < min->maxIter) && ((min->lastDelta > min->minTol) || !isfinite(min->lastDelta))) {
-        psTrace("psphot", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psphot", 5, "Last delta is %f.  Min->minTol is %f.\n", min->lastDelta, min->minTol);
-
+    bool done = (min->iter >= min->maxIter);
+    while (!done) {
+        psTrace("psModules.objects", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+        psTrace("psModules.objects", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
 
         // set a new guess for Alpha, Beta, Params
         if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
             min->iter ++;
+            if (min->iter >=  min->maxIter) break;
             lambda *= 10.0;
             continue;
@@ -86 +116 @@
 
         // dump some useful info if trace is defined
-        if (psTraceGetLevel("psphot") >= 6) {
+        if (psTraceGetLevel(FACILITY) >= 6) {
             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) {
+        if (psTraceGetLevel(FACILITY) >= 5) {
             p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
         }
 
         // calculate Chisq for new guess, update Alpha & Beta
-        Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, pcm, source, psf, func);
+        Chisq = pmPCM_SetABX(Alpha, Beta, Params, paramMask, pcm, source);
         if (isnan(Chisq)) {
             min->iter ++;
+            if (min->iter >=  min->maxIter) break;
             lambda *= 10.0;
             continue;
@@ -109 +140 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        psTrace("psphot", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
-                Chisq, min->lastDelta, rho);
+        psTrace(FACILITY, 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
 
         // dump some useful info if trace is defined
-        if (psTraceGetLevel("psphot") >= 6) {
+        if (psTraceGetLevel(FACILITY) >= 6) {
             p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
             p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
@@ -119 +149 @@
 
         /* if (Chisq < min->value) {  */
-        if (rho > 0.0) {
+        if (rho >= -1e-6) {
             min->lastDelta = (min->value - Chisq) / (source->pixels->numCols*source->pixels->numRows - params->n);
             min->value = Chisq;
@@ -129 +159 @@
             // save the new convolved model image
             psFree (source->modelFlux);
-            source->modelFlux = pmPCMDataSaveImage(pcm);
+            source->modelFlux = pmPCMdataSaveImage(pcm);
         } else {
             lambda *= 10.0;
         }
         min->iter++;
-    }
-    psTrace("psphot", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+
+        done = (min->iter >= min->maxIter);
+        
+        // check for convergence:
+        float chisqDOF = Chisq / pcm->nDOF;
+        if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) {
+            done |= (min->lastDelta < min->minTol);
+        }
+    }
+    psTrace(FACILITY, 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
 
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
         if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
-            psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
+            psTrace (FACILITY, 5, "failure to calculate covariance matrix\n");
         }
         // set covar values which are not masked
@@ -164 +202 @@
     psFree(Beta);
     psFree(Params);
-    psFree(pcm);
 
     // if the last improvement was at least as good as maxTol, accept the fit:
     if (min->lastDelta <= min->maxTol) {
-        psTrace("psphot", 6, "---- end (true) ----\n");
+        psTrace(FACILITY, 6, "---- end (true) ----\n");
         return(true);
     }
-    psTrace("psphot", 6, "---- end (false) ----\n");
+    psTrace(FACILITY, 6, "---- end (false) ----\n");
     return(false);
 }
 
-psF32 psphotModelWithPSF_SetABX(
+psF32 pmPCM_SetABX(
     psImage  *alpha,
     psVector *beta,
     const psVector *params,
     const psVector *paramMask,
-    pmPCMData *pcm,
-    const pmSource *source,
-    const psKernel *psf,
-    psMinimizeLMChi2Func func)
+    pmPCMdata *pcm,
+    const pmSource *source)
 {
     // XXX: Check vector sizes.
@@ -208 +243 @@
     psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
 
-    psImageInit (pcm->model, 0.0);
+    psImageInit (pcm->modelFlux, 0.0);
     for (int n = 0; n < params->n; n++) {
-        if (!pcm->dmodels->data[n]) continue;
-        psImageInit (pcm->dmodels->data[n], 0.0);
+        if (!pcm->dmodelsFlux->data[n]) continue;
+        psImageInit (pcm->dmodelsFlux->data[n], 0.0);
     }
 
@@ -243 +278 @@
             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
 
-            pcm->model->data.F32[i][j] = func (deriv, params, coord);
+            pcm->modelFlux->data.F32[i][j] = pcm->modelConv->modelFunc (deriv, params, coord);
 
             for (int n = 0; n < params->n; n++) {
                 if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
-                psImage *dmodel = pcm->dmodels->data[n];
+                psImage *dmodel = pcm->dmodelsFlux->data[n];
                 dmodel->data.F32[i][j] = deriv->data.F32[n];
             }
@@ -256 +291 @@
 
     // 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);
+    psImageConvolveDirect (pcm->modelConvFlux, pcm->modelFlux, pcm->psf);
+    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
+        if (pcm->dmodelsFlux->data[n] == NULL) continue;
+        psImage *dmodel = pcm->dmodelsFlux->data[n];
+        psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
+        psImageConvolveDirect (dmodelConv, dmodel, pcm->psf);
     }
 
     // XXX TEST : SAVE IMAGES
 # if (SAVE_IMAGES)
-    psphotSaveImage (NULL, psf->image, "psf.fits");
-    psphotSaveImage (NULL, pcm->model, "model.fits");
-    psphotSaveImage (NULL, pcm->modelConv, "modelConv.fits");
+    psphotSaveImage (NULL, pcm->psf->image, "psf.fits");
+    psphotSaveImage (NULL, pcm->modelFlux, "model.fits");
+    psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits");
     psphotSaveImage (NULL, source->pixels, "obj.fits");
     psphotSaveImage (NULL, source->maskObj, "mask.fits");
@@ -297 +332 @@
             }
 
-            float ymodel  = pcm->modelConv->data.F32[i][j];
+            float ymodel  = pcm->modelConvFlux->data.F32[i][j];
             float yweight = 1.0 / source->variance->data.F32[i][j];
             float delta = ymodel - source->pixels->data.F32[i][j];
@@ -309 +344 @@
             for (int n1 = 0, N1 = 0; n1 < params->n; n1++) {
                 if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n1])) continue;
-                psImage *dmodel = pcm->dmodelsConv->data[n1];
+                psImage *dmodel = pcm->dmodelsConvFlux->data[n1];
                 float weight = dmodel->data.F32[i][j] * yweight;
                 for (int n2 = 0, N2 = 0; n2 <= n1; n2++) {
                     if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n2])) continue;
-                    dmodel = pcm->dmodelsConv->data[n2];
+                    dmodel = pcm->dmodelsConvFlux->data[n2];
                     alpha->data.F32[N1][N2] += weight * dmodel->data.F32[i][j];
                     N2++;