Index: trunk/psModules/src/objects/pmSourceFitSet.c =================================================================== --- trunk/psModules/src/objects/pmSourceFitSet.c (revision 27903) +++ trunk/psModules/src/objects/pmSourceFitSet.c (revision 29004) @@ -22,23 +22,30 @@ #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,5 +436,5 @@ bool pmSourceFitSet (pmSource *source, psArray *modelSet, - pmSourceFitMode mode, + pmSourceFitOptions *options, psImageMaskType maskVal) { @@ -478,11 +485,11 @@ // 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,49 +497,49 @@ 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,35 +547,35 @@ 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); }