- Timestamp:
- Aug 19, 2007, 3:58:16 PM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branch_20070817/psModules/src/objects/pmSourceFitModel.c
r14544 r14546 6 6 * @author GLG, MHPCC 7 7 * 8 * @version $Revision: 1.24.4. 1$ $Name: not supported by cvs2svn $9 * @date $Date: 2007-08- 17 21:01:59$8 * @version $Revision: 1.24.4.2 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2007-08-20 01:58:16 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 230 230 } 231 231 232 # define SKIP_FIT_SET 1233 # if (SKIP_FIT_SET)234 235 bool pmSourceFitSet (pmSource *source,236 psArray *modelSet,237 pmSourceFitMode mode,238 psMaskType maskVal)239 {240 psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);241 PS_ASSERT_PTR_NON_NULL(source, false);242 PS_ASSERT_PTR_NON_NULL(source->pixels, false);243 PS_ASSERT_PTR_NON_NULL(source->maskObj, false);244 PS_ASSERT_PTR_NON_NULL(source->weight, false);245 246 return true;247 }248 249 # else250 251 /********************* Source Model Set Functions ***************************/252 253 // these static variables are used by one pass of pmSourceFitSet to store254 // data for a model set. If we are going to make psphot thread-safe, these255 // will have to go in a structure of their own or be allocated once per thread256 // sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2,257 // nPar = nSrc*(nOnePar - 1) + 1258 static pmModelFunc oneModelFunc;259 static pmModelLimits oneCheckLimits;260 static psVector *onePar;261 static psVector *oneDeriv;262 static int nOnePar;263 264 bool pmSourceFitSetInit (pmModelType type)265 {266 267 oneModelFunc = pmModelFunc_GetFunction (type);268 oneCheckLimits = pmModelLimits_GetFunction (type);269 nOnePar = pmModelClassParameterCount (type);270 271 onePar = psVectorAlloc (nOnePar, PS_DATA_F32);272 oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32);273 274 return true;275 }276 277 void pmSourceFitSetClear (void)278 {279 280 psFree (onePar);281 psFree (oneDeriv);282 return;283 }284 285 bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)286 {287 // convert the value of nParam into corresponding single model parameter entry288 // convert params into corresponding single model parameter pointer289 290 int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1;291 float *paramSingle = params + nParam - nParamSingle;292 float *betaSingle = betas + nParam - nParamSingle;293 bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle);294 return status;295 }296 297 psF32 pmSourceFitSet_Function(psVector *deriv,298 const psVector *params,299 const psVector *x)300 {301 302 psF32 value;303 psF32 model;304 305 psF32 *PAR = onePar->data.F32;306 psF32 *dPAR = oneDeriv->data.F32;307 308 psF32 *pars = params->data.F32;309 psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;310 311 int nSrc = (params->n - 1) / (nOnePar - 1);312 313 PAR[0] = model = pars[0];314 for (int i = 0; i < nSrc; i++) {315 int nOff = i*nOnePar - i;316 for (int n = 1; n < nOnePar; n++) {317 PAR[n] = pars[nOff + n];318 }319 if (deriv == NULL) {320 value = oneModelFunc (NULL, onePar, x);321 } else {322 value = oneModelFunc (oneDeriv, onePar, x);323 for (int n = 1; n < nOnePar; n++) {324 dpars[nOff + n] = dPAR[n];325 }326 }327 model += value;328 }329 if (deriv != NULL) {330 dpars[0] = dPAR[0]*2.0;331 }332 return (model);333 }334 335 /*336 i: 0 1 2337 n: 1 2 3 4 5 6 1 2 3 4 5 6 1 2 3 4 5 6338 i*6 + n: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18339 */340 341 bool pmSourceFitSet (pmSource *source,342 psArray *modelSet,343 pmSourceFitMode mode,344 psMaskType maskVal)345 {346 psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);347 PS_ASSERT_PTR_NON_NULL(source, false);348 PS_ASSERT_PTR_NON_NULL(source->pixels, false);349 PS_ASSERT_PTR_NON_NULL(source->maskObj, false);350 PS_ASSERT_PTR_NON_NULL(source->weight, false);351 352 psBool fitStatus = true;353 psBool onPic = true;354 psBool rc = true;355 356 // maximum number of valid pixels357 psS32 nPix = source->pixels->numRows * source->pixels->numCols;358 359 // construct the coordinate and value entries360 psArray *x = psArrayAllocEmpty(nPix);361 psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);362 psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);363 364 // fill in the coordinate and value entries365 nPix = 0;366 for (psS32 i = 0; i < source->pixels->numRows; i++) {367 for (psS32 j = 0; j < source->pixels->numCols; j++) {368 // skip masked points369 if (source->maskObj->data.U8[i][j] & maskVal) {370 continue;371 }372 // skip zero-weight points373 if (source->weight->data.F32[i][j] == 0) {374 continue;375 }376 // skip nan values in image377 if (!isfinite(source->pixels->data.F32[i][j])) {378 continue;379 }380 381 psVector *coord = psVectorAlloc(2, PS_TYPE_F32);382 383 // Convert i/j to image space:384 coord->data.F32[0] = (psF32) (j + source->pixels->col0);385 coord->data.F32[1] = (psF32) (i + source->pixels->row0);386 x->data[nPix] = (psPtr *) coord;387 y->data.F32[nPix] = source->pixels->data.F32[i][j];388 389 // psMinimizeLMChi2 takes wt = 1/dY^2. suggestion from RHL is to use the local sky390 // as weight to avoid the bias from systematic errors here we would just use the391 // source sky variance392 if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {393 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];394 } else {395 yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;396 }397 nPix++;398 }399 }400 x->n = nPix;401 y->n = nPix;402 yErr->n = nPix;403 404 // base values on first model (all models must be identical)405 pmModel *model = modelSet->data[0];406 407 // determine number of model parameters408 int nSrc = modelSet->n;409 int nPar = model->params->n - 1; // number of object parameters (excluding sky)410 411 // define parameter vectors for source set412 psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);413 psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);414 415 // set the static variables416 pmSourceFitSetInit (model->type);417 418 // create the minimization constraints419 psMinConstraint *constraint = psMinConstraintAlloc();420 constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);421 constraint->checkLimits = pmSourceFitSet_CheckLimits;422 423 // set the parameter guesses for the multiple models424 // first the value for the single sky parameter425 params->data.F32[0] = model->params->data.F32[0];426 427 // next, the values for the source parameters428 for (int i = 0; i < nSrc; i++) {429 model = modelSet->data[i];430 for (int n = 1; n < nPar + 1; n++) {431 params->data.F32[i*nPar + n] = model->params->data.F32[n];432 }433 }434 435 if (psTraceGetLevel("psModules.objects") >= 5) {436 for (int i = 0; i < params->n; i++) {437 fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);438 }439 }440 441 // set the parameter masks based on the fitting mode442 int nParams = 0;443 switch (mode) {444 case PM_SOURCE_FIT_NORM:445 // NORM-only model fits only source normalization (Io)446 nParams = nSrc;447 psVectorInit (constraint->paramMask, 1);448 for (int i = 0; i < nSrc; i++) {449 constraint->paramMask->data.U8[1 + i*nPar] = 0;450 }451 break;452 case PM_SOURCE_FIT_PSF:453 // PSF model only fits x,y,Io454 nParams = 3*nSrc;455 psVectorInit (constraint->paramMask, 1);456 for (int i = 0; i < nSrc; i++) {457 constraint->paramMask->data.U8[1 + i*nPar] = 0;458 constraint->paramMask->data.U8[2 + i*nPar] = 0;459 constraint->paramMask->data.U8[3 + i*nPar] = 0;460 }461 break;462 case PM_SOURCE_FIT_EXT:463 // EXT model fits all params (except sky)464 nParams = nPar*nSrc;465 psVectorInit (constraint->paramMask, 0);466 constraint->paramMask->data.U8[0] = 1;467 break;468 default:469 psAbort("invalid fitting mode");470 }471 472 // force the floating parameters to fall within the contraint ranges473 for (int i = 0; i < params->n; i++) {474 pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);475 pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);476 }477 478 if (nPix < nParams + 1) {479 psTrace (__func__, 4, "insufficient valid pixels\n");480 psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__);481 model->flags |= PM_MODEL_STATUS_BADARGS;482 psFree (x);483 psFree (y);484 psFree (yErr);485 psFree (params);486 psFree (dparams);487 psFree(constraint);488 pmSourceFitSetClear ();489 return(false);490 }491 492 psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);493 494 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);495 496 fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function);497 if (!fitStatus) {498 psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);499 }500 501 // parameter errors from the covariance matrix502 for (int i = 0; i < dparams->n; i++) {503 if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])504 continue;505 dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);506 }507 508 // get the Gauss-Newton distance for fixed model parameters509 if (constraint->paramMask != NULL) {510 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);511 psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);512 altmask->data.U8[0] = 1;513 for (int i = 1; i < dparams->n; i++) {514 altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;515 }516 psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function);517 for (int i = 0; i < dparams->n; i++) {518 if (!constraint->paramMask->data.U8[i])519 continue;520 // note that delta is the value *subtracted* from the parameter521 // to get the new guess. for dparams to represent the direction522 // of motion, we need to take -delta523 dparams->data.F32[i] = -delta->data.F32[i];524 }525 psFree (delta);526 psFree (altmask);527 }528 529 // assign back the parameters to the models530 for (int i = 0; i < nSrc; i++) {531 model = modelSet->data[i];532 model->params->data.F32[0] = params->data.F32[0];533 for (int n = 1; n < nPar + 1; n++) {534 if (psTraceGetLevel("psModules.objects") >= 4) {535 fprintf (stderr, "%f ", params->data.F32[i*nPar + n]);536 }537 model->params->data.F32[n] = params->data.F32[i*nPar + n];538 model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];539 }540 psTrace ("psModules.objects", 4, " src %d", i);541 542 // save the resulting chisq, nDOF, nIter543 // these are not unique for any one source544 model->chisq = myMin->value;545 model->nIter = myMin->iter;546 model->nDOF = y->n - nParams;547 548 // set the model success or failure status549 model->flags |= PM_MODEL_STATUS_FITTED;550 if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;551 552 // models can go insane: reject these553 onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);554 onPic &= (model->params->data.F32[PM_PAR_XPOS] < source->pixels->col0 + source->pixels->numCols);555 onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);556 onPic &= (model->params->data.F32[PM_PAR_XPOS] < source->pixels->row0 + source->pixels->numRows);557 if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE;558 }559 psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);560 561 source->mode |= PM_SOURCE_MODE_FITTED;562 563 psFree(x);564 psFree(y);565 psFree(yErr);566 psFree(myMin);567 psFree(covar);568 psFree(constraint);569 psFree(params);570 psFree(dparams);571 572 // free static memory used by pmSourceFitSet573 pmSourceFitSetClear ();574 575 rc = (onPic && fitStatus);576 psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);577 psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);578 return(rc);579 }580 581 # endif
Note:
See TracChangeset
for help on using the changeset viewer.
