IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 30784


Ignore:
Timestamp:
Mar 3, 2011, 3:11:09 PM (15 years ago)
Author:
eugene
Message:

fix the calculation of chisq values for constant errors (re-calculated, do not attempt to correct raw value)

Location:
branches/eam_branches/ipp-20110213
Files:
15 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.c

    r30705 r30784  
    6767    tmp->chisqNorm = NAN;
    6868    tmp->nDOF  = 0;
     69    tmp->nPar  = 0;
    6970    tmp->nPix  = 0;
    7071    tmp->nIter = 0;

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.h

    r30705 r30784  
    3939    float magErr;                       ///< integrated model magnitude error
    4040    int nPix;                           ///< number of pixels used for fit
    41     int nDOF;                           ///< number of degrees of freedom
     41    int nPar;                           ///< number of parameters in fit
     42    int nDOF;                           ///< number of degrees of freedom (nDOF = nPix - nPar)
    4243    int nIter;                          ///< number of iterations to reach min
    4344    pmModelStatus flags;                ///< model status flags

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c

    r30763 r30784  
    254254
    255255    pcm->nPix = nPix;
    256     pcm->nDOF = nPix - nParams - 1;
     256    pcm->nPar = nParams;
     257    pcm->nDOF = nPix - nParams;
    257258
    258259    return pcm;
    … …  
    341342        return false;
    342343    }
     344    pcm->nPar = nParams;
    343345    pcm->nDOF = pcm->nPix - nParams;
    344346

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.h

    r30621 r30784  
    3333    psMinConstraint *constraint;
    3434    int nPix;
     35    int nPar;
    3536    int nDOF;
    3637} pmPCMdata;

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c

    r30763 r30784  
    236236        model->covar = psMemIncrRefCounter(covar);
    237237    }
     238    model->nIter = myMin->iter;
     239    model->nPar = nParams;
     240
    238241    psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
    239242
    240243    // save the resulting chisq, nDOF, nIter
     244    // NOTE: if (!options->poissonErrors) chisq will be wrong : recalculate
    241245    if (options->poissonErrors) {
    242246        model->chisq = myMin->value;
    243247        model->nPix  = y->n;
    244         model->nDOF  = y->n - nParams;
     248        model->nDOF  = y->n - model->nPar;
    245249        model->chisqNorm = model->chisq / model->nDOF;
    246250    } else {
    247         pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, nParams);
    248     }
    249     model->nIter = myMin->iter;
     251        pmSourceChisqUnsubtracted (source, model, maskVal);
     252    }
    250253
    251254    // set the model success or failure status

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c

    r30763 r30784  
    8484        }
    8585    }
     86    pcm->modelConv->nIter = myMin->iter;
     87    pcm->modelConv->nPar = pcm->nPar;
    8688
    8789    // save the resulting chisq, nDOF, nIter
    … …  
    9294        pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF;
    9395    } else {
    94         pmSourceChisq (pcm->modelConv, source->pixels, source->maskObj, source->variance, maskVal, fitOptions->covarFactor, pcm->nPix - pcm->nDOF - 1);
     96        // xxx this is wrong because it does not convolve with the psf
     97        pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal);
    9598    }
    96     pcm->modelConv->nIter = myMin->iter;
    9799
    98100    // set the model success or failure status

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c

    r30763 r30784  
    335335        psTrace ("psModules.objects", 4, " src %d", i);
    336336
     337        model->nIter = myMin->iter;
     338        // model->nPar is set by pmSourceFitSetMasks
     339
    337340        // save the resulting chisq, nDOF, nIter
    338341        // these are not unique for any one source
    … …  
    340343            model->chisq = myMin->value;
    341344            model->nPix  = nPix;
    342             model->nDOF  = nPix - model->params->n;
     345            model->nDOF  = nPix - model->nPar;
    343346            model->chisqNorm = model->chisq / model->nDOF;
    344347        } else {
    345             pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, model->params->n);
     348            pmSourceChisqUnsubtracted (source, model, maskVal);
    346349        }
    347         model->nIter = myMin->iter;
    348350
    349351        // set the model success or failure status
    … …  
    399401    for (int i = 0; i < set->paramSet->n; i++) {
    400402        psVector *paramOne = set->paramSet->data[i];
     403        pmModel  *modelOne = set->modelSet->data[i];
    401404
    402405        switch (mode) {
    … …  
    406409                if (j == PM_PAR_I0) continue;
    407410                constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1;
     411                modelOne->nPar = 1;
    408412            }
    409413            break;
    … …  
    415419                if (j == PM_PAR_I0) continue;
    416420                constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1;
     421                modelOne->nPar = 3;
    417422            }
    418423            break;
    … …  
    420425            // EXT model fits all params (except sky)
    421426            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + PM_PAR_SKY] = 1;
     427            modelOne->nPar = paramOne->n - 1;
    422428            break;
    423429          default:

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c

    r30763 r30784  
    7272    }
    7373
    74     *nImageOverlap = 1;
     74    *nImageOverlap = psMetadataLookupS32 (&status2, header, "NINPUTS");
    7575    return true;
    7676

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c

    r30772 r30784  
    741741# endif
    742742
    743 // determine chisq, etc for linear normalization-only fit
    744 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor, int nParams)
     743// determine chisq, nPix, nDOF, chisqNorm : model->nPar must be set
     744bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal)
    745745{
    746746    PS_ASSERT_PTR_NON_NULL(model, false);
    … …  
    757757            if (variance->data.F32[j][i] <= 0)
    758758                continue;
    759             // dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]);
    760759            dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i];
    761760            Npix ++;
    762761        }
    763762    }
    764 
    765763    model->nPix = Npix;
    766     model->nDOF = Npix - nParams - 1;
     764    model->nDOF = Npix - model->nPar;
    767765    model->chisq = dC;
    768766    model->chisqNorm = dC / model->nDOF;
    … …  
    771769}
    772770
     771
     772// return source aperture magnitude
     773bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal)
     774{
     775    PS_ASSERT_PTR_NON_NULL(source, false);
     776    PS_ASSERT_PTR_NON_NULL(model, false);
     777
     778    float dC = 0.0;
     779    int Npix = 0;
     780
     781    // the model function returns the source flux at a position
     782    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
     783
     784    psVector *params = model->params;
     785    psImage  *image = source->pixels;
     786    psImage  *mask = source->maskObj;
     787    psImage  *variance = source->variance;
     788
     789    int dX = image->col0;
     790    int dY = image->row0;
     791
     792    for (int iy = 0; iy < image->numRows; iy++) {
     793        for (int ix = 0; ix < image->numCols; ix++) {
     794
     795            // skip pixels which are masked
     796            if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue;
     797
     798            if (variance->data.F32[iy][ix] <= 0) continue;
     799
     800            coord->data.F32[0] = (psF32) ix + dX + 0.5;
     801            coord->data.F32[1] = (psF32) iy + dY + 0.5;
     802
     803            // for the full model, add all points
     804            float value = model->modelFunc (NULL, params, coord);
     805
     806            // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n",
     807            // ix, iy, coord->data.F32[0], coord->data.F32[1], image->data.F32[iy][ix], value, dC);
     808
     809            dC += PS_SQR (image->data.F32[iy][ix] - value) / variance->data.F32[iy][ix];
     810            Npix ++;
     811        }
     812    }
     813    model->nPix = Npix;
     814    model->nDOF = Npix - model->nPar;
     815    model->chisq = dC;
     816    model->chisqNorm = dC / model->nDOF;
     817
     818    psFree (coord);
     819    return (true);
     820}
    773821
    774822double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal)

  • branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h

    r30772 r30784  
    6969bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius);
    7070
    71 bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor, int nParams);
     71bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal);
     72bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal);
    7273
    7374bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal);

  • branches/eam_branches/ipp-20110213/psphot/doc/stack.txt

    r30779 r30784  
    4949     1652559         2000000    pmSubtractionMatch.c:189
    5050
    51  
    52      2132413         4000000    pmFPAfileDefine.c:1275
     51  chisq:
    5352     2133004         4000000    pmFPACopy.c:71
    5453     2133010         4000000    pmFPACopy.c:71
    5554     2133007         2000000    pmFPACopy.c:71
     55
     56  backmdl: 
     57     2132413         4000000    pmFPAfileDefine.c:1275
     58
     59  ???
    5660     8614548         2000000    psBinaryOp.c:502
    5761     8617313         2000000    psBinaryOp.c:502
     62
     63  pmSource (modelFlx):
     64     2775 * 6k = 16.9M
     65
     66  pmSource (maskObj):
     67     2775 * 3k =  8.4M
     68
     69  (span + footprints account for ~5 - 10 M, rest in little things)
    5870
    597120101221

  • branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinear.c

    r30764 r30784  
    312312    for (int i = 0; i < fitSources->n; i++) {
    313313        pmSource *source = fitSources->data[i];
    314         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
    315314        pmModel *model = pmSourceGetModel (NULL, source);
    316         pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor, 1);
     315        if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
     316            model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
     317        }
     318        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
    317319    }
    318320    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);

  • branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinearStack.c

    r30764 r30784  
    163163    for (int i = 0; i < fitSources->n; i++) {
    164164        pmSource *source = fitSources->data[i];
    165         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
    166165        pmModel *model = pmSourceGetModel (NULL, source);
    167         pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR, 1);
     166        if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) {
     167            model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value
     168        }
     169        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
    168170    }
    169171    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);

  • branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c

    r30749 r30784  
    99}
    1010
     11// for now, let's store the detections on the readout->analysis for each readout
     12bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule)
     13{
     14    static int npass = 0;
     15    char filename[64];
     16
     17    bool status = true;
     18
     19    int num = psphotFileruleCount(config, filerule);
     20
     21    snprintf (filename, 64, "chisq.%02d.dat", npass);
     22    FILE *f = fopen (filename, "w");
     23
     24    // loop over the available readouts
     25    for (int i = 0; i < num; i++) {
     26
     27        // find the currently selected readout
     28        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
     29        psAssert (file, "missing file?");
     30
     31        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
     32        psAssert (readout, "missing readout?");
     33
     34        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
     35        psAssert (detections, "missing detections?");
     36
     37        psArray *sources = detections->allSources;
     38        psAssert (sources, "missing sources?");
     39
     40        for (int i = 0; i < sources->n; i++) {
     41            pmSource *source = sources->data[i];
     42            if (!source) continue;
     43
     44            pmModel *model = pmSourceGetModel (NULL, source);
     45            if (!model) continue;
     46       
     47            if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) {
     48                fprintf (f, "%f %f %f %d %d %f  1 NONLINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
     49            } else {
     50                fprintf (f, "%f %f %f %d %d %f  0 LINEAR\n", model->mag, model->params->data.F32[1], model->chisq, model->nDOF, model->nPix, model->chisqNorm);
     51            }
     52        }
     53    }
     54    fclose (f);
     55    npass ++;
     56
     57    return true;
     58}
     59
    1160bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) {
    1261
    … …  
    147196    // linear PSF fit to source peaks, subtract the models from the image (in PSF mask)
    148197    psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources)
     198    psphotDumpChisqs (config, view, filerule);
    149199
    150200    // identify CRs and extended sources (only unmeasured sources are measured)
    … …  
    157207    // replace model flux, adjust mask as needed, fit, subtract the models (full stamp)
    158208    psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources)
     209    psphotDumpChisqs (config, view, filerule);
    159210
    160211    // replace all sources
    … …  
    164215    // NOTE : apply to ALL sources (extended + psf)
    165216    psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources)
     217    psphotDumpChisqs (config, view, filerule);
    166218
    167219    // if we only do one pass, skip to extended source analysis
    … …  
    209261        // NOTE: apply to ALL sources
    210262        psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources)
     263        psphotDumpChisqs (config, view, filerule);
    211264    }
    212265

  • branches/eam_branches/ipp-20110213/psphot/src/psphotSourceFits.c

    r30624 r30784  
    101101    if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    102102
    103     // correct model chisq for flux trend
    104     double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]);
    105     PSF->chisqNorm = PSF->chisq / chiTrend;
    106 
    107103    // evaluate the blend objects, subtract if good, free otherwise
    108104    for (int i = 1; i < modelSet->n; i++) {
    … …  
    111107
    112108        if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    113 
    114         // correct model chisq for flux trend
    115         chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[PM_PAR_I0]);
    116         model->chisqNorm = model->chisq / chiTrend;
    117109
    118110        // if this one failed, skip it
    … …  
    159151bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) {
    160152
    161     double chiTrend;
    162153    pmSourceFitOptions options = *fitOptions;
    163154
    … …  
    182173    // clear the circular mask
    183174    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
    184 
    185     // correct model chisq for flux trend
    186     chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]);
    187     PSF->chisqNorm = PSF->chisq / chiTrend;
    188175
    189176    // does the PSF model succeed?
    … …  
    225212    bool okEXT, okDBL;
    226213    float chiEXT, chiDBL;
    227     double chiTrend;
    228214    pmModel *ONE = NULL;
    229215    pmSource *tmpSrc = NULL;
    … …  
    271257
    272258        // correct first model chisqs for flux trend
    273         chiDBL = NAN;
    274259        ONE = DBL->data[0];
    275260        if (ONE) {
    276261            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    277             chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    278             ONE->chisqNorm = ONE->chisq / chiTrend;
    279             chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison
     262            chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison
    280263            ONE->fitRadius = radius;
    281264        }
    … …  
    285268        if (ONE) {
    286269            if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
    287             chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]);
    288             ONE->chisqNorm = ONE->chisq / chiTrend;
    289270            ONE->fitRadius = radius;
    290271        }
    … …  
    298279
    299280        okEXT = psphotEvalEXT (tmpSrc, EXT);
    300         chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN;
     281        chiEXT = EXT ? EXT->chisqNorm : NAN;
    301282    }
    302283
Note: See TracChangeset for help on using the changeset viewer.