Changeset 29004 for trunk/psModules/src/objects/pmSourceFitSet.c

Timestamp:

Aug 20, 2010, 1:14:11 PM (16 years ago)

Author:

eugene

Message:

Many changes:

• psphot-related issues

1) added a new feature to psLib/src/fft/psImageFFT to allow for a pre-FFTed kernel generated for a specific image dimension which can then be applied to an arbitrary number of images for convolution -- this avoids the kernel FFT for every image convolution.
2) updated the psMinimizeLMM to include 2 levels of tolerance: minTol and maxTol: if the minimization reaches minTol, it stops. if it only reaches maxTol within maxIter, it stops and is considered successful, if it fails to reach maxTol, then it fails. This allows us to accept fits that are actually acceptable (within error) even if they are not as close as ideally possible.
3) add new stat: psfWeightNotPoor (vs psfWeightNotBad) : the first gives the fraction of psf-weighted unmasked pixels considering any mask bits (except the internal 'mark' bit), while the former considers only 'bad' mask bits -- these are written to QF_PSF and QF_PSF_PERFECT in the CMF files.
4) define user-set parameters for max and min valid flux in the linear fit analysis: Note this was the cause of the non-negative fluxes in forced photometry -- the min limit was hard-wired to 0.0.
5) significance image is now constructed as (image + (image/1000)^{2) / sqrt(variance) so that bright sources have well-defined peaks (other wise, they become somewhat flat-topped as image = sqrt(variance), leading to ill-defined peaks).
6) modification of the visualization functions to accept facility and level values akin to psTrace
7) modification of the source fitting APIs to allow thread- and source- independent fit options (iterations, tolerances, etc)
8) use Kron magnitude as test for source size (CR vs EXT) instead of PSF-based aperture
9) set a min systematic error in the aperture mags when used for size classification significance
10) threaded the psf-convolved model (PCM) fitting process (extended source fits)
11) for extended sources, adjust the radius based on the footprint and re-calculate moments for guess; save the psf-based moments for output in psf table
12) for Sersic fitting, choose the best index with a grid search using only a few iterations; iterate fully on the selected value.
13) same for Sersic fitting using the PCM fitting process
14) fixed a bug in which the PSF candidate stars which failed in the psf fitting were not correctly unmarked as being PSF stars
15) define galaxy fit radius based on sky stdev (global, not local)
16) move the PCM code to psModules
17) save and report the raw aperture mag in addition to the curve-of-growth corrected one (PS1_V3)
18) save and report the Kron parameters and higher-order moments
19) some psModule header reorganization for code clarity
20) fixed a bug in which the diff stats were counting 'marked' pixels (outside area of interest) as 'masked'.
21) save the radial profile aperture sizes in the headers
22) better guess for Sersic parameters based on moments (requires setting the functional form so that the scale length is right)}

• ppSub-related:

1) ensure masked pixels are NANed in output diff image
2) add code to flag detections if they have bright positive neighbors (+ output of these in PS1_DV2)
3) define separate penalties for each image (based on their fwhm values) (this requires the penalties to be calculated later in the code).
4) define separate apertures for each image for flux normalization
5) choose aperture based on curve-of-growth (was based on fixed fraction of full aperture flux, and thus noisy)
6) some fine tuning of the penalty factor (this still seems arbitrary, and results are somewhat sensitive to the right value)

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/pmSourceFitSet.c (modified) (5 diffs)

trunk/psModules/src/objects/pmSourceFitSet.c

@@ -22 +22 @@
 #include "pmHDU.h"
 #include "pmFPA.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 "pmGrowthCurve.h"
-#include "pmResiduals.h"
-#include "pmTrend2D.h"
-#include "pmPSF.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 "pmModelClass.h"
+
 #include "pmSourceFitModel.h"
 #include "pmSourceFitSet.h"
 
 // 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_WEIGHT = 1.0;
-static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
+// static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
+// static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
+// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
+// static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
 
 /********************* Source Model Set Functions ***************************/
@@ -429 +436 @@
 bool pmSourceFitSet (pmSource *source,
                      psArray *modelSet,
-                     pmSourceFitMode mode,
+                     pmSourceFitOptions *options,
                      psImageMaskType maskVal)
 {
@@ -478 +485 @@
             // as variance to avoid the bias from systematic errors here we would just use the
             // source sky variance
-            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
-                yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
-            } else {
-                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
-            }
-            nPix++;
-        }
+            if (options->poissonErrors) {
+                yErr->data.F32[nPix] = 1.0 / source->variance->data.F32[i][j];
+            } else {
+                yErr->data.F32[nPix] = 1.0 / options->weight;
+            }
+            nPix++;
+        }
     }
     x->n = nPix;
@@ -490 +497 @@
     yErr->n = nPix;
 
