IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 23, 2007, 2:11:02 PM (19 years ago)
Author:
magnier
Message:
  • added function pointers to pmModel to provide class-dependent functions
  • dropped pmModel*_GetFunction functions (use pmModel->func functions instead)
  • added modelParamsFromPSF functions to model classes
  • changed pmModelGroup to pmModelClass
  • added pmSourceFitSet.[ch]
  • updated pmSourceFitSet to allow variable model classes
  • added functions to add/sub and eval models & sources with an offset between image and chip
  • added function to set a pmModel flux
  • added function to instatiate a pmModel from a pmPSF at a coordinate
  • moved pmPSF I/O to chip->analysis from readout->analysis
  • changed pmPSF I/O function names from *ForPSFmodel to pmPSFmodel*
  • changed pmModel.params_NEW back to pmModel.params
  • changed pmFind*Peaks to pmPeaksIn* (* = Image,Vector)
  • dropped pmCullPeaks (deprecated)
  • changed pmModelSetType to pmModelClassGetType
  • changed pmModelGetType to pmModelClassGetName
  • fixed PGAUSS implementation of modelRadius function
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psModules/src/objects/pmSourceFitModel.c

    r13932 r14652  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2007-06-21 22:58:11 $
     8 *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2007-08-24 00:11:02 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    2727#include "pmGrowthCurve.h"
    2828#include "pmResiduals.h"
     29#include "pmPSF.h"
    2930#include "pmModel.h"
    30 #include "pmPSF.h"
    3131#include "pmSource.h"
    32 #include "pmModelGroup.h"
     32#include "pmModelClass.h"
    3333#include "pmSourceFitModel.h"
    3434
    … …  
    116116    psVector *dparams = model->dparams;
    117117
    118     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    119     if (!modelFunc)
    120         psAbort("invalid model function");
    121     pmModelLimits checkLimits = pmModelLimits_GetFunction (model->type);
    122     if (!checkLimits)
    123         psAbort("invalid model limits function");
    124 
    125118    // create the minimization constraints
    126119    psMinConstraint *constraint = psMinConstraintAlloc();
    127120    constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    128     constraint->checkLimits = checkLimits;
     121    constraint->checkLimits = model->modelLimits;
    129122
    130123    // set parameter mask based on fitting mode
    … …  
    156149    // force the floating parameters to fall within the contraint ranges
    157150    for (int i = 0; i < params->n; i++) {
    158         checkLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
    159         checkLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     151        model->modelLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
     152        model->modelLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
    160153    }
    161154
    … …  
    174167    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    175168
    176     fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, modelFunc);
     169    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, model->modelFunc);
    177170    for (int i = 0; i < dparams->n; i++) {
    178171        if (psTraceGetLevel("psModules.objects") >= 4) {
    … …  
    201194            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
    202195        }
    203         psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, modelFunc);
     196        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, model->modelFunc);
    204197        for (int i = 0; i < dparams->n; i++) {
    205198            if (!constraint->paramMask->data.U8[i])
    … …  
    237230}
    238231
    239 /********************* Source Model Set Functions ***************************/
    240 
    241 // these static variables are used by one pass of pmSourceFitSet to store
    242 // data for a model set.  If we are going to make psphot thread-safe, these
    243 // will have to go in a structure of their own or be allocated once per thread
    244 // sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2,
    245 // nPar = nSrc*(nOnePar - 1) + 1
    246 static pmModelFunc oneModelFunc;
    247 static pmModelLimits oneCheckLimits;
    248 static psVector *onePar;
    249 static psVector *oneDeriv;
    250 static int nOnePar;
    251 
    252 bool pmSourceFitSetInit (pmModelType type)
    253 {
    254 
    255     oneModelFunc = pmModelFunc_GetFunction (type);
    256     oneCheckLimits = pmModelLimits_GetFunction (type);
    257     nOnePar = pmModelParameterCount (type);
    258 
    259     onePar = psVectorAlloc (nOnePar, PS_DATA_F32);
    260     oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32);
    261 
    262     return true;
    263 }
    264 
    265 void pmSourceFitSetClear (void)
    266 {
    267 
    268     psFree (onePar);
    269     psFree (oneDeriv);
    270     return;
    271 }
    272 
    273 bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)
    274 {
    275     // convert the value of nParam into corresponding single model parameter entry
    276     // convert params into corresponding single model parameter pointer
    277 
    278     int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1;
    279     float *paramSingle = params + nParam - nParamSingle;
    280     float *betaSingle = betas + nParam - nParamSingle;
    281     bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle);
    282     return status;
    283 }
    284 
    285 psF32 pmSourceFitSet_Function(psVector *deriv,
    286                               const psVector *params,
    287                               const psVector *x)
    288 {
    289 
    290     psF32 value;
    291     psF32 model;
    292 
    293     psF32 *PAR = onePar->data.F32;
    294     psF32 *dPAR = oneDeriv->data.F32;
    295 
    296     psF32 *pars = params->data.F32;
    297     psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;
    298 
    299     int nSrc = (params->n - 1) / (nOnePar - 1);
    300 
    301     PAR[0] = model = pars[0];
    302     for (int i = 0; i < nSrc; i++) {
    303         int nOff = i*nOnePar - i;
    304         for (int n = 1; n < nOnePar; n++) {
    305             PAR[n] = pars[nOff + n];
    306         }
    307         if (deriv == NULL) {
    308             value = oneModelFunc (NULL, onePar, x);
    309         } else {
    310             value = oneModelFunc (oneDeriv, onePar, x);
    311             for (int n = 1; n < nOnePar; n++) {
    312                 dpars[nOff + n] = dPAR[n];
    313             }
    314         }
    315         model += value;
    316     }
    317     if (deriv != NULL) {
    318         dpars[0] = dPAR[0]*2.0;
    319     }
    320     return (model);
    321 }
    322 
    323 /*
    324 i:       0                   1                 2
    325 n:         1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6
    326 i*6 + n: 0 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
    327 */
    328 
    329 bool pmSourceFitSet (pmSource *source,
    330                      psArray *modelSet,
    331                      pmSourceFitMode mode,
    332                      psMaskType maskVal)
    333 {
    334     psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
    335     PS_ASSERT_PTR_NON_NULL(source, false);
    336     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    337     PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
    338     PS_ASSERT_PTR_NON_NULL(source->weight, false);
    339 
    340     psBool fitStatus = true;
    341     psBool onPic     = true;
    342     psBool rc        = true;
    343 
    344     // maximum number of valid pixels
    345     psS32 nPix = source->pixels->numRows * source->pixels->numCols;
    346 
    347     // construct the coordinate and value entries
    348     psArray *x = psArrayAllocEmpty(nPix);
    349     psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    350     psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    351 
    352     // fill in the coordinate and value entries
    353     nPix = 0;
    354     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    355         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    356             // skip masked points
    357             if (source->maskObj->data.U8[i][j] & maskVal) {
    358                 continue;
    359             }
    360             // skip zero-weight points
    361             if (source->weight->data.F32[i][j] == 0) {
    362                 continue;
    363             }
    364             // skip nan values in image
    365             if (!isfinite(source->pixels->data.F32[i][j])) {
    366                 continue;
    367             }
    368 
    369             psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    370 
    371             // Convert i/j to image space:
    372             coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    373             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    374             x->data[nPix] = (psPtr *) coord;
    375             y->data.F32[nPix] = source->pixels->data.F32[i][j];
    376 
    377             // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
    378             // as weight to avoid the bias from systematic errors here we would just use the
    379             // source sky variance
    380             if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
    381                 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
    382             } else {
    383                 yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
    384             }
    385             nPix++;
    386         }
    387     }
    388     x->n = nPix;
    389     y->n = nPix;
    390     yErr->n = nPix;
    391 
    392     // base values on first model (all models must be identical)
    393     pmModel *model = modelSet->data[0];
    394 
    395     // determine number of model parameters
    396     int nSrc = modelSet->n;
    397     int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
    398 
    399     // define parameter vectors for source set
    400     psVector *params = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    401     psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    402 
    403     // set the static variables
    404     pmSourceFitSetInit (model->type);
    405 
    406     // create the minimization constraints
    407     psMinConstraint *constraint = psMinConstraintAlloc();
    408     constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
    409     constraint->checkLimits = pmSourceFitSet_CheckLimits;
    410 
    411     // set the parameter guesses for the multiple models
    412     // first the value for the single sky parameter
    413     params->data.F32[0] = model->params->data.F32[0];
    414 
    415     // next, the values for the source parameters
    416     for (int i = 0; i < nSrc; i++) {
    417         model = modelSet->data[i];
    418         for (int n = 1; n < nPar + 1; n++) {
    419             params->data.F32[i*nPar + n] = model->params->data.F32[n];
    420         }
    421     }
    422 
    423     if (psTraceGetLevel("psModules.objects") >= 5) {
    424         for (int i = 0; i < params->n; i++) {
    425             fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
    426         }
    427     }
    428 
    429     // set the parameter masks based on the fitting mode
    430     int nParams = 0;
    431     switch (mode) {
    432     case PM_SOURCE_FIT_NORM:
    433         // NORM-only model fits only source normalization (Io)
    434         nParams = nSrc;
    435         psVectorInit (constraint->paramMask, 1);
    436         for (int i = 0; i < nSrc; i++) {
    437             constraint->paramMask->data.U8[1 + i*nPar] = 0;
    438         }
    439         break;
    440     case PM_SOURCE_FIT_PSF:
    441         // PSF model only fits x,y,Io
    442         nParams = 3*nSrc;
    443         psVectorInit (constraint->paramMask, 1);
    444         for (int i = 0; i < nSrc; i++) {
    445             constraint->paramMask->data.U8[1 + i*nPar] = 0;
    446             constraint->paramMask->data.U8[2 + i*nPar] = 0;
    447             constraint->paramMask->data.U8[3 + i*nPar] = 0;
    448         }
    449         break;
    450     case PM_SOURCE_FIT_EXT:
    451         // EXT model fits all params (except sky)
    452         nParams = nPar*nSrc;
    453         psVectorInit (constraint->paramMask, 0);
    454         constraint->paramMask->data.U8[0] = 1;
    455         break;
    456     default:
    457         psAbort("invalid fitting mode");
    458     }
    459 
    460     // force the floating parameters to fall within the contraint ranges
    461     for (int i = 0; i < params->n; i++) {
    462         pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
    463         pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
    464     }
    465 
    466     if (nPix <  nParams + 1) {
    467         psTrace (__func__, 4, "insufficient valid pixels\n");
    468         psTrace("psModules.objects", 3, "---- %s() end : fail pixels ----\n", __func__);
    469         model->flags |= PM_MODEL_STATUS_BADARGS;
    470         psFree (x);
    471         psFree (y);
    472         psFree (yErr);
    473         psFree (params);
    474         psFree (dparams);
    475         psFree(constraint);
    476         pmSourceFitSetClear ();
    477         return(false);
    478     }
    479 
    480     psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
    481 
    482     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    483 
    484     fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function);
    485     if (!fitStatus) {
    486         psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);
    487     }
    488 
    489     // parameter errors from the covariance matrix
    490     for (int i = 0; i < dparams->n; i++) {
    491         if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
    492             continue;
    493         dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    494     }
    495 
    496     // get the Gauss-Newton distance for fixed model parameters
    497     if (constraint->paramMask != NULL) {
    498         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
    499         psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
    500         altmask->data.U8[0] = 1;
    501         for (int i = 1; i < dparams->n; i++) {
    502             altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
    503         }
    504         psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function);
    505         for (int i = 0; i < dparams->n; i++) {
    506             if (!constraint->paramMask->data.U8[i])
    507                 continue;
    508             // note that delta is the value *subtracted* from the parameter
    509             // to get the new guess.  for dparams to represent the direction
    510             // of motion, we need to take -delta
    511             dparams->data.F32[i] = -delta->data.F32[i];
    512         }
    513         psFree (delta);
    514         psFree (altmask);
    515     }
    516 
    517     // assign back the parameters to the models
    518     for (int i = 0; i < nSrc; i++) {
    519         model = modelSet->data[i];
    520         model->params->data.F32[0] = params->data.F32[0];
    521         for (int n = 1; n < nPar + 1; n++) {
    522             if (psTraceGetLevel("psModules.objects") >= 4) {
    523                 fprintf (stderr, "%f ", params->data.F32[i*nPar + n]);
    524             }
    525             model->params->data.F32[n] = params->data.F32[i*nPar + n];
    526             model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
    527         }
    528         psTrace ("psModules.objects", 4, " src %d", i);
    529 
    530         // save the resulting chisq, nDOF, nIter
    531         // these are not unique for any one source
    532         model->chisq = myMin->value;
    533         model->nIter = myMin->iter;
    534         model->nDOF  = y->n - nParams;
    535 
    536         // set the model success or failure status
    537         model->flags |= PM_MODEL_STATUS_FITTED;
    538         if (!fitStatus) model->flags |= PM_MODEL_STATUS_NONCONVERGE;
    539 
    540         // models can go insane: reject these
    541         onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
    542         onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
    543         onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);
    544         onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->row0 + source->pixels->numRows);
    545         if (!onPic) model->flags |= PM_MODEL_STATUS_OFFIMAGE;
    546     }
    547     psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    548 
    549     source->mode |= PM_SOURCE_MODE_FITTED;
    550 
    551     psFree(x);
    552     psFree(y);
    553     psFree(yErr);
    554     psFree(myMin);
    555     psFree(covar);
    556     psFree(constraint);
    557     psFree(params);
    558     psFree(dparams);
    559 
    560     // free static memory used by pmSourceFitSet
    561     pmSourceFitSetClear ();
    562 
    563     rc = (onPic && fitStatus);
    564     psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
    565     psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);
    566     return(rc);
    567 }
    568 
Note: See TracChangeset for help on using the changeset viewer.