IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 28, 2006, 4:41:31 PM (19 years ago)
Author:
magnier
Message:

removed dSky derived from moments: this should be supplied externally
PM_PSF_POISSON_WEIGHTS to PM_SOURCE_FIT_MODEL_PIX_WEIGHTS
changed pmSourceModelFit to use the new dynamic parameter limits for psMinimizeLMM

changed pmSourceFitSet to use the new dynamic parameter limits for psMinimizeLMM

  • this required a new internal function pmSourceFitSet_CheckLimits

changed old functions with names like pmModelFitSet to pmSourceFitSet

changed all psMinimizeLMM related data to F32

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r10199 r10259  
    66 *  @author GLG, MHPCC
    77 *
    8  *  @version $Revision: 1.17 $ $Name: not supported by cvs2svn $
    9  *  @date $Date: 2006-11-26 22:26:51 $
     8 *  @version $Revision: 1.18 $ $Name: not supported by cvs2svn $
     9 *  @date $Date: 2006-11-29 02:41:31 $
    1010 *
    1111 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    3636static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
    3737static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
    38 static bool PM_PSF_POISSON_WEIGHTS = true;
     38static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
    3939
    4040bool pmSourceFitModelInit (float nIter, float tol, float weight, bool poissonErrors)
    … …  
    4444    PM_SOURCE_FIT_MODEL_TOLERANCE = tol;
    4545    PM_SOURCE_FIT_MODEL_WEIGHT = weight;
    46     PM_PSF_POISSON_WEIGHTS = poissonErrors;
     46    PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = poissonErrors;
    4747
    4848    return true;
    … …  
    5353                       pmSourceFitMode mode)
    5454{
    55     psTrace("psModules.objects", 5, "---- %s() begin ----\n", __func__);
     55    psTrace("psModules.objects", 5, "---- %s begin ----\n", __func__);
    5656    PS_ASSERT_PTR_NON_NULL(source, false);
    57     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    5857    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    5958    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    6059    PS_ASSERT_PTR_NON_NULL(source->weight, false);
    6160
    62     // XXX EAM : is it necessary for the mask & weight to exist?  the
    63     //           tests below could be conditions (!NULL)
    64 
    6561    psBool fitStatus = true;
    6662    psBool onPic     = true;
    6763    psBool rc        = true;
    6864
    69     psVector *params = model->params;
    70     psVector *dparams = model->dparams;
    71     psVector *paramMask = NULL;
    72 
    73     pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
    74 
    75     psF32 dSky = source->moments->dSky;
    76 
    7765    // maximum number of valid pixels
    7866    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
    7967
    80     // construct the coordinate and value entries
     68    // arrays to hold the data to be fitted
    8169    psArray *x = psArrayAllocEmpty(nPix);
    8270    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    8371    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    8472
     73    // fill in the coordinate and value entries
    8574    nPix = 0;
    8675    for (psS32 i = 0; i < source->pixels->numRows; i++) {
    … …  
    10291            x->data[nPix] = (psPtr *) coord;
    10392            y->data.F32[nPix] = source->pixels->data.F32[i][j];
    104             // psMinimizeLMChi2 takes wt = 1/dY^2
    105             if (PM_PSF_POISSON_WEIGHTS) {
     93
     94            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
     95            // as weight to avoid the bias from systematic errors here we would just use the
     96            // source sky variance
     97            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
    10698                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
    10799            } else {
    108                 yErr->data.F32[nPix] = 1.0 / dSky;
    109             }
    110             // XXX EAM : suggestion from RHL is to use the local sky as weight
    111             // to avoid the bias from systematic errors
    112             // here we would just use the source sky variance
     100                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
     101            }
    113102            nPix++;
    114103        }
    … …  
    118107    yErr->n = nPix;
    119108
    120     // XXX EAM : the new minimization API supplies the constraints as a struct
    121     // XXX the chisq if a fcn of source flux. the minimization criterion should
    122     // probably  be scaled somehow by flux.  measure the trend?  it depends on
    123     // the about of systematic error in the models themselves...
    124     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    125                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    126     psMinConstrain *constrain = psMinConstrainAlloc();
     109    psVector *params = model->params;
     110    psVector *dparams = model->dparams;
     111
     112    pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type);
     113    if (!modelFunc)
     114        psAbort ("pmSourceFitModel", "invalid model function");
     115    pmModelLimits checkLimits = pmModelLimits_GetFunction (model->type);
     116    if (!checkLimits)
     117        psAbort ("pmSourceFitModel", "invalid model limits function");
     118
     119    // create the minimization constraints
     120    psMinConstraint *constraint = psMinConstraintAlloc();
     121    constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
     122    constraint->checkLimits = checkLimits;
    127123
    128124    // set parameter mask based on fitting mode
    129     paramMask = psVectorAlloc (params->n, PS_TYPE_U8);
    130     psVectorInit (paramMask, 1);
    131 
    132125    int nParams = 0;
    133126    switch (mode) {
    … …  
    135128        // NORM-only model fits only source normalization (Io)
    136129        nParams = 1;
    137         paramMask->data.U8[PM_PAR_I0] = 0;
     130        psVectorInit (constraint->paramMask, 1);
     131        constraint->paramMask->data.U8[PM_PAR_I0] = 0;
    138132        break;
    139133    case PM_SOURCE_FIT_PSF:
    140134        // PSF model only fits x,y,Io
    141135        nParams = 3;
    142         paramMask->data.U8[PM_PAR_I0] = 0;
    143         paramMask->data.U8[PM_PAR_XPOS] = 0;
    144         paramMask->data.U8[PM_PAR_YPOS] = 0;
     136        psVectorInit (constraint->paramMask, 1);
     137        constraint->paramMask->data.U8[PM_PAR_I0] = 0;
     138        constraint->paramMask->data.U8[PM_PAR_XPOS] = 0;
     139        constraint->paramMask->data.U8[PM_PAR_YPOS] = 0;
    145140        break;
    146141    case PM_SOURCE_FIT_EXT:
    147142        // EXT model fits all params (except sky)
    148143        nParams = params->n - 1;
    149         psVectorInit (paramMask, 0);
    150         paramMask->data.U8[PM_PAR_SKY] = 1;
    151         break;
    152     case PM_SOURCE_FIT_PSF_AND_SKY:
    153         nParams = 4;
    154         psAbort ("pmSourceFitModel", "PSF + SKY not currently available");
    155         break;
    156     case PM_SOURCE_FIT_EXT_AND_SKY:
    157         nParams = params->n;
    158         psAbort ("pmSourceFitModel", "EXT + SKY not currently available");
     144        psVectorInit (constraint->paramMask, 0);
     145        constraint->paramMask->data.U8[PM_PAR_SKY] = 1;
    159146        break;
    160147    default:
    161148        psAbort ("pmSourceFitModel", "invalid fitting mode");
    162149    }
    163     constrain->paramMask = paramMask;
     150    // force the floating parameters to fall within the contraint ranges
     151    for (int i = 0; i < params->n; i++) {
     152        checkLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
     153        checkLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     154    }
    164155
    165156    if (nPix <  nParams + 1) {
    166157        psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
    167         psTrace ("psModules.objects", 5, "---- %s(false) end ----\n", __func__);
    168158        model->status = PM_MODEL_BADARGS;
    169159        psFree (x);
    170160        psFree (y);
    171161        psFree (yErr);
    172         psFree (myMin);
    173         psFree (constrain);
    174         psFree (paramMask);
     162        psFree (constraint);
    175163        return(false);
    176164    }
    177165
    178     // Set the parameter range checks
    179     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    180     modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
    181 
    182     // force the floating parameters to fall within the contraint ranges
    183     for (int i = 0; i < params->n; i++) {
    184         if (constrain->paramMask->data.U8[i])
     166    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
     167
     168    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     169
     170    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, modelFunc);
     171    for (int i = 0; i < dparams->n; i++) {
     172        if (psTraceGetLevel("psModules.objects") >= 4) {
     173            fprintf (stderr, "%f ", params->data.F32[i]);
     174        }
     175        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
    185176            continue;
    186         params->data.F32[i] = PS_MIN(constrain->paramMax->data.F32[i], params->data.F32[i]);
    187         params->data.F32[i] = PS_MAX(constrain->paramMin->data.F32[i], params->data.F32[i]);
    188     }
    189 
    190     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    191 
    192     psTrace (__func__, 5, "fitting function\n");
    193 
    194     fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc);
    195     for (int i = 0; i < dparams->n; i++) {
    196         psTrace ("psModules.objects", 4, "%f ", params->data.F32[i]);
    197         if ((paramMask != NULL) && paramMask->data.U8[i])
    198             continue;
    199         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     177        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    200178    }
    201179    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    … …  
    208186    // get the Gauss-Newton distance for fixed model parameters
    209187    // hold the fitted parameters fixed; mask sky which is not fitted at all
    210     if (paramMask != NULL) {
    211         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
     188    if (constraint->paramMask != NULL) {
     189        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
    212190        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
    213191        altmask->data.U8[0] = 1;
    214192        for (int i = 1; i < dparams->n; i++) {
    215             altmask->data.U8[i] = (paramMask->data.U8[i]) ? 0 : 1;
     193            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
    216194        }
    217195        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, modelFunc);
    218196        for (int i = 0; i < dparams->n; i++) {
    219             if (!paramMask->data.U8[i])
     197            if (!constraint->paramMask->data.U8[i])
    220198                continue;
    221199            // note that delta is the value *subtracted* from the parameter
    222200            // to get the new guess.  for dparams to represent the direction
    223201            // of motion, we need to take -delta
    224             dparams->data.F32[i] = -delta->data.F64[i];
     202            dparams->data.F32[i] = -delta->data.F32[i];
    225203        }
    226204        psFree (delta);
    … …  
    251229    psFree(myMin);
    252230    psFree(covar);
    253     psFree(constrain->paramMask);
    254     psFree(constrain->paramMin);
    255     psFree(constrain->paramMax);
    256     psFree(constrain->paramDelta);
    257     psFree(constrain);
     231    psFree(constraint);
    258232
    259233    rc = (onPic && fitStatus);
    … …  
    266240// these static variables are used by one pass of pmSourceFitSet to store
    267241// data for a model set.  If we are going to make psphot thread-safe, these
    268 // will have to go in a structure of their own.
     242// will have to go in a structure of their own or be allocated once per thread
    269243// sky, p1.1, p1.2, p1.3,... p1.n, p2.1, p2.2,
    270244// nPar = nSrc*(nOnePar - 1) + 1
    271 static pmModelFunc mFunc;
    272 static int nPar;
     245static pmModelFunc oneModelFunc;
     246static pmModelLimits oneCheckLimits;
    273247static psVector *onePar;
    274248static psVector *oneDeriv;
    275 
    276 bool pmModelFitSetInit (pmModelType type)
    277 {
    278 
    279     mFunc = pmModelFunc_GetFunction (type);
    280     nPar  = pmModelParameterCount (type);
    281 
    282     onePar = psVectorAlloc (nPar, PS_DATA_F32);
    283     oneDeriv = psVectorAlloc (nPar, PS_DATA_F32);
     249static int nOnePar;
     250
     251bool pmSourceFitSetInit (pmModelType type)
     252{
     253
     254    oneModelFunc = pmModelFunc_GetFunction (type);
     255    oneCheckLimits = pmModelLimits_GetFunction (type);
     256    nOnePar = pmModelParameterCount (type);
     257
     258    onePar = psVectorAlloc (nOnePar, PS_DATA_F32);
     259    oneDeriv = psVectorAlloc (nOnePar, PS_DATA_F32);
    284260
    285261    return true;
    286262}
    287263
    288 void pmModelFitSetClear (void)
     264void pmSourceFitSetClear (void)
    289265{
    290266
    … …  
    294270}
    295271
    296 psF32 pmModelFitSet(psVector *deriv,
    297                     const psVector *params,
    298                     const psVector *x)
     272bool pmSourceFitSet_CheckLimits (psMinConstraintMode mode, int nParam, float *params, float *betas)
     273{
     274    // convert the value of nParam into corresponding single model parameter entry
     275    // convert params into corresponding single model parameter pointer
     276
     277    int nParamSingle = (nParam - 1) % (nOnePar - 1) + 1;
     278    float *paramSingle = params + nParam - nParamSingle;
     279    float *betaSingle = betas + nParam - nParamSingle;
     280    bool status = oneCheckLimits (mode, nParamSingle, paramSingle, betaSingle);
     281    return status;
     282}
     283
     284psF32 pmSourceFitSet_Function(psVector *deriv,
     285                              const psVector *params,
     286                              const psVector *x)
    299287{
    300288
    … …  
    308296    psF32 *dpars = (deriv == NULL) ? NULL : deriv->data.F32;
    309297
    310     int nSrc = (params->n - 1) / (nPar - 1);
     298    int nSrc = (params->n - 1) / (nOnePar - 1);
    311299
    312300    PAR[0] = model = pars[0];
    313301    for (int i = 0; i < nSrc; i++) {
    314         int nOff = i*nPar - i;
    315         for (int n = 1; n < nPar; n++) {
     302        int nOff = i*nOnePar - i;
     303        for (int n = 1; n < nOnePar; n++) {
    316304            PAR[n] = pars[nOff + n];
    317305        }
    318306        if (deriv == NULL) {
    319             value = mFunc (NULL, onePar, x);
     307            value = oneModelFunc (NULL, onePar, x);
    320308        } else {
    321             value = mFunc (oneDeriv, onePar, x);
    322             for (int n = 1; n < nPar; n++) {
     309            value = oneModelFunc (oneDeriv, onePar, x);
     310            for (int n = 1; n < nOnePar; n++) {
    323311                dpars[nOff + n] = dPAR[n];
    324312            }
    … …  
    333321
    334322/*
    335 i:         0           1               2
     323i:       0                   1                 2
    336324n:         1  2  3  4  5  6  1  2  3  4  5  6  1  2  3  4  5  6
    337325i*6 + n: 0 1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18
    … …  
    342330                     pmSourceFitMode mode)
    343331{
    344     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
     332    psTrace("psModules.objects", 3, "---- %s begin ----\n", __func__);
    345333    PS_ASSERT_PTR_NON_NULL(source, false);
    346     PS_ASSERT_PTR_NON_NULL(source->moments, false);
    347334    PS_ASSERT_PTR_NON_NULL(source->pixels, false);
    348335    PS_ASSERT_PTR_NON_NULL(source->mask, false);
    … …  
    353340    psBool rc        = true;
    354341
    355     // base values on first model
     342    // maximum number of valid pixels
     343    psS32 nPix = source->pixels->numRows * source->pixels->numCols;
     344
     345    // construct the coordinate and value entries
     346    psArray *x = psArrayAllocEmpty(nPix);
     347    psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
     348    psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
     349
     350    // fill in the coordinate and value entries
     351    nPix = 0;
     352    for (psS32 i = 0; i < source->pixels->numRows; i++) {
     353        for (psS32 j = 0; j < source->pixels->numCols; j++) {
     354            // skip masked points
     355            if (source->mask->data.U8[i][j]) {
     356                continue;
     357            }
     358            // skip zero-weight points
     359            if (source->weight->data.F32[i][j] == 0) {
     360                continue;
     361            }
     362            psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     363
     364            // Convert i/j to image space:
     365            coord->data.F32[0] = (psF32) (j + source->pixels->col0);
     366            coord->data.F32[1] = (psF32) (i + source->pixels->row0);
     367            x->data[nPix] = (psPtr *) coord;
     368            y->data.F32[nPix] = source->pixels->data.F32[i][j];
     369
     370            // psMinimizeLMChi2 takes wt = 1/dY^2.  suggestion from RHL is to use the local sky
     371            // as weight to avoid the bias from systematic errors here we would just use the
     372            // source sky variance
     373            if (PM_SOURCE_FIT_MODEL_PIX_WEIGHTS) {
     374                yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
     375            } else {
     376                yErr->data.F32[nPix] = 1.0 / PM_SOURCE_FIT_MODEL_WEIGHT;
     377            }
     378            nPix++;
     379        }
     380    }
     381    x->n = nPix;
     382    y->n = nPix;
     383    yErr->n = nPix;
     384
     385    // base values on first model (all models must be identical)
    356386    pmModel *model = modelSet->data[0];
    357387
    358     // set the static variables
    359     pmModelFitSetInit (model->type);
    360 
     388    // determine number of model parameters
    361389    int nSrc = modelSet->n;
    362390    int nPar = model->params->n - 1;  // number of object parameters (excluding sky)
    … …  
    366394    psVector *dparams = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    367395
    368     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
    369 
    370     // define limits for a single-source model
    371     psVector *oneDelta;
    372     psVector *oneParMin;
    373     psVector *oneParMax;
    374     modelLimits (&oneDelta, &oneParMin, &oneParMax);
    375 
    376     psMinConstrain *constrain = psMinConstrainAlloc();
    377     constrain->paramMin = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    378     constrain->paramMax = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    379     constrain->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
    380     constrain->paramDelta = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_F32);
    381 
    382     // set the parameter guesses and constraints for the multiple models
    383     // first the values for the single sky parameter
     396    // set the static variables
     397    pmSourceFitSetInit (model->type);
     398
     399    // create the minimization constraints
     400    psMinConstraint *constraint = psMinConstraintAlloc();
     401    constraint->paramMask = psVectorAlloc (nSrc*nPar + 1, PS_TYPE_U8);
     402    constraint->checkLimits = pmSourceFitSet_CheckLimits;
     403
     404    // set the parameter guesses for the multiple models
     405    // first the value for the single sky parameter
    384406    params->data.F32[0] = model->params->data.F32[0];
    385     constrain->paramMin->data.F32[0]   = oneParMin->data.F32[0];
    386     constrain->paramMax->data.F32[0]   = oneParMax->data.F32[0];
    387     constrain->paramDelta->data.F32[0] = oneDelta->data.F32[0];
    388407
    389408    // next, the values for the source parameters
    … …  
    392411        for (int n = 1; n < nPar + 1; n++) {
    393412            params->data.F32[i*nPar + n] = model->params->data.F32[n];
    394             constrain->paramDelta->data.F32[i*nPar + n] = oneDelta->data.F32[n];
    395             constrain->paramMin->data.F32[i*nPar + n]   = oneParMin->data.F32[n];
    396             constrain->paramMax->data.F32[i*nPar + n]   = oneParMax->data.F32[n];
    397         }
    398     }
    399     psFree (oneDelta);
    400     psFree (oneParMin);
    401     psFree (oneParMax);
     413        }
     414    }
    402415
    403416    if (psTraceGetLevel("psModules.objects") >= 5) {
    404417        for (int i = 0; i < params->n; i++) {
    405             fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constrain->paramMask->data.U8[i]);
     418            fprintf (stderr, "%d %f %d\n", i, params->data.F32[i], constraint->paramMask->data.U8[i]);
    406419        }
    407420    }
    … …  
    413426        // NORM-only model fits only source normalization (Io)
    414427        nParams = nSrc;
    415         psVectorInit (constrain->paramMask, 1);
     428        psVectorInit (constraint->paramMask, 1);
    416429        for (int i = 0; i < nSrc; i++) {
    417             constrain->paramMask->data.U8[1 + i*nPar] = 0;
     430            constraint->paramMask->data.U8[1 + i*nPar] = 0;
    418431        }
    419432        break;
    … …  
    421434        // PSF model only fits x,y,Io
    422435        nParams = 3*nSrc;
    423         psVectorInit (constrain->paramMask, 1);
     436        psVectorInit (constraint->paramMask, 1);
    424437        for (int i = 0; i < nSrc; i++) {
    425             constrain->paramMask->data.U8[1 + i*nPar] = 0;
    426             constrain->paramMask->data.U8[2 + i*nPar] = 0;
    427             constrain->paramMask->data.U8[3 + i*nPar] = 0;
     438            constraint->paramMask->data.U8[1 + i*nPar] = 0;
     439            constraint->paramMask->data.U8[2 + i*nPar] = 0;
     440            constraint->paramMask->data.U8[3 + i*nPar] = 0;
    428441        }
    429442        break;
    … …  
    431444        // EXT model fits all params (except sky)
    432445        nParams = nPar*nSrc;
    433         psVectorInit (constrain->paramMask, 0);
    434         constrain->paramMask->data.U8[0] = 1;
    435         break;
    436     case PM_SOURCE_FIT_PSF_AND_SKY:
    437         nParams = 1 + 3*nSrc;
    438         psAbort ("pmSourceFitModel", "PSF + SKY not currently available");
    439         break;
    440     case PM_SOURCE_FIT_EXT_AND_SKY:
    441         nParams = 1 + nPar*nSrc;
    442         psAbort ("pmSourceFitModel", "EXT + SKY not currently available");
     446        psVectorInit (constraint->paramMask, 0);
     447        constraint->paramMask->data.U8[0] = 1;
    443448        break;
    444449    default:
    … …  
    448453    // force the floating parameters to fall within the contraint ranges
    449454    for (int i = 0; i < params->n; i++) {
    450         if (constrain->paramMask->data.U8[i])
    451             continue;
    452         params->data.F32[i] = PS_MIN(constrain->paramMax->data.F32[i], params->data.F32[i]);
    453         params->data.F32[i] = PS_MAX(constrain->paramMin->data.F32[i], params->data.F32[i]);
    454     }
    455 
    456     // maximum number of valid pixels
    457     psS32 nPix = source->pixels->numRows * source->pixels->numCols;
    458 
    459     // local sky variance
    460     psF32 dSky = source->moments->dSky;
    461 
    462     // construct the coordinate and value entries
    463     psArray *x = psArrayAllocEmpty(nPix);
    464     psVector *y = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    465     psVector *yErr = psVectorAllocEmpty(nPix, PS_TYPE_F32);
    466 
    467     nPix = 0;
    468     for (psS32 i = 0; i < source->pixels->numRows; i++) {
    469         for (psS32 j = 0; j < source->pixels->numCols; j++) {
    470             // skip masked points
    471             if (source->mask->data.U8[i][j]) {
    472                 continue;
    473             }
    474             // skip zero-weight points
    475             if (source->weight->data.F32[i][j] == 0) {
    476                 continue;
    477             }
    478             psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
    479 
    480             // Convert i/j to image space:
    481             coord->data.F32[0] = (psF32) (j + source->pixels->col0);
    482             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
    483             x->data[nPix] = (psPtr *) coord;
    484             y->data.F32[nPix] = source->pixels->data.F32[i][j];
    485             // psMinimizeLMChi2 takes wt = 1/dY^2
    486             if (PM_PSF_POISSON_WEIGHTS) {
    487                 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j];
    488             } else {
    489                 yErr->data.F32[nPix] = 1.0 / dSky;
    490             }
    491             nPix++;
    492         }
    493     }
     455        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MIN, i, params->data.F32, NULL);
     456        pmSourceFitSet_CheckLimits (PS_MINIMIZE_PARAM_MAX, i, params->data.F32, NULL);
     457    }
     458
    494459    if (nPix <  nParams + 1) {
    495460        psTrace (__func__, 4, "insufficient valid pixels\n");
    … …  
    501466        psFree (params);
    502467        psFree (dparams);
    503         psFree(constrain->paramMask);
    504         psFree(constrain->paramMin);
    505         psFree(constrain->paramMax);
    506         psFree(constrain->paramDelta);
    507         psFree(constrain);
    508         pmModelFitSetClear ();
     468        psFree(constraint);
     469        pmSourceFitSetClear ();
    509470        return(false);
    510471    }
    511     x->n = nPix;
    512     y->n = nPix;
    513     yErr->n = nPix;
    514 
    515     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
    516                             PM_SOURCE_FIT_MODEL_TOLERANCE);
    517 
    518     psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    519 
    520     psTrace (__func__, 5, "fitting function\n");
    521     fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, pmModelFitSet);
     472
     473    psMinimization *myMin = psMinimizationAlloc (PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, PM_SOURCE_FIT_MODEL_TOLERANCE);
     474
     475    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     476
     477    fitStatus = psMinimizeLMChi2(myMin, covar, params, constraint, x, y, yErr, pmSourceFitSet_Function);
    522478    if (!fitStatus) {
    523         // psError(PS_ERR_UNKNOWN, false, "Failed to fit model (%d)\n", nSrc);
    524479        psTrace("psModules.objects", 4, "Failed to fit model (%d)\n", nSrc);
    525480    }
    … …  
    527482    // parameter errors from the covariance matrix
    528483    for (int i = 0; i < dparams->n; i++) {
    529         if ((constrain->paramMask != NULL) && constrain->paramMask->data.U8[i])
     484        if ((constraint->paramMask != NULL) && constraint->paramMask->data.U8[i])
    530485            continue;
    531         dparams->data.F32[i] = sqrt(covar->data.F64[i][i]);
     486        dparams->data.F32[i] = sqrt(covar->data.F32[i][i]);
    532487    }
    533488
    534489    // get the Gauss-Newton distance for fixed model parameters
    535     if (constrain->paramMask != NULL) {
    536         psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64);
     490    if (constraint->paramMask != NULL) {
     491        psVector *delta = psVectorAlloc (params->n, PS_TYPE_F32);
    537492        psVector *altmask = psVectorAlloc (params->n, PS_TYPE_U8);
    538493        altmask->data.U8[0] = 1;
    539494        for (int i = 1; i < dparams->n; i++) {
    540             altmask->data.U8[i] = (constrain->paramMask->data.U8[i]) ? 0 : 1;
    541         }
    542         psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmModelFitSet);
     495            altmask->data.U8[i] = (constraint->paramMask->data.U8[i]) ? 0 : 1;
     496        }
     497        psMinimizeGaussNewtonDelta(delta, params, altmask, x, y, yErr, pmSourceFitSet_Function);
    543498        for (int i = 0; i < dparams->n; i++) {
    544             if (!constrain->paramMask->data.U8[i])
     499            if (!constraint->paramMask->data.U8[i])
    545500                continue;
    546501            // note that delta is the value *subtracted* from the parameter
    547502            // to get the new guess.  for dparams to represent the direction
    548503            // of motion, we need to take -delta
    549             dparams->data.F32[i] = -delta->data.F64[i];
     504            dparams->data.F32[i] = -delta->data.F32[i];
    550505        }
    551506        psFree (delta);
    … …  
    558513        model->params->data.F32[0] = params->data.F32[0];
    559514        for (int n = 1; n < nPar + 1; n++) {
     515            if (psTraceGetLevel("psModules.objects") >= 4) {
     516                fprintf (stderr, "%f ", params->data.F32[i*nPar + n]);
     517            }
    560518            model->params->data.F32[n] = params->data.F32[i*nPar + n];
    561519            model->dparams->data.F32[n] = dparams->data.F32[i*nPar + n];
    562520        }
     521        psTrace ("psModules.objects", 4, " src %d", i);
     522
    563523        // save the resulting chisq, nDOF, nIter
    564524        // these are not unique for any one source
    … …  
    571531
    572532        // models can go insane: reject these
    573         onPic &= (model->params->data.F32[2] >= source->pixels->col0);
    574         onPic &= (model->params->data.F32[2] <  source->pixels->col0 + source->pixels->numCols);
    575         onPic &= (model->params->data.F32[3] >= source->pixels->row0);
    576         onPic &= (model->params->data.F32[3] <  source->pixels->row0 + source->pixels->numRows);
     533        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->col0);
     534        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->col0 + source->pixels->numCols);
     535        onPic &= (model->params->data.F32[PM_PAR_XPOS] >= source->pixels->row0);
     536        onPic &= (model->params->data.F32[PM_PAR_XPOS] <  source->pixels->row0 + source->pixels->numRows);
    577537        if (!onPic) {
    578538            model->status = PM_MODEL_OFFIMAGE;
    579539        }
    580540    }
     541    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    581542
    582543    source->mode |= PM_SOURCE_MODE_FITTED;
    … …  
    587548    psFree(myMin);
    588549    psFree(covar);
    589     psFree(constrain->paramMask);
    590     psFree(constrain->paramMin);
    591     psFree(constrain->paramMax);
    592     psFree(constrain->paramDelta);
    593     psFree(constrain);
     550    psFree(constraint);
    594551    psFree(params);
    595552    psFree(dparams);
    596553
    597     // free static memory used by pmModelFitSet
    598     pmModelFitSetClear ();
     554    // free static memory used by pmSourceFitSet
     555    pmSourceFitSetClear ();
    599556
    600557    rc = (onPic && fitStatus);
    601558    psTrace (__func__, 5, "onPic: %d, fitStatus: %d, nIter: %d, chisq: %f, nDof: %d\n", onPic, fitStatus, model->nIter, model->chisq, model->nDOF);
    602     psTrace("psModules.objects", 3, "---- %s end : status %d ----\n", __func__, rc);
     559    psTrace("psModules.objects", 5, "---- %s end (%d) ----\n", __func__, rc);
    603560    return(rc);
    604561}
Note: See TracChangeset for help on using the changeset viewer.