Index: /branches/rel10_ifa/psModules/src/objects/Makefile.am =================================================================== --- /branches/rel10_ifa/psModules/src/objects/Makefile.am (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/Makefile.am (revision 6848) @@ -12,5 +12,4 @@ pmSourceContour.c \ pmSourceFitModel.c \ - pmSourceFitSet.c \ pmSourcePhotometry.c \ pmSourceIO.c \ @@ -42,5 +41,4 @@ pmSourceContour.h \ pmSourceFitModel.h \ - pmSourceFitSet.h \ pmSourcePhotometry.h \ pmSourceIO.h \ Index: /branches/rel10_ifa/psModules/src/objects/pmPSF.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSF.c (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmPSF.c (revision 6848) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-10 20:21:36 $ + * @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -80,5 +80,5 @@ Object Normalization *****************************************************************************/ -pmPSF *pmPSFAlloc (pmModelType type) +pmPSF *pmPSFAlloc (pmModelType type, bool poissonErrors) { int Nparams; @@ -92,4 +92,5 @@ psf->skyBias = 0.0; psf->skySat = 0.0; + psf->poissonErrors = poissonErrors; // the ApTrend components are (x, y, r2rflux, flux) @@ -97,5 +98,5 @@ pmPSF_MaskApTrend (psf, PM_PSF_SKYBIAS); - if (PM_PSF_POISSON_WEIGHTS) { + if (psf->poissonErrors) { psf->ChiTrend = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); } else { @@ -319,4 +320,6 @@ psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NAME", PS_DATA_STRING, "PSF model name", modelName); + psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_POISSON_ERRORS", PS_DATA_BOOL, "Poisson errors for fits", psf->poissonErrors); + int nPar = pmModelParameterCount (psf->type) ; psMetadataAdd (metadata, PS_LIST_TAIL, "PSF_MODEL_NPAR", PS_DATA_S32, "PSF model parameter count", nPar); @@ -347,5 +350,9 @@ pmModelType type = pmModelSetType (modelName); - pmPSF *psf = pmPSFAlloc (type); + bool poissonErrors = psMetadataLookupPtr (&status, metadata, "PSF_POISSON_ERRORS"); + if (!status) + poissonErrors = true; + + pmPSF *psf = pmPSFAlloc (type, poissonErrors); int nPar = psMetadataLookupS32 (&status, metadata, "PSF_MODEL_NPAR"); Index: /branches/rel10_ifa/psModules/src/objects/pmPSF.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSF.h (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmPSF.h (revision 6848) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1.22.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-10 20:21:36 $ + * @version $Revision: 1.1.22.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -15,6 +15,4 @@ # ifndef PM_PSF_H # define PM_PSF_H - -# define PM_PSF_POISSON_WEIGHTS 1 /** pmPSF data structure @@ -44,4 +42,5 @@ int nPSFstars; ///< number of stars used to measure PSF int nApResid; ///< number of stars used to measure ApResid + bool poissonErrors; } pmPSF; @@ -65,5 +64,6 @@ */ pmPSF *pmPSFAlloc( - pmModelType type ///< Add comment + pmModelType type, // type of model for PSF + bool poissonErrors ///< use poissonian errors or not? ); Index: /branches/rel10_ifa/psModules/src/objects/pmPSFtry.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSFtry.c (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmPSFtry.c (revision 6848) @@ -5,6 +5,6 @@ * @author EAM, IfA * - * @version $Revision: 1.4.4.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-10 20:21:36 $ + * @version $Revision: 1.4.4.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -50,5 +50,5 @@ // allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName) +pmPSFtry *pmPSFtryAlloc (psArray *sources, char *modelName, bool poissonErrors) { @@ -59,5 +59,5 @@ // XXX probably need to increment ref counter type = pmModelSetType (modelName); - test->psf = pmPSFAlloc (type); + test->psf = pmPSFAlloc (type, poissonErrors); test->sources = psMemIncrRefCounter(sources); test->modelEXT = psArrayAlloc (sources->n); @@ -88,5 +88,5 @@ // mask values indicate the reason the source was rejected: -pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS) +pmPSFtry *pmPSFtryModel (psArray *sources, char *modelName, float RADIUS, bool poissonErrors) { bool status; @@ -98,5 +98,5 @@ int Npsf = 0; - pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName); + pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName, poissonErrors); // stage 1: fit an independent model (freeModel) to all sources @@ -113,5 +113,5 @@ psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED); - status = pmSourceFitModel (source, model, false); + status = pmSourceFitModel (source, model, PM_SOURCE_FIT_EXT); psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PM_SOURCE_MASK_MARKED); @@ -147,5 +147,5 @@ psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PM_SOURCE_MASK_MARKED); - status = pmSourceFitModel (source, modelPSF, true); + status = pmSourceFitModel (source, modelPSF, PM_SOURCE_FIT_PSF); // skip poor fits Index: /branches/rel10_ifa/psModules/src/objects/pmPSFtry.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmPSFtry.h (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmPSFtry.h (revision 6848) @@ -6,6 +6,6 @@ * @author EAM, IfA * - * @version $Revision: 1.3.4.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-02-17 17:13:42 $ + * @version $Revision: 1.3.4.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -80,5 +80,6 @@ pmPSFtry *pmPSFtryAlloc( psArray *stars, ///< Add comment. - char *modelName ///< Add comment. + char *modelName, ///< Add comment. + bool poissonErrors // use poissonian or constant errors? ); @@ -94,5 +95,6 @@ psArray *sources, ///< Add comment. char *modelName, ///< Add comment. - float radius ///< Add comment. + float radius, ///< Add comment. + bool poissonErrors // use poissonian or constant errors? ); Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.c (revision 6848) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-10 20:21:36 $ + * @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -28,13 +28,16 @@ #include "pmSourceFitModel.h" -// save a static values so they may be set externally +// 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; - -bool pmSourceFitModelInit (float nIter, float tol) +static bool PM_PSF_POISSON_WEIGHTS = true; + +bool pmSourceFitModelInit (float nIter, float tol, bool poissonErrors) { PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; PM_SOURCE_FIT_MODEL_TOLERANCE = tol; + PM_PSF_POISSON_WEIGHTS = poissonErrors; + return true; } @@ -42,5 +45,5 @@ bool pmSourceFitModel (pmSource *source, pmModel *model, - const bool PSF) + pmSourceFitMode mode) { psTrace(__func__, 3, "---- %s() begin ----\n", __func__); @@ -65,5 +68,4 @@ pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); - int nParams = PSF ? 4 : params->n; psF32 dSky = source->moments->dSky; @@ -110,4 +112,50 @@ y->n = nPix; yErr->n = nPix; + + // XXX EAM : the new minimization API supplies the constraints as a struct + // XXX the chisq if a fcn of source flux. the minimization criterion should + // probably be scaled somehow by flux. measure the trend? it depends on + // the about of systematic error in the models themselves... + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + psMinConstrain *constrain = psMinConstrainAlloc(); + + // set parameter mask based on fitting mode + paramMask = psVectorAlloc (params->n, PS_TYPE_U8); + psVectorInit (paramMask, 1); + + int nParams = 0; + switch (mode) { + case PM_SOURCE_FIT_NORM: + // NORM-only model fits only source normalization (Io) + nParams = 1; + paramMask->data.U8[1] = 0; + break; + case PM_SOURCE_FIT_PSF: + // PSF model only fits x,y,Io + nParams = 3; + paramMask->data.U8[1] = 0; + paramMask->data.U8[2] = 0; + paramMask->data.U8[3] = 0; + break; + case PM_SOURCE_FIT_EXT: + // EXT model fits all params (except sky) + nParams = params->n - 1; + psVectorInit (paramMask, 0); + paramMask->data.U8[0] = 1; + break; + case PM_SOURCE_FIT_PSF_AND_SKY: + nParams = 4; + psAbort ("pmSourceFitModel", "PSF + SKY not currently available"); + break; + case PM_SOURCE_FIT_EXT_AND_SKY: + nParams = params->n; + psAbort ("pmSourceFitModel", "EXT + SKY not currently available"); + break; + default: + psAbort ("pmSourceFitModel", "invalid fitting mode"); + } + constrain->paramMask = paramMask; + if (nPix < nParams + 1) { psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); @@ -117,31 +165,9 @@ psFree (y); psFree (yErr); + psFree (myMin); + psFree (constrain); + psFree (paramMask); return(false); } - - // XXX EAM : the new minimization API supplies the constraints as a struct - // XXX the chisq if a fcn of source flux. the minimization criterion should - // probably be scaled somehow by flux. measure the trend? it depends on - // the about of systematic error in the models themselves... - psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, - PM_SOURCE_FIT_MODEL_TOLERANCE); - psMinConstrain *constrain = psMinConstrainAlloc(); - - // PSF model only fits x,y,Io, EXT model fits all (both skip sky) - paramMask = psVectorAlloc (params->n, PS_TYPE_U8); - paramMask->data.U8[0] = 1; - if (PSF) { - for (int i = 1; i < 4; i++) { - paramMask->data.U8[i] = 0; - } - for (int i = 4; i < paramMask->n; i++) { - paramMask->data.U8[i] = 1; - } - } else { - for (int i = 1; i < paramMask->n; i++) { - paramMask->data.U8[i] = 0; - } - } - constrain->paramMask = paramMask; // Set the parameter range checks @@ -210,2 +236,324 @@ return(rc); } + +/********************* Source Model Set Functions ***************************/ + +// these static variables are used by one pass of pmSourceFitSet to store +// data for a model set. If we are going to make psphot thread-safe, these +// will have to go in a structure of their own. +// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2, +// nPar = nSrc*(nOnePar - 1) + 1 +static pmModelFunc mFunc; +static int nPar; +static psVector *onePar; +static psVector *oneDeriv; + +bool pmModelFitSetInit (pmModelType type) +{ + + mFunc = pmModelFunc_GetFunction (type); + nPar = pmModelParameterCount (type); + + onePar = psVectorAlloc (nPar, PS_DATA_F32); + oneDeriv = psVectorAlloc (nPar, PS_DATA_F32); + + return true; +} + +void pmModelFitSetClear (void) +{ + + psFree (onePar); + psFree (oneDeriv); + return; +} + +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x) +{ + + psF32 value; + psF32 model; + + psF32 *PAR = onePar->data.F32; + psF32 *dPAR = oneDeriv->data.F32; + + psF32 *pars = params->data.F32; + psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32; + + int nSrc = (params->n - 1) / (nPar - 1); + + PAR[0] = model = pars[0]; + for (int i = 0; i < nSrc; i++) { + int nOff = i*nPar - i; + for (int n = 1; n < nPar; n++) { + PAR[n] = pars[nOff + n]; + } + if (deriv == NULL) { + value = mFunc (NULL, onePar, x); + } else { + value = mFunc (oneDeriv, onePar, x); + for (int n = 1; n < nPar; n++) { + dpars[nOff + n] = dPAR[n]; + } + } + model += value; + } + if (deriv != NULL) { + dpars[0] = dPAR[0]*2.0; + } + return (model); +} + +/* +i: 0 1 2 +n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6 +i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 +*/ + +bool pmSourceFitSet (pmSource *source, + psArray *modelSet, + pmSourceFitMode mode) +{ + psTrace(__func__, 3, "---- %s() begin ----\n", __func__); + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->moments, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + PS_ASSERT_PTR_NON_NULL(source->mask, false); + PS_ASSERT_PTR_NON_NULL(source->weight, false); + + psBool fitStatus = true; + psBool onPic = true; + psBool rc = true; + + // base values on first model + pmModel *model = modelSet->data[0]; + + // set the static variables + pmModelFitSetInit (model->type); + + int nSrc = modelSet->n; + int nPar = model->params->n - 1; // number of object parameters (excluding sky) + + // define parameter vectors for source set + psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + + pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); + + // define limits for a single-source model + psVector *oneDelta; + psVector *oneParMin; + psVector *oneParMax; + modelLimits (&oneDelta, &oneParMin, &oneParMax); + + psMinConstrain *constrain = psMinConstrainAlloc(); + constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + constrain->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8); + constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32); + + // set the parameter guesses and constraints for the multiple models + // first the values for the single sky parameter + params->data.F32[0] = model->params->data.F32[0]; + constrain->paramMin->data.F32[0] = oneParMin->data.F32[0]; + constrain->paramMax->data.F32[0] = oneParMax->data.F32[0]; + constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0]; + + // next, the values for the source parameters + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + for (int n = 1; n < nPar + 1; n++) { + params->data.F32[i*nPar + n] = model->params->data.F32[n]; + constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n]; + constrain->paramMin->data.F32[i*nPar + n] = oneParMin->data.F32[n]; + constrain->paramMax->data.F32[i*nPar + n] = oneParMax->data.F32[n]; + } + } + psFree (oneDelta); + psFree (oneParMin); + psFree (oneParMax); + + if (psTraceGetLevel(__func__) >= 5) { + for (int i = 0; i < params->n; i++) { + fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]); + } + } + + // set the parameter masks based on the fitting mode + int nParams = 0; + switch (mode) { + case PM_SOURCE_FIT_NORM: + // NORM-only model fits only source normalization (Io) + nParams = nSrc; + psVectorInit (constrain->paramMask, 1); + for (int i = 0; i < nSrc; i++) { + constrain->paramMask->data.U8[1 + i*nPar] = 0; + } + break; + case PM_SOURCE_FIT_PSF: + // PSF model only fits x,y,Io + nParams = 3*nSrc; + psVectorInit (constrain->paramMask, 1); + for (int i = 0; i < nSrc; i++) { + constrain->paramMask->data.U8[1 + i*nPar] = 0; + constrain->paramMask->data.U8[2 + i*nPar] = 0; + constrain->paramMask->data.U8[3 + i*nPar] = 0; + } + break; + case PM_SOURCE_FIT_EXT: + // EXT model fits all params (except sky) + nParams = nPar*nSrc; + psVectorInit (constrain->paramMask, 0); + constrain->paramMask->data.U8[0] = 1; + break; + case PM_SOURCE_FIT_PSF_AND_SKY: + nParams = 1 + 3*nSrc; + psAbort ("pmSourceFitModel", "PSF + SKY not currently available"); + break; + case PM_SOURCE_FIT_EXT_AND_SKY: + nParams = 1 + nPar*nSrc; + psAbort ("pmSourceFitModel", "EXT + SKY not currently available"); + break; + default: + psAbort ("pmSourceFitModel", "invalid fitting mode"); + } + + // maximum number of valid pixels + psS32 nPix = source->pixels->numRows * source->pixels->numCols; + + // local sky variance + psF32 dSky = source->moments->dSky; + + // construct the coordinate and value entries + psArray *x = psArrayAlloc(nPix); + psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); + psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); + + nPix = 0; + for (psS32 i = 0; i < source->pixels->numRows; i++) { + for (psS32 j = 0; j < source->pixels->numCols; j++) { + // skip masked points + if (source->mask->data.U8[i][j]) { + continue; + } + // skip zero-weight points + if (source->weight->data.F32[i][j] == 0) { + continue; + } + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + // Convert i/j to image space: + coord->data.F32[0] = (psF32) (j + source->pixels->col0); + coord->data.F32[1] = (psF32) (i + source->pixels->row0); + x->data[nPix] = (psPtr *) coord; + y->data.F32[nPix] = source->pixels->data.F32[i][j]; + // psMinimizeLMChi2 takes wt = 1/dY^2 + if (PM_PSF_POISSON_WEIGHTS) { + yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; + } else { + yErr->data.F32[nPix] = 1.0 / dSky; + } + nPix++; + } + } + if (nPix < nParams + 1) { + psTrace (__func__, 4, "insufficient valid pixels\n"); + psTrace(__func__, 3, "---- %s() end : fail pixels ----\n", __func__); + model->status = PM_MODEL_BADARGS; + psFree (x); + psFree (y); + psFree (yErr); + psFree (params); + psFree (dparams); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + return(false); + } + x->n = nPix; + y->n = nPix; + yErr->n = nPix; + + psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, + PM_SOURCE_FIT_MODEL_TOLERANCE); + + psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64); + + psTrace (__func__, 5, "fitting function\n"); + fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet); + + // parameter errors from the covariance matrix + for (int i = 0; i < dparams->n; i++) { + if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); + } + + // get the Gauss-Newton distance for fixed model parameters + if (constrain->paramMask != NULL) { + psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); + psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, pmModelFitSet); + for (int i = 0; i < dparams->n; i++) { + if (!constrain->paramMask->data.U8[i]) + continue; + dparams->data.F32[i] = delta->data.F64[i]; + } + psFree (delta); + } + + // assign back the parameters to the models + for (int i = 0; i < nSrc; i++) { + model = modelSet->data[i]; + model->params->data.F32[0] = params->data.F32[0]; + for (int n = 1; n < nPar + 1; n++) { + model->params->data.F32[n] = params->data.F32[i*nPar + n]; + model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n]; + } + // save the resulting chisq, nDOF, nIter + // these are not unique for any one source + model->chisq = myMin->value; + model->nIter = myMin->iter; + model->nDOF = y->n - nParams; + + // set the model success or failure status + model->status = fitStatus ? PM_MODEL_SUCCESS : PM_MODEL_NONCONVERGE; + + // models can go insane: reject these + onPic &= (model->params->data.F32[2] >= source->pixels->col0); + onPic &= (model->params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); + onPic &= (model->params->data.F32[3] >= source->pixels->row0); + onPic &= (model->params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); + if (!onPic) { + model->status = PM_MODEL_OFFIMAGE; + } + } + + source->mode |= PM_SOURCE_MODE_FITTED; + + psFree(x); + psFree(y); + psFree(yErr); + psFree(myMin); + psFree(covar); + psFree(constrain->paramMask); + psFree(constrain->paramMin); + psFree(constrain->paramMax); + psFree(constrain->paramDelta); + psFree(constrain); + psFree(params); + psFree(dparams); + + // free static memory used by pmModelFitSet + pmModelFitSetClear (); + + rc = (onPic && fitStatus); + psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF); + psTrace(__func__, 3, "---- %s end : status %d ----\n", __func__, rc); + return(rc); +} + Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitModel.h (revision 6848) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-09 03:14:23 $ + * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -13,16 +13,24 @@ # define PM_SOURCE_FIT_MODEL_H +typedef enum { + PM_SOURCE_FIT_NORM, + PM_SOURCE_FIT_PSF, + PM_SOURCE_FIT_EXT, + PM_SOURCE_FIT_PSF_AND_SKY, + PM_SOURCE_FIT_EXT_AND_SKY +} pmSourceFitMode; + bool pmSourceFitModelInit( float nIter, ///< max number of allowed iterations - float tol ///< convergence criterion + float tol, ///< convergence criterion + bool poissonErrors // use poisson errors for fits? ); /** pmSourceFitModel() * - * Fit the requested model to the specified source. The starting guess for the - * model is given by the input source.model parameter values. The pixels of - * interest are specified by the source.pixelsand source.maskentries. This - * function calls psMinimizeLMChi2() on the image data. The function returns TRUE - * on success or FALSE on failure. + * Fit the requested model to the specified source. The starting guess for the model is given + * by the input source.model parameter values. The pixels of interest are specified by the + * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image + * data. The function returns TRUE on success or FALSE on failure. * */ @@ -30,5 +38,31 @@ pmSource *source, ///< The input pmSource pmModel *model, ///< model to be fitted - const bool PSF ///< Treat model as PSF or EXT? + pmSourceFitMode mode ///< define parameters to be fitted +); + + +// initialize data for a group of object models +bool pmModelFitSetInit (pmModelType type); + +// clear data for a group of object models +void pmModelFitSetClear (void); + +// function used to fit a group of object models +psF32 pmModelFitSet(psVector *deriv, + const psVector *params, + const psVector *x); + +/** pmSourceFitSet() + * + * Fit the requested model to the specified source. The starting guess for the model is given + * by the input source.model parameter values. The pixels of interest are specified by the + * source.pixels and source.mask entries. This function calls psMinimizeLMChi2() on the image + * data. The function returns TRUE on success or FALSE on failure. + * + */ +bool pmSourceFitSet( + pmSource *source, ///< The input pmSource + psArray *modelSet, ///< model to be fitted + pmSourceFitMode mode ///< define parameters to be fitted ); Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.c (revision 6848) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-10 20:21:36 $ + * @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -33,4 +33,20 @@ static psVector *oneDeriv; +// 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 bool PM_PSF_POISSON_WEIGHTS = true; + +bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors) +{ + + PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; + PM_SOURCE_FIT_MODEL_TOLERANCE = tol; + + PM_PSF_POISSON_WEIGHTS = poissonErrors; + + return true; +} + bool pmModelFitSetInit (pmModelType type) { @@ -97,6 +113,8 @@ */ +/* # define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 # define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 +*/ bool pmSourceFitSet (pmSource *source, Index: /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceFitSet.h (revision 6848) @@ -3,6 +3,6 @@ * @author EAM, IfA; GLG, MHPCC * - * @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-03-07 06:33:35 $ + * @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -20,4 +20,6 @@ void pmModelFitSetClear (void); + +bool pmSourceFitSetInit (float nIter, float tol, bool poissonErrors); /** pmSourceFitSet() Index: /branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c =================================================================== --- /branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c (revision 6847) +++ /branches/rel10_ifa/psModules/src/objects/pmSourceIO_RAW.c (revision 6848) @@ -3,6 +3,6 @@ * @author EAM, IfA * - * @version $Revision: 1.1.2.3 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-04-01 02:47:20 $ + * @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-04-13 06:18:46 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -91,7 +91,8 @@ fprintf (f, "%9.6f ", dPAR[j]); } - fprintf (f, ": %8.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + fprintf (f, ": %8.4f %2d %#5x %7.3f %7.3f %7.1f %7.2f %4d %2d\n", source[0].apMag, source[0].type, source[0].mode, log10(model[0].chisq/model[0].nDOF), + log10(model[0].chisqNorm/model[0].nDOF), source[0].moments->SN, model[0].radiusTMP, @@ -146,8 +147,9 @@ fprintf (f, "%9.6f ", dPAR[j]); } - fprintf (f, ": %7.4f %2d %#5x %7.3f %7.1f %7.2f %4d %2d\n", + fprintf (f, ": %7.4f %2d %#5x %7.3f %7.3f %7.1f %7.2f %4d %2d\n", source->apMag, source[0].type, source[0].mode, log10(model[0].chisq/model[0].nDOF), + log10(model[0].chisqNorm/model[0].nDOF), source[0].moments->SN, model[0].radiusTMP, Index: /branches/rel10_ifa/psModules/src/pslib/psAdditionals.h =================================================================== --- /branches/rel10_ifa/psModules/src/pslib/psAdditionals.h (revision 6847) +++ /branches/rel10_ifa/psModules/src/pslib/psAdditionals.h (revision 6848) @@ -31,5 +31,5 @@ // strip whitespace from head and tail of string -int psStringStrip (char *string); +int psStringStrip (char *string); // write out header with NAXIS=0 Index: /branches/rel10_ifa/psModules/src/psmodules.h =================================================================== --- /branches/rel10_ifa/psModules/src/psmodules.h (revision 6847) +++ /branches/rel10_ifa/psModules/src/psmodules.h (revision 6848) @@ -63,5 +63,4 @@ # include # include -# include # include # include