-    // create the FitSet for this thread and set the initial parameter guesses
+// create the FitSet for this thread and set the initial parameter guesses
     pmSourceFitSetData *thisSet = pmSourceFitSetDataSet(modelSet);
 
-    // define param and deriv vectors for complete set of parameters
+// define param and deriv vectors for complete set of parameters
     psVector *params = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
 
-    // set the param and deriv vectors based on the curent values
+// set the param and deriv vectors based on the curent values
     pmSourceFitSetJoin (NULL, params, thisSet);
 
-    // create the minimization constraints
+// create the minimization constraints
     psMinConstraint *constraint = psMinConstraintAlloc();
     constraint->paramMask = psVectorAlloc (thisSet->nParamSet, PS_TYPE_VECTOR_MASK);
     constraint->checkLimits = pmSourceFitSetCheckLimits;
 
-    pmSourceFitSetMasks (constraint, thisSet, mode);
-
-    // force the floating parameters to fall within the contraint ranges
+    pmSourceFitSetMasks (constraint, thisSet, options->mode);
+
+// force the floating parameters to fall within the contraint ranges
     for (int i = 0; i < params->n; i++) {
-        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
-        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
+        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
+        pmSourceFitSetCheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     }
 
     if (psTraceGetLevel("psModules.objects") >= 5) {
-        for (int i = 0; i < params->n; i++) {
-            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
-        }
+        for (int i = 0; i < params->n; i++) {
+            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
+        }
     }
 
     if (nPix <  thisSet->nParamSet + 1) {
-        psTrace (__func__, 4, "insufficient valid pixels\n");
-        psTrace("psModules.objects", 10, "---- %s() end : fail pixels ----\n", __func__);
-        for (int i = 0; i < modelSet->n; i++) {
-            pmModel *model = modelSet->data[i];
-            model->flags |= PM_MODEL_STATUS_BADARGS;
-        }
-        psFree (x);
-        psFree (y);
-        psFree (yErr);
-        psFree (params);
-        psFree(constraint);
-        pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
-        return(false);
-    }
-
-    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
+        psTrace (__func__, 4, "insufficient valid pixels\n");
+        psTrace("psModules.objects", 10, "---- %s() end : fail pixels ----\n", __func__);
+        for (int i = 0; i < modelSet->n; i++) {
+            pmModel *model = modelSet->data[i];
+            model->flags |= PM_MODEL_STATUS_BADARGS;
+        }
+        psFree (x);
+        psFree (y);
+        psFree (yErr);
+        psFree (params);
+        psFree(constraint);
+        pmSourceFitSetDataClear(); // frees thisSet and removes if from the array of fitSets
+        return(false);
+    }
+
+    psMinimization *myMin = psMinimizationAlloc (options->nIter, options->minTol, options->maxTol);
 
     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
@@ -540 +547 @@
     fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSetFunction);
     if (!fitStatus) {
-        psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n);
-    }
-
-    // parameter errors from the covariance matrix
+        psTrace("psModules.objects", 4, "Failed to fit model (%ld components)\n", modelSet->n);
+    }
+
+// parameter errors from the covariance matrix
     psVector *dparams = psVectorAlloc (thisSet->nParamSet, PS_TYPE_F32);
     for (int i = 0; i < dparams->n; i++) {
-        if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
-            continue;
-        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
-    }
-
-    // get the Gauss-Newton distance for fixed model parameters
+        if ((constraint->paramMask != NULL) && constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
+            continue;
+        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
+    }
+
+// get the Gauss-Newton distance for fixed model parameters
     if (constraint->paramMask != NULL) {
-        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
-        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
-        altmask->data.PS_TYPE_VECTOR_MASK_DATA[0] = 1;
-        for (int i = 1; i < dparams->n; i++) {
-            altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1;
-        }
-        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction);
-
-        for (int i = 0; i < dparams->n; i++) {
-            if (!constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
-                continue;
-            // note that delta is the value *subtracted* from the parameter
-            // to get the new guess.  for dparams to represent the direction
-            // of motion, we need to take -delta
-            dparams->data.F32[i] = -delta->data.F32[i];
-        }
-        psFree (delta);
-        psFree (altmask);
+        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
+        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
+        altmask->data.PS_TYPE_VECTOR_MASK_DATA[0] = 1;
+        for (int i = 1; i < dparams->n; i++) {
+            altmask->data.PS_TYPE_VECTOR_MASK_DATA[i] = (constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) ? 0 : 1;
+        }
+        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSetFunction);
+
+        for (int i = 0; i < dparams->n; i++) {
+            if (!constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[i])
+                continue;
+            // note that delta is the value *subtracted* from the parameter
+            // to get the new guess.  for dparams to represent the direction
+            // of motion, we need to take -delta
+            dparams->data.F32[i] = -delta->data.F32[i];
+        }
+        psFree (delta);
+        psFree (altmask);
     }