Index: branches/eam_branches/ipp-20110213/psphot/doc/stack.txt =================================================================== --- branches/eam_branches/ipp-20110213/psphot/doc/stack.txt (revision 30780) +++ branches/eam_branches/ipp-20110213/psphot/doc/stack.txt (revision 30784) @@ -49,11 +49,23 @@ 1652559 2000000 pmSubtractionMatch.c:189 - - 2132413 4000000 pmFPAfileDefine.c:1275 + chisq: 2133004 4000000 pmFPACopy.c:71 2133010 4000000 pmFPACopy.c:71 2133007 2000000 pmFPACopy.c:71 + + backmdl: + 2132413 4000000 pmFPAfileDefine.c:1275 + + ??? 8614548 2000000 psBinaryOp.c:502 8617313 2000000 psBinaryOp.c:502 + + pmSource (modelFlx): + 2775 * 6k = 16.9M + + pmSource (maskObj): + 2775 * 3k = 8.4M + + (span + footprints account for ~5 - 10 M, rest in little things) 20101221 Index: branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinear.c =================================================================== --- branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinear.c (revision 30780) +++ branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinear.c (revision 30784) @@ -312,7 +312,9 @@ for (int i = 0; i < fitSources->n; i++) { pmSource *source = fitSources->data[i]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor, 1); + if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) { + model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value + } + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); Index: branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinearStack.c =================================================================== --- branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinearStack.c (revision 30780) +++ branches/eam_branches/ipp-20110213/psphot/src/psphotFitSourcesLinearStack.c (revision 30784) @@ -163,7 +163,9 @@ for (int i = 0; i < fitSources->n; i++) { pmSource *source = fitSources->data[i]; - if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue; pmModel *model = pmSourceGetModel (NULL, source); - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR, 1); + if (!(source->mode & PM_SOURCE_MODE_NONLINEAR_FIT)) { + model->nPar = 1; // LINEAR-only sources have 1 parameter; NONLINEAR sources have their original value + } + pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal); } psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem); Index: branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c =================================================================== --- branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c (revision 30780) +++ branches/eam_branches/ipp-20110213/psphot/src/psphotReadout.c (revision 30784) @@ -9,4 +9,53 @@ } +// for now, let's store the detections on the readout->analysis for each readout +bool psphotDumpChisqs (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + static int npass = 0; + char filename[64]; + + bool status = true; + + int num = psphotFileruleCount(config, filerule); + + snprintf (filename, 64, "chisq.%02d.dat", npass); + FILE *f = fopen (filename, "w"); + + // loop over the available readouts + for (int i = 0; i < num; i++) { + + // find the currently selected readout + pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest + psAssert (file, "missing file?"); + + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); + psAssert (readout, "missing readout?"); + + pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS"); + psAssert (detections, "missing detections?"); + + psArray *sources = detections->allSources; + psAssert (sources, "missing sources?"); + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source) continue; + + pmModel *model = pmSourceGetModel (NULL, source); + if (!model) continue; + + if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) { + 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); + } else { + 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); + } + } + } + fclose (f); + npass ++; + + return true; +} + bool psphotReadout(pmConfig *config, const pmFPAview *view, const char *filerule) { @@ -147,4 +196,5 @@ // linear PSF fit to source peaks, subtract the models from the image (in PSF mask) psphotFitSourcesLinear (config, view, filerule, false); // pass 1 (detections->allSources) + psphotDumpChisqs (config, view, filerule); // identify CRs and extended sources (only unmeasured sources are measured) @@ -157,4 +207,5 @@ // replace model flux, adjust mask as needed, fit, subtract the models (full stamp) psphotBlendFit (config, view, filerule); // pass 1 (detections->allSources) + psphotDumpChisqs (config, view, filerule); // replace all sources @@ -164,4 +215,5 @@ // NOTE : apply to ALL sources (extended + psf) psphotFitSourcesLinear (config, view, filerule, true); // pass 2 (detections->allSources) + psphotDumpChisqs (config, view, filerule); // if we only do one pass, skip to extended source analysis @@ -209,4 +261,5 @@ // NOTE: apply to ALL sources psphotFitSourcesLinear (config, view, filerule, true); // pass 3 (detections->allSources) + psphotDumpChisqs (config, view, filerule); } Index: branches/eam_branches/ipp-20110213/psphot/src/psphotSourceFits.c =================================================================== --- branches/eam_branches/ipp-20110213/psphot/src/psphotSourceFits.c (revision 30780) +++ branches/eam_branches/ipp-20110213/psphot/src/psphotSourceFits.c (revision 30784) @@ -101,8 +101,4 @@ if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - // correct model chisq for flux trend - double chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]); - PSF->chisqNorm = PSF->chisq / chiTrend; - // evaluate the blend objects, subtract if good, free otherwise for (int i = 1; i < modelSet->n; i++) { @@ -111,8 +107,4 @@ if (!isfinite(model->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - - // correct model chisq for flux trend - chiTrend = psPolynomial1DEval (psf->ChiTrend, model->params->data.F32[PM_PAR_I0]); - model->chisqNorm = model->chisq / chiTrend; // if this one failed, skip it @@ -159,5 +151,4 @@ bool psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal) { - double chiTrend; pmSourceFitOptions options = *fitOptions; @@ -182,8 +173,4 @@ // clear the circular mask psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // correct model chisq for flux trend - chiTrend = psPolynomial1DEval (psf->ChiTrend, PSF->params->data.F32[PM_PAR_I0]); - PSF->chisqNorm = PSF->chisq / chiTrend; // does the PSF model succeed? @@ -225,5 +212,4 @@ bool okEXT, okDBL; float chiEXT, chiDBL; - double chiTrend; pmModel *ONE = NULL; pmSource *tmpSrc = NULL; @@ -271,11 +257,8 @@ // correct first model chisqs for flux trend - chiDBL = NAN; ONE = DBL->data[0]; if (ONE) { if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); - ONE->chisqNorm = ONE->chisq / chiTrend; - chiDBL = ONE->chisq / ONE->nDOF; // save chisq for double-star/galaxy comparison + chiDBL = ONE->chisqNorm; // save chisq for double-star/galaxy comparison ONE->fitRadius = radius; } @@ -285,6 +268,4 @@ if (ONE) { if (!isfinite(ONE->params->data.F32[PM_PAR_I0])) psAbort("nan in fit"); - chiTrend = psPolynomial1DEval (psf->ChiTrend, ONE->params->data.F32[1]); - ONE->chisqNorm = ONE->chisq / chiTrend; ONE->fitRadius = radius; } @@ -298,5 +279,5 @@ okEXT = psphotEvalEXT (tmpSrc, EXT); - chiEXT = EXT ? EXT->chisq / EXT->nDOF : NAN; + chiEXT = EXT ? EXT->chisqNorm : NAN; }