Changeset 29060 for branches/sc_branches/trunkTest/psModules/src/objects/pmSourceFitSet.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

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

branches/sc_branches/trunkTest/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);
     }