Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.c (revision 30784) @@ -67,4 +67,5 @@ tmp->chisqNorm = NAN; tmp->nDOF = 0; + tmp->nPar = 0; tmp->nPix = 0; tmp->nIter = 0; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.h (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmModel.h (revision 30784) @@ -39,5 +39,6 @@ float magErr; ///< integrated model magnitude error int nPix; ///< number of pixels used for fit - int nDOF; ///< number of degrees of freedom + int nPar; ///< number of parameters in fit + int nDOF; ///< number of degrees of freedom (nDOF = nPix - nPar) int nIter; ///< number of iterations to reach min pmModelStatus flags; ///< model status flags Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.c (revision 30784) @@ -254,5 +254,6 @@ pcm->nPix = nPix; - pcm->nDOF = nPix - nParams - 1; + pcm->nPar = nParams; + pcm->nDOF = nPix - nParams; return pcm; @@ -341,4 +342,5 @@ return false; } + pcm->nPar = nParams; pcm->nDOF = pcm->nPix - nParams; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.h (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmPCMdata.h (revision 30784) @@ -33,4 +33,5 @@ psMinConstraint *constraint; int nPix; + int nPar; int nDOF; } pmPCMdata; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitModel.c (revision 30784) @@ -236,16 +236,19 @@ model->covar = psMemIncrRefCounter(covar); } + model->nIter = myMin->iter; + model->nPar = nParams; + psTrace ("psModules.objects", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value); // save the resulting chisq, nDOF, nIter + // NOTE: if (!options->poissonErrors) chisq will be wrong : recalculate if (options->poissonErrors) { model->chisq = myMin->value; model->nPix = y->n; - model->nDOF = y->n - nParams; + model->nDOF = y->n - model->nPar; model->chisqNorm = model->chisq / model->nDOF; } else { - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, nParams); - } - model->nIter = myMin->iter; + pmSourceChisqUnsubtracted (source, model, maskVal); + } // set the model success or failure status Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitPCM.c (revision 30784) @@ -84,4 +84,6 @@ } } + pcm->modelConv->nIter = myMin->iter; + pcm->modelConv->nPar = pcm->nPar; // save the resulting chisq, nDOF, nIter @@ -92,7 +94,7 @@ pcm->modelConv->chisqNorm = pcm->modelConv->chisq / pcm->modelConv->nDOF; } else { - pmSourceChisq (pcm->modelConv, source->pixels, source->maskObj, source->variance, maskVal, fitOptions->covarFactor, pcm->nPix - pcm->nDOF - 1); + // xxx this is wrong because it does not convolve with the psf + pmSourceChisqUnsubtracted (source, pcm->modelConv, maskVal); } - pcm->modelConv->nIter = myMin->iter; // set the model success or failure status Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceFitSet.c (revision 30784) @@ -335,4 +335,7 @@ psTrace ("psModules.objects", 4, " src %d", i); + model->nIter = myMin->iter; + // model->nPar is set by pmSourceFitSetMasks + // save the resulting chisq, nDOF, nIter // these are not unique for any one source @@ -340,10 +343,9 @@ model->chisq = myMin->value; model->nPix = nPix; - model->nDOF = nPix - model->params->n; + model->nDOF = nPix - model->nPar; model->chisqNorm = model->chisq / model->nDOF; } else { - pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, options->covarFactor, model->params->n); + pmSourceChisqUnsubtracted (source, model, maskVal); } - model->nIter = myMin->iter; // set the model success or failure status @@ -399,4 +401,5 @@ for (int i = 0; i < set->paramSet->n; i++) { psVector *paramOne = set->paramSet->data[i]; + pmModel *modelOne = set->modelSet->data[i]; switch (mode) { @@ -406,4 +409,5 @@ if (j == PM_PAR_I0) continue; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1; + modelOne->nPar = 1; } break; @@ -415,4 +419,5 @@ if (j == PM_PAR_I0) continue; constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + j] = 1; + modelOne->nPar = 3; } break; @@ -420,4 +425,5 @@ // EXT model fits all params (except sky) constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n + PM_PAR_SKY] = 1; + modelOne->nPar = paramOne->n - 1; break; default: Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourceOutputs.c (revision 30784) @@ -72,5 +72,5 @@ } - *nImageOverlap = 1; + *nImageOverlap = psMetadataLookupS32 (&status2, header, "NINPUTS"); return true; Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.c (revision 30784) @@ -741,6 +741,6 @@ # endif -// determine chisq, etc for linear normalization-only fit -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal, const float covarFactor, int nParams) +// determine chisq, nPix, nDOF, chisqNorm : model->nPar must be set +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *variance, psImageMaskType maskVal) { PS_ASSERT_PTR_NON_NULL(model, false); @@ -757,12 +757,10 @@ if (variance->data.F32[j][i] <= 0) continue; - // dC += PS_SQR (image->data.F32[j][i]) / (covarFactor * variance->data.F32[j][i]); dC += PS_SQR (image->data.F32[j][i]) / variance->data.F32[j][i]; Npix ++; } } - model->nPix = Npix; - model->nDOF = Npix - nParams - 1; + model->nDOF = Npix - model->nPar; model->chisq = dC; model->chisqNorm = dC / model->nDOF; @@ -771,4 +769,54 @@ } + +// return source aperture magnitude +bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal) +{ + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(model, false); + + float dC = 0.0; + int Npix = 0; + + // the model function returns the source flux at a position + psVector *coord = psVectorAlloc(2, PS_TYPE_F32); + + psVector *params = model->params; + psImage *image = source->pixels; + psImage *mask = source->maskObj; + psImage *variance = source->variance; + + int dX = image->col0; + int dY = image->row0; + + for (int iy = 0; iy < image->numRows; iy++) { + for (int ix = 0; ix < image->numCols; ix++) { + + // skip pixels which are masked + if (mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal) continue; + + if (variance->data.F32[iy][ix] <= 0) continue; + + coord->data.F32[0] = (psF32) ix + dX + 0.5; + coord->data.F32[1] = (psF32) iy + dY + 0.5; + + // for the full model, add all points + float value = model->modelFunc (NULL, params, coord); + + // fprintf (stderr, "%d, %d : %f, %f : %f - %f : %f\n", + // ix, iy, coord->data.F32[0], coord->data.F32[1], image->data.F32[iy][ix], value, dC); + + dC += PS_SQR (image->data.F32[iy][ix] - value) / variance->data.F32[iy][ix]; + Npix ++; + } + } + model->nPix = Npix; + model->nDOF = Npix - model->nPar; + model->chisq = dC; + model->chisqNorm = dC / model->nDOF; + + psFree (coord); + return (true); +} double pmSourceModelWeight(const pmSource *Mi, int term, const bool unweighted_sum, const float covarFactor, psImageMaskType maskVal) Index: /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h =================================================================== --- /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h (revision 30783) +++ /branches/eam_branches/ipp-20110213/psModules/src/objects/pmSourcePhotometry.h (revision 30784) @@ -69,5 +69,6 @@ bool pmSourcePixelWeight (pmSource *source, pmModel *model, psImage *mask, psImageMaskType maskVal, float radius); -bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal, const float covarFactor, int nParams); +bool pmSourceChisq (pmModel *model, psImage *image, psImage *mask, psImage *weight, psImageMaskType maskVal); +bool pmSourceChisqUnsubtracted (pmSource *source, pmModel *model, psImageMaskType maskVal); bool pmSourceMeasureDiffStats (pmSource *source, psImageMaskType maskVal, psImageMaskType markVal); Index: /branches/eam_branches/ipp-20110213/psphot/doc/stack.txt =================================================================== --- /branches/eam_branches/ipp-20110213/psphot/doc/stack.txt (revision 30783) +++ /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 30783) +++ /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 30783) +++ /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 30783) +++ /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 30783) +++ /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; }