Changeset 28657

Timestamp:

Jul 11, 2010, 3:21:36 PM (16 years ago)

Author:

eugene

Message:

add maxTol and maxChisqDOF to psMin to allow for better control over fit success and convergence; for SERSIC model, first do a grid search for the index

Location:

branches/eam_branches/ipp-20100621/psphot/src

Files:

• psphotBlendFit.c (modified) (2 diffs)
• psphotChoosePSF.c (modified) (1 diff)
• psphotEllipticalContour.c (modified) (1 diff)
• psphotModelWithPSF.c (modified) (2 diffs)
• psphotPSFConvModel.c (modified) (2 diffs)
• psphotSourceFits.c (modified) (3 diffs)

branches/eam_branches/ipp-20100621/psphot/src/psphotBlendFit.c

@@ -65 +65 @@
     assert (status && fitIter > 0);
 
-    float fitTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_TOL"); // Fit tolerance
-    assert (status && isfinite(fitTol) && fitTol > 0);
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance
+    assert (status && isfinite(fitMinTol) && fitMinTol > 0);
+
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
+    assert (status && isfinite(fitMaxTol) && fitMaxTol > 0);
 
     bool poisson = psMetadataLookupBool(&status, recipe, "POISSON.ERRORS.PHOT.LMM"); // Poisson errors?
     assert (status);
+
+    float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_CHISQ"); // Fit tolerance
 
     float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
@@ -77 +82 @@
     pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
     fitOptions->nIter         = fitIter;
-    fitOptions->tol           = fitTol;
+    fitOptions->minTol        = fitMinTol;
+    fitOptions->maxTol        = fitMaxTol;
+    fitOptions->maxChisqDOF   = maxChisqDOF;
     fitOptions->poissonErrors = poisson;
     fitOptions->weight        = PS_SQR(skySig);

branches/eam_branches/ipp-20100621/psphot/src/psphotChoosePSF.c

@@ -138 +138 @@
         return false;
     }
-    float fitTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
-    if (!status || !isfinite(fitTol) || fitTol <= 0) {
-        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_TOL is not positive");
-        return false;
-    }
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MIN_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
+            return false;
+        }
+    }
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
+        fitMaxTol = 1.0;
+    }
+    float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance
 
     // options which modify the behavior of the model fitting
     options->fitOptions                = pmSourceFitOptionsAlloc();
     options->fitOptions->nIter         = fitIter;
-    options->fitOptions->tol           = fitTol;
+    options->fitOptions->minTol        = fitMinTol;
+    options->fitOptions->maxTol        = fitMaxTol;
+    options->fitOptions->maxChisqDOF   = maxChisqDOF;
     options->fitOptions->poissonErrors = options->poissonErrorsPhotLMM;
     options->fitOptions->weight        = PS_SQR(SKY_SIG);

branches/eam_branches/ipp-20100621/psphot/src/psphotEllipticalContour.c

@@ -82 +82 @@
     params->data.F32[PAR_RMIN]    = Rmin;
 
-    psMinimization *myMin = psMinimizationAlloc (25, 0.001);
+    psMinimization *myMin = psMinimizationAlloc (25, 0.01, 1.00);
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     

branches/eam_branches/ipp-20100621/psphot/src/psphotModelWithPSF.c

@@ -73 +73 @@
 
     // iterate until the tolerance is reached, or give up
-    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
+    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->tol is %f.\n", min->lastDelta, min->tol);
+        psTrace("psphot", 5, "Last delta is %f.  Min->minTol is %f.\n", min->lastDelta, min->minTol);
 
 
@@ -166 +166 @@
     psFree(pcm);
 
-    if (min->iter == min->maxIter) {
-        psTrace("psphot", 3, "---- end (false) ----\n");
-        return(false);
-    }
-
-    psTrace("psphot", 3, "---- end (true) ----\n");
-    return(true);
+    // 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");
+        return(true);
+    }
+    psTrace("psphot", 6, "---- end (false) ----\n");
+    return(false);
 }
 

branches/eam_branches/ipp-20100621/psphot/src/psphotPSFConvModel.c

@@ -4 +4 @@
 // save as static values so they may be set externally
 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
-static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
+static psF32 PM_SOURCE_FIT_MODEL_MIN_TOL = 0.1;
+static psF32 PM_SOURCE_FIT_MODEL_MAX_TOL = 2.0;
 
 // input source has both modelPSF and modelEXT.  on successful exit, we set the
@@ -90 +91 @@
 
     // set up the minimization process
-    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_MIN_TOL, PM_SOURCE_FIT_MODEL_MAX_TOL);
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);

branches/eam_branches/ipp-20100621/psphot/src/psphotSourceFits.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal);
 
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
@@ -460 +461 @@
 
     // for sersic models, get the initial guess more carefully
-    // if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
-    //  psphotGuessSersic ();
-    //  // test and return NULL on failure?
-    // }
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        psphotFitSersicIndex (source, EXT, fitOptions, maskVal);
+    }
 
     // fit EXT (not PSF) model (set/unset the pixel mask)
@@ -473 +473 @@
 }
 
+float indexGuess[] = {0.5, 0.25, 0.125};
+
 // A sersic model is very sensitive to the index.  attempt to find the index first by grid search in just the index
-// bool psphotFitSersic (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
-// 
-//     assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
-// 
-// 
-//     
-//     if (!model->modelGuess(model, source)) {
-//     }
-//     
-// 
-// }
+// for a sersic model, attempt to fit just the index and normalization with a modest number of iterations
+bool psphotFitSersicIndex (pmSource *source, pmModel *model, pmSourceFitOptions *fitOptions, psImageMaskType maskVal) {
+
+    assert (model->type == pmModelClassGetType("PS_MODEL_SERSIC"));
+
+    pmSourceFitOptions options = *fitOptions;
+    
+    // fit EXT (not PSF) model (set/unset the pixel mask)
+    options.mode = PM_SOURCE_FIT_NO_INDEX;
+    options.nIter = 3;
+
+    float xMin, chiSquare[3];
+    int iMin;
+
+    for (int i = 0; i < 3; i++) {
+        model->params->data.F32[PM_PAR_7] = indexGuess[i];
+        model->modelGuess(model, source);
+        pmSourceFitModel (source, model, &options, maskVal);
+        chiSquare[i] = model->chisq;
+        if (i == 0) {
+            xMin = chiSquare[i];
+            iMin = i;
+        } else {
+            if (chiSquare[i] < xMin) {
+                xMin = chiSquare[i];
+                iMin = i;
+            }
+        }
+    }
+
+    model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
+    model->params->data.F32[PM_PAR_7] = indexGuess[iMin];
+    model->modelGuess(model, source);
+
+    return true;
+}