Index: anches/eam_branches/20090715/psModules/src/objects/pmPSFtry.old.c =================================================================== --- /branches/eam_branches/20090715/psModules/src/objects/pmPSFtry.old.c (revision 25586) +++ (revision ) @@ -1,1142 +1,0 @@ -/** @file pmPSFtry.c - * - * XXX: need description of file purpose - * - * @author EAM, IfA - * - * @version $Revision: 1.69 $ $Name: not supported by cvs2svn $ - * @date $Date: 2009-01-27 06:39:38 $ - * Copyright 2004 Maui High Performance Computing Center, University of Hawaii - * - */ - -#ifdef HAVE_CONFIG_H -#include -#endif - -# include -#include "pmHDU.h" -#include "pmFPA.h" -#include "pmFPAMaskWeight.h" -#include "pmSpan.h" -#include "pmFootprint.h" -#include "pmPeaks.h" -#include "pmMoments.h" -#include "pmResiduals.h" -#include "pmGrowthCurve.h" -#include "pmTrend2D.h" -#include "pmPSF.h" -#include "pmModel.h" -#include "pmSource.h" -#include "pmSourceUtils.h" -#include "pmPSFtry.h" -#include "pmModelClass.h" -#include "pmModelUtils.h" -#include "pmSourceFitModel.h" -#include "pmSourcePhotometry.h" -#include "pmSourceVisual.h" - -// ******** pmPSFtry functions ************************************************** -// * pmPSFtry holds a single pmPSF model test, with the input sources, the freely -// * fitted version of the model, the pmPSF fit to the fitted model parameters, -// * and the PSF fits to the source. It also includes the statistics from the -// * fits, both the individual sources, and the collection - -// free a pmPSFtry structure -static void pmPSFtryFree (pmPSFtry *test) -{ - if (test == NULL) return; - - psFree (test->psf); - psFree (test->sources); - psFree (test->metric); - psFree (test->metricErr); - psFree (test->fitMag); - psFree (test->mask); - return; -} - -// allocate a pmPSFtry based on the desired sources and the model (identified by name) -pmPSFtry *pmPSFtryAlloc (const psArray *sources, const pmPSFOptions *options) -{ - pmPSFtry *test = (pmPSFtry *) psAlloc(sizeof(pmPSFtry)); - psMemSetDeallocator(test, (psFreeFunc) pmPSFtryFree); - - test->psf = pmPSFAlloc (options); - test->metric = psVectorAlloc (sources->n, PS_TYPE_F32); - test->metricErr = psVectorAlloc (sources->n, PS_TYPE_F32); - test->fitMag = psVectorAlloc (sources->n, PS_TYPE_F32); - test->mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK); - - psVectorInit (test->mask, 0); - psVectorInit (test->metric, 0.0); - psVectorInit (test->metricErr, 0.0); - psVectorInit (test->fitMag, 0.0); - - test->sources = psArrayAlloc (sources->n); - - for (int i = 0; i < sources->n; i++) { - test->sources->data[i] = pmSourceCopy (sources->data[i]); - } - - return (test); -} - -bool psMemCheckPSFtry(psPtr ptr) -{ - PS_ASSERT_PTR(ptr, false); - return ( psMemGetDeallocator(ptr) == (psFreeFunc) pmPSFtryFree); -} - - -// build a pmPSFtry for the given model: -// - fit each source with the free-floating model -// - construct the pmPSF from the collection of models -// - fit each source with the PSF-parameter models -// - measure the pmPSF quality metric (dApResid) - -// sources used in for pmPSFtry may be masked by the analysis -// mask values indicate the reason the source was rejected: - -// generate a pmPSFtry with a copy of the test PSF sources -pmPSFtry *pmPSFtryModel (const psArray *sources, const char *modelName, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) -{ - bool status; - int Next = 0; - int Npsf = 0; - - // validate the requested model name - options->type = pmModelClassGetType (modelName); - if (options->type == -1) { - psError (PS_ERR_UNKNOWN, true, "invalid model name %s", modelName); - return NULL; - } - - pmPSFtry *psfTry = pmPSFtryAlloc (sources, options); - if (psfTry == NULL) { - psError (PS_ERR_UNKNOWN, false, "failed to allocate psf model"); - return NULL; - } - - // maskVal is used to test for rejected pixels, and must include markVal - maskVal |= markVal; - - // stage 1: fit an EXT model to all candidates PSF sources -- this is independent of the modeled 2D variations in the PSF - psTimerStart ("psf.fit"); - for (int i = 0; i < psfTry->sources->n; i++) { - - pmSource *source = psfTry->sources->data[i]; - if (!source->moments) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - if (!source->moments->nPixels) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - continue; - } - - source->modelEXT = pmSourceModelGuess (source, psfTry->psf->type); - if (source->modelEXT == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : failed to generate model guess\n", i, source->peak->x, source->peak->y); - continue; - } - - // set object mask to define valid pixels - // XXX 0.5 PIX: is the circle symmetric about the peak coordinate (given 0.5,0.5 center)? - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); - - // fit model as EXT, not PSF - status = pmSourceFitModel (source, source->modelEXT, PM_SOURCE_FIT_EXT, maskVal); - - // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - - // exclude the poor fits - if (!status) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_EXT_FAIL; - psTrace ("psModules.objects", 4, "masking %d (%d,%d) : status is poor\n", i, source->peak->x, source->peak->y); - continue; - } - Next ++; - } - psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit ext: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Next, sources->n); - psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (EXT)\n", Next, sources->n); - - if (Next == 0) { - psError(PS_ERR_UNKNOWN, false, "No sources with good extended fits from which to determine PSF."); - psFree(psfTry); - return NULL; - } - - // stage 2: construct a psf (pmPSF) from this collection of model fits, including the 2D variation - if (!pmPSFFromPSFtry (psfTry)) { - psError(PS_ERR_UNKNOWN, false, "failed to construct a psf model from collection of sources"); - psFree(psfTry); - return NULL; - } - - // stage 3: refit with fixed shape parameters - psTimerStart ("psf.fit"); - for (int i = 0; i < psfTry->sources->n; i++) { - - pmSource *source = psfTry->sources->data[i]; - - // masked for: bad model fit, outlier in parameters - if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) { - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : source is masked\n", i, source->peak->x, source->peak->y); - continue; - } - - // set shape for this model based on PSF - source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf); - if (source->modelPSF == NULL) { - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL; - abort(); - continue; - } - source->modelPSF->radiusFit = options->radius; - - // XXXX use a different radius for the aperture magnitude than for the PSF fit? - - // set object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "OR", markVal); - - // fit the PSF model to the source - status = pmSourceFitModel (source, source->modelPSF, PM_SOURCE_FIT_PSF, maskVal); - - // skip poor fits - if (!status) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y); - continue; - } - - status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal); - if (!status || isnan(source->apMag)) { - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT; - psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y); - continue; - } - - // clear object mask to define valid pixels - psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->radius, "AND", PS_NOT_IMAGE_MASK(markVal)); - - psfTry->fitMag->data.F32[i] = source->psfMag; - psfTry->metric->data.F32[i] = source->apMag - source->psfMag; - psfTry->metricErr->data.F32[i] = source->errMag; - - psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n); - Npsf ++; - } - psfTry->psf->nPSFstars = Npsf; - - psLogMsg ("psphot.psftry", PS_LOG_MINUTIA, "fit psf: %f sec for %d of %ld sources\n", psTimerMark ("psf.fit"), Npsf, sources->n); - psTrace ("psModules.object", 3, "keeping %d of %ld PSF candidates (PSF)\n", Npsf, sources->n); - - if (Npsf == 0) { - psError(PS_ERR_UNKNOWN, false, "No sources with good PSF fits after model is built."); - psFree(psfTry); - return NULL; - } - - // measure the chi-square trend as a function of flux (PAR[PM_PAR_I0]) - // this should be linear for Poisson errors and quadratic for constant sky errors - psVector *flux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *chisq = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *mask = psVectorAlloc (psfTry->sources->n, PS_TYPE_VECTOR_MASK); - - // generate the x and y vectors, and mask missing models - for (int i = 0; i < psfTry->sources->n; i++) { - pmSource *source = psfTry->sources->data[i]; - if (source->modelPSF == NULL) { - flux->data.F32[i] = 0.0; - chisq->data.F32[i] = 0.0; - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0xff; - } else { - flux->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_I0]; - chisq->data.F32[i] = source->modelPSF->chisq / source->modelPSF->nDOF; - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = 0; - } - } - - // use 3hi/3lo sigma clipping on the chisq fit - psStats *stats = options->stats; - - // linear clipped fit of chisq trend vs flux - if (options->chiFluxTrend) { - bool result = psVectorClipFitPolynomial1D(psfTry->psf->ChiTrend, options->stats, - mask, 0xff, chisq, NULL, flux); - psStatsOptions meanStat = psStatsMeanOption(options->stats->options); // Statistic for mean - psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev - - psLogMsg ("pmPSFtry", 4, "chisq vs flux fit: %f +/- %f\n", - psStatsGetValue(stats, meanStat), psStatsGetValue(stats, stdevStat)); - - psFree(flux); - psFree(mask); - psFree(chisq); - - if (!result) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit psf->ChiTrend"); - psFree(psfTry); - return NULL; - } - } - - for (int i = 0; i < psfTry->psf->ChiTrend->nX + 1; i++) { - psLogMsg ("pmPSFtry", 4, "chisq vs flux fit term %d: %f +/- %f\n", i, - psfTry->psf->ChiTrend->coeff[i]*pow(10000, i), - psfTry->psf->ChiTrend->coeffErr[i]*pow(10000,i)); - } - - // XXX this function wants aperture radius for pmSourcePhotometry - if (!pmPSFtryMetric (psfTry, options)) { - psError(PS_ERR_UNKNOWN, false, "Attempt to fit PSF with model %s failed.", modelName); - psFree (psfTry); - return NULL; - } - - psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f : sky bias: %f\n", - modelName, psfTry->psf->ApResid, psfTry->psf->dApResid, psfTry->psf->skyBias); - - return (psfTry); -} - -bool pmPSFtryMetric (pmPSFtry *psfTry, pmPSFOptions *options) -{ - PS_ASSERT_PTR_NON_NULL(psfTry, false); - PS_ASSERT_PTR_NON_NULL(options, false); - PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); - - float RADIUS = options->radius; - - // the measured (aperture - fit) magnitudes (dA == psfTry->metric) - // depend on both the true ap-fit (dAo) and the bias in the sky measurement: - // dA = dAo + dsky/flux - // where flux is the flux of the star - // we fit this trend to find the infinite flux aperture correction (dAo), - // the nominal sky bias (dsky), and the error on dAo - // the values of dA are contaminated by stars with close neighbors in the aperture - // we use an outlier rejection to avoid this bias - - // r2rflux = radius^2 * ten(0.4*fitMag); - psVector *r2rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - - for (int i = 0; i < psfTry->sources->n; i++) { - if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) - continue; - r2rflux->data.F32[i] = PS_SQR(RADIUS) * pow(10.0, 0.4*psfTry->fitMag->data.F32[i]); - } - - // XXX test dump of aperture residual data - if (psTraceGetLevel("psModules.objects") >= 5) { - FILE *f = fopen ("apresid.dat", "w"); - for (int i = 0; i < psfTry->sources->n; i++) { - int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL); - - pmSource *source = psfTry->sources->data[i]; - float x = source->peak->x; - float y = source->peak->y; - - fprintf (f, "%d %d, %f %f %f %f %f %f \n", - i, keep, x, y, - psfTry->fitMag->data.F32[i], - r2rflux->data.F32[i], - psfTry->metric->data.F32[i], - psfTry->metricErr->data.F32[i]); - } - fclose (f); - } - - // This analysis of the apResid statistics is only approximate. The fitted magnitudes - // measured at this point (in the PSF fit) use Poisson errors, and are thus biased as a - // function of magnitude. We re-measure the apResid statistics later in psphot using the - // linear, constant-error fitting. Do not reject outliers with excessive vigor here. - - // fit ApTrend only to r2rflux, ignore x,y,flux variations for now - // linear clipped fit of ApResid to r2rflux - psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); - poly->coeffMask[1] = PS_POLY_MASK_SET; // fit only a constant offset (no SKYBIAS) - - // XXX replace this with a psVectorStats call? since we are not fitting the trend - bool result = psVectorClipFitPolynomial1D(poly, options->stats, psfTry->mask, PSFTRY_MASK_ALL, - psfTry->metric, psfTry->metricErr, r2rflux); - if (!result) { - psError(PS_ERR_UNKNOWN, false, "Failed to fit clipped poly"); - - psFree(poly); - psFree(r2rflux); - - return false; - } - psStatsOptions stdevStat = psStatsStdevOption(options->stats->options); // Statistic for stdev - psLogMsg ("pmPSFtryMetric", 4, "apresid: %f +/- %f; from statistics of %ld psf stars\n", poly->coeff[0], - psStatsGetValue(options->stats, stdevStat), psfTry->sources->n); - - float dSys = psVectorSystematicError (psfTry->metric, psfTry->metricErr, 0.1); - fprintf (stderr, "systematic error: %f\n", dSys); - - int n = 0; - psVector *bright = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); - psVector *brightErr = psVectorAllocEmpty (psfTry->metric->n, PS_TYPE_F32); - for (int i = 0; i < psfTry->metric->n; i++) { - if (!isfinite(psfTry->metric->data.F32[i])) continue; - if (!isfinite(psfTry->metricErr->data.F32[i])) continue; - if (psfTry->metricErr->data.F32[i] <= 0.0) continue; - if (psfTry->metricErr->data.F32[i] > 0.005) continue; - bright->data.F32[n] = psfTry->metric->data.F32[i]; - brightErr->data.F32[n] = psfTry->metricErr->data.F32[i]; - n++; - } - bright->n = brightErr->n = n; - - float dSysBright = psVectorSystematicError (bright, brightErr, 0.1); - fprintf (stderr, "bright systematic error: %f\n", dSysBright); - psFree(bright); - psFree(brightErr); - - // XXX test dump of fitted model (dump when tracing?) - if (psTraceGetLevel("psModules.objects") >= 4) { - FILE *f = fopen ("resid.dat", "w"); - psVector *apfit = psPolynomial1DEvalVector (poly, r2rflux); - for (int i = 0; i < psfTry->sources->n; i++) { - int keep = (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL); - - pmSource *source = psfTry->sources->data[i]; - float x = source->peak->x; - float y = source->peak->y; - - fprintf (f, "%d %d, %f %f %f %f %f %f %f\n", - i, keep, x, y, - psfTry->fitMag->data.F32[i], - r2rflux->data.F32[i], - psfTry->metric->data.F32[i], - psfTry->metricErr->data.F32[i], - apfit->data.F32[i]); - } - fclose (f); - psFree (apfit); - } - - // XXX drop the skyBias value, or include above?? - psfTry->psf->skyBias = poly->coeff[1]; - psfTry->psf->ApResid = poly->coeff[0]; - psfTry->psf->dApResid = psStatsGetValue(options->stats, stdevStat); - - psFree (r2rflux); - psFree (poly); - - return true; -} - -/* - (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y) - (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y) - (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y) -*/ - -/***************************************************************************** -pmPSFFromPSFtry (psfTry): build a PSF model from a collection of -source->modelEXT entries. The PSF ignores the first 4 (independent) model -parameters and constructs a polynomial fit to the remaining as a function of -image coordinate. - input: psfTry with fitted source->modelEXT collection, pre-allocated psf -Note: some of the array entries may be NULL (failed fits); ignore them. - *****************************************************************************/ -bool pmPSFFromPSFtry (pmPSFtry *psfTry) -{ - PS_ASSERT_PTR_NON_NULL(psfTry, false); - PS_ASSERT_PTR_NON_NULL(psfTry->sources, false); - - pmPSF *psf = psfTry->psf; - - // construct the fit vectors from the collection of objects - psVector *x = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *y = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - psVector *z = psVectorAlloc (psfTry->sources->n, PS_TYPE_F32); - - // construct the x,y terms - for (int i = 0; i < psfTry->sources->n; i++) { - pmSource *source = psfTry->sources->data[i]; - if (source->modelEXT == NULL) - continue; - - x->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_XPOS]; - y->data.F32[i] = source->modelEXT->params->data.F32[PM_PAR_YPOS]; - } - - if (!pmPSFFitShapeParams (psf, psfTry->sources, x, y, psfTry->mask)) { - psFree(x); - psFree(y); - psFree(z); - return false; - } - - // skip the unfitted parameters (X, Y, Io, Sky) and the shape parameters (SXX, SYY, SXY) - // apply the values of Nx, Ny determined above for E0,E1,E2 to the remaining parameters - for (int i = 0; i < psf->params->n; i++) { - switch (i) { - case PM_PAR_SKY: - case PM_PAR_I0: - case PM_PAR_XPOS: - case PM_PAR_YPOS: - case PM_PAR_SXX: - case PM_PAR_SYY: - case PM_PAR_SXY: - continue; - default: - break; - } - - // select the per-object fitted data for this parameter - for (int j = 0; j < psfTry->sources->n; j++) { - pmSource *source = psfTry->sources->data[j]; - if (source->modelEXT == NULL) continue; - z->data.F32[j] = source->modelEXT->params->data.F32[i]; - } - - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = psf->trendNx; - binning->nYruff = psf->trendNy; - binning->nXfine = psf->fieldNx; - binning->nYfine = psf->fieldNy; - - if (psf->psfTrendMode == PM_TREND_MAP) { - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - } - - // free existing trend, re-alloc - psFree (psf->params->data[i]); - psf->params->data[i] = pmTrend2DNoImageAlloc (psf->psfTrendMode, binning, psf->psfTrendStats); - psFree (binning); - - // fit the collection of measured parameters to the PSF 2D model - // the mask is carried from previous steps and updated with this operation - // the weight is either the flux error or NULL, depending on 'psf->poissonErrorParams' - if (!pmTrend2DFit (psf->params->data[i], psfTry->mask, 0xff, x, y, z, NULL)) { - psError(PS_ERR_UNKNOWN, false, "failed to build psf model for parameter %d", i); - psFree(x); - psFree(y); - psFree(z); - return false; - } - } - - // test dump of star parameters vs position (compare with fitted values) - if (psTraceGetLevel("psModules.objects") >= 4) { - FILE *f = fopen ("params.dat", "w"); - - for (int j = 0; j < psfTry->sources->n; j++) { - pmSource *source = psfTry->sources->data[j]; - if (source == NULL) continue; - if (source->modelEXT == NULL) continue; - - pmModel *modelPSF = pmModelFromPSF (source->modelEXT, psf); - - fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[PM_PAR_XPOS], source->modelEXT->params->data.F32[PM_PAR_YPOS]); - - for (int i = 0; i < psf->params->n; i++) { - if (psf->params->data[i] == NULL) continue; - fprintf (f, "%f %f : ", source->modelEXT->params->data.F32[i], modelPSF->params->data.F32[i]); - } - fprintf (f, "%f %d\n", source->modelEXT->chisq, source->modelEXT->nIter); - - psFree(modelPSF); - } - fclose (f); - } - - psFree (x); - psFree (y); - psFree (z); - return true; -} - -bool pmPSFFitShapeParams (pmPSF *psf, psArray *sources, psVector *x, psVector *y, psVector *srcMask) { - - // we are doing a robust fit. after each pass, we drop points which are more deviant than - // three sigma. mask is currently updated for each pass. - - // The shape parameters (SXX, SXY, SYY) are strongly coupled. We have to handle them very - // carefully. First, we convert them to the Ellipse Polarization terms (E0, E1, E2) for - // each source and fit this set of parameters. These values are less tightly coupled, but - // are still inter-related. The fitted values do a good job of constraining the major axis - // and the position angle, but the minor axis is weakly measured. When we apply the PSF - // model to construct a source model, we convert the fitted values of E0,E1,E2 to the shape - // parameters, with the constraint that the minor axis must be greater than a minimum - // threshold. - - // XXX re-read the sextractor manual on handling 'infinitely thin' sources... - - // convert the measured source shape paramters to polarization terms - psVector *e0 = psVectorAlloc (sources->n, PS_TYPE_F32); - psVector *e1 = psVectorAlloc (sources->n, PS_TYPE_F32); - psVector *e2 = psVectorAlloc (sources->n, PS_TYPE_F32); - psVector *mag = psVectorAlloc (sources->n, PS_TYPE_F32); - - for (int i = 0; i < sources->n; i++) { - // skip any masked sources (failed to fit one of the model steps or get a magnitude) - if (srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - - pmSource *source = sources->data[i]; - assert (source->modelEXT); // all unmasked sources should have modelEXT - - psEllipsePol pol = pmPSF_ModelToFit (source->modelEXT->params->data.F32); - - e0->data.F32[i] = pol.e0; - e1->data.F32[i] = pol.e1; - e2->data.F32[i] = pol.e2; - - float flux = source->modelEXT->params->data.F32[PM_PAR_I0]; - mag->data.F32[i] = (flux > 0.0) ? -2.5*log(flux) : -100.0; - } - - if (psf->psfTrendMode == PM_TREND_MAP) { - float scatterTotal = 0.0; - float scatterTotalMin = FLT_MAX; - int entryMin = -1; - - psVector *dz = NULL; - psVector *mask = psVectorAlloc (sources->n, PS_TYPE_VECTOR_MASK); - - // check the fit residuals and increase Nx,Ny until the error is minimized - // pmPSFParamTrend increases the number along the longer of x or y - for (int i = 1; i <= PS_MAX (psf->trendNx, psf->trendNy); i++) { - - // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) - for (int i = 0; i < mask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - if (!pmPSFFitShapeParamsMap (psf, i, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { - break; - } - - // store the resulting scatterTotal values and the scales, redo the best - if (scatterTotal < scatterTotalMin) { - scatterTotalMin = scatterTotal; - entryMin = i; - } - } - if (entryMin == -1) { - psError (PS_ERR_UNKNOWN, false, "failed to find image map for shape params"); - return false; - } - - // copy srcMask to mask (we do not want the mask values set in pmPSFFitShapeParamsMap to be sticky) - for (int i = 0; i < mask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - if (!pmPSFFitShapeParamsMap (psf, entryMin, &scatterTotal, mask, x, y, mag, e0, e1, e2, dz)) { - psAbort ("failed pmPSFFitShapeParamsMap on second pass?"); - } - - pmTrend2D *trend = psf->params->data[PM_PAR_E0]; - psf->trendNx = trend->map->map->numCols; - psf->trendNy = trend->map->map->numRows; - - // copy mask back to srcMask - for (int i = 0; i < mask->n; i++) { - srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - - psFree (mask); - psFree (dz); - } else { - - // XXX iterate Nx, Ny based on scatter? - // XXX we force the x & y order to be the same - // modify the order to correspond to the actual number of matched stars: - int order = PS_MAX (psf->trendNx, psf->trendNy); - if ((sources->n < 15) && (order >= 3)) order = 2; - if ((sources->n < 11) && (order >= 2)) order = 1; - if ((sources->n < 8) && (order >= 1)) order = 0; - if ((sources->n < 3)) { - psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); - return false; - } - psf->trendNx = order; - psf->trendNy = order; - - // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. - // This way, the parameters masked by one of the fits will be applied to the others - for (int i = 0; i < psf->psfTrendStats->clipIter; i++) { - - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); - psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); - - pmTrend2D *trend = NULL; - float mean, stdev; - - // XXX we are using the same stats structure on each pass: do we need to re-init it? - bool status = true; - - trend = psf->params->data[PM_PAR_E0]; - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e0, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0->n); - pmSourceVisualPSFModelResid (trend, x, y, e0, srcMask); - - trend = psf->params->data[PM_PAR_E1]; - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e1, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1->n); - pmSourceVisualPSFModelResid (trend, x, y, e1, srcMask); - - trend = psf->params->data[PM_PAR_E2]; - status &= pmTrend2DFit (trend, srcMask, 0xff, x, y, e2, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2->n); - pmSourceVisualPSFModelResid (trend, x, y, e2, srcMask); - - if (!status) { - psError (PS_ERR_UNKNOWN, true, "failed to fit polynomial to shape params"); - return false; - } - } - } - - // test dump of the psf parameters - if (psTraceGetLevel("psModules.objects") >= 4) { - FILE *f = fopen ("pol.dat", "w"); - fprintf (f, "# x y : e0obs e1obs e2obs : e0fit e1fit e2fit : mask\n"); - for (int i = 0; i < e0->n; i++) { - fprintf (f, "%f %f : %f %f %f : %f %f %f : %d\n", - x->data.F32[i], y->data.F32[i], - e0->data.F32[i], e1->data.F32[i], e2->data.F32[i], - pmTrend2DEval (psf->params->data[PM_PAR_E0], x->data.F32[i], y->data.F32[i]), - pmTrend2DEval (psf->params->data[PM_PAR_E1], x->data.F32[i], y->data.F32[i]), - pmTrend2DEval (psf->params->data[PM_PAR_E2], x->data.F32[i], y->data.F32[i]), - srcMask->data.PS_TYPE_VECTOR_MASK_DATA[i]); - } - fclose (f); - } - - psFree (e0); - psFree (e1); - psFree (e2); - psFree (mag); - return true; -} - -// fit the shape variations as a psImageMap for the given scale factor -bool pmPSFFitShapeParamsMap (pmPSF *psf, int scale, float *scatterTotal, psVector *mask, psVector *x, psVector *y, psVector *mag, psVector *e0obs, psVector *e1obs, psVector *e2obs, psVector *dz) { - - int Nx, Ny; - - // set the map scale to match the aspect ratio : for a scale of 1, we guarantee - // that we have a single cell - if (psf->fieldNx > psf->fieldNy) { - Nx = scale; - float AR = psf->fieldNy / (float) psf->fieldNx; - Ny = (int) (Nx * AR + 0.5); - Ny = PS_MAX (1, Ny); - } else { - Ny = scale; - float AR = psf->fieldNx / (float) psf->fieldNy; - Nx = (int) (Ny * AR + 0.5); - Nx = PS_MAX (1, Nx); - } - - // do we have enough sources for this fine of a grid? - if (x->n < 10*Nx*Ny) { - return false; - } - - // XXX check this against the exising type - pmTrend2DMode psfTrendMode = PM_TREND_MAP; - - psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = Nx; - binning->nYruff = Ny; - binning->nXfine = psf->fieldNx; - binning->nYfine = psf->fieldNy; - psImageBinningSetScale (binning, PS_IMAGE_BINNING_CENTER); - psImageBinningSetSkipByOffset (binning, psf->fieldXo, psf->fieldYo); - - psFree (psf->params->data[PM_PAR_E0]); - psFree (psf->params->data[PM_PAR_E1]); - psFree (psf->params->data[PM_PAR_E2]); - - int nIter = psf->psfTrendStats->clipIter; - psf->psfTrendStats->clipIter = 1; - psf->params->data[PM_PAR_E0] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); - psf->params->data[PM_PAR_E1] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); - psf->params->data[PM_PAR_E2] = pmTrend2DNoImageAlloc (psfTrendMode, binning, psf->psfTrendStats); - psFree (binning); - - // if the map is 1x1 (a single value), we measure the resulting ensemble scatter - - // if the map is more finely sampled, divide the values into two sets: measure the fit from - // one set and the scatter from the other set. - psVector *x_fit = NULL; - psVector *y_fit = NULL; - psVector *x_tst = NULL; - psVector *y_tst = NULL; - - psVector *e0obs_fit = NULL; - psVector *e1obs_fit = NULL; - psVector *e2obs_fit = NULL; - psVector *e0obs_tst = NULL; - psVector *e1obs_tst = NULL; - psVector *e2obs_tst = NULL; - - if (scale == 1) { - x_fit = psMemIncrRefCounter (x); - y_fit = psMemIncrRefCounter (y); - x_tst = psMemIncrRefCounter (x); - y_tst = psMemIncrRefCounter (y); - e0obs_fit = psMemIncrRefCounter (e0obs); - e1obs_fit = psMemIncrRefCounter (e1obs); - e2obs_fit = psMemIncrRefCounter (e2obs); - e0obs_tst = psMemIncrRefCounter (e0obs); - e1obs_tst = psMemIncrRefCounter (e1obs); - e2obs_tst = psMemIncrRefCounter (e2obs); - } else { - x_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - y_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - x_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - y_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e0obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e1obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e2obs_fit = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e0obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e1obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - e2obs_tst = psVectorAlloc (e0obs->n/2, PS_TYPE_F32); - for (int i = 0; i < e0obs_fit->n; i++) { - // e0obs->n == 8 or 9: - // i = 0, 1, 2, 3 : 2i+0 = 0, 2, 4, 6 - // i = 0, 1, 2, 3 : 2i+1 = 1, 3, 5, 7 - x_fit->data.F32[i] = x->data.F32[2*i+0]; - x_tst->data.F32[i] = x->data.F32[2*i+1]; - y_fit->data.F32[i] = y->data.F32[2*i+0]; - y_tst->data.F32[i] = y->data.F32[2*i+1]; - - e0obs_fit->data.F32[i] = e0obs->data.F32[2*i+0]; - e0obs_tst->data.F32[i] = e0obs->data.F32[2*i+1]; - e1obs_fit->data.F32[i] = e1obs->data.F32[2*i+0]; - e1obs_tst->data.F32[i] = e1obs->data.F32[2*i+1]; - e2obs_fit->data.F32[i] = e2obs->data.F32[2*i+0]; - e2obs_tst->data.F32[i] = e2obs->data.F32[2*i+1]; - } - } - - // the mask marks the values not used to calculate the ApTrend - psVector *fitMask = psVectorAlloc (x_fit->n, PS_TYPE_VECTOR_MASK); - // copy mask values to fitMask as a starting point - for (int i = 0; i < fitMask->n; i++) { - fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = mask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - - // we run 'clipIter' cycles clipping in each of x and y, with only one iteration each. - // This way, the parameters masked by one of the fits will be applied to the others - for (int i = 0; i < nIter; i++) { - // XXX we are using the same stats structure on each pass: do we need to re-init it? - psStatsOptions meanOption = psStatsMeanOption(psf->psfTrendStats->options); - psStatsOptions stdevOption = psStatsStdevOption(psf->psfTrendStats->options); - - pmTrend2D *trend = NULL; - float mean, stdev; - - // XXX we are using the same stats structure on each pass: do we need to re-init it? - bool status = true; - - trend = psf->params->data[PM_PAR_E0]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e0obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E0 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e0obs_fit->n); - // pmTrend2DPrintMap (trend); - psImageMapCleanup (trend->map); - // pmTrend2DPrintMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e0obs, mask); - - trend = psf->params->data[PM_PAR_E1]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e1obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E1 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e1obs_fit->n); - // pmTrend2DPrintMap (trend); - psImageMapCleanup (trend->map); - // pmTrend2DPrintMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e1obs, mask); - - trend = psf->params->data[PM_PAR_E2]; - status &= pmTrend2DFit (trend, fitMask, 0xff, x_fit, y_fit, e2obs_fit, NULL); - mean = psStatsGetValue (trend->stats, meanOption); - stdev = psStatsGetValue (trend->stats, stdevOption); - psTrace ("psModules.objects", 4, "clipped E2 : %f +/- %f keeping %ld of %ld\n", mean, stdev, psf->psfTrendStats->clippedNvalues, e2obs->n); - // pmTrend2DPrintMap (trend); - psImageMapCleanup (trend->map); - // pmTrend2DPrintMap (trend); - pmSourceVisualPSFModelResid (trend, x, y, e2obs, mask); - } - psf->psfTrendStats->clipIter = nIter; // restore default setting - - // construct the fitted values and the residuals - psVector *e0fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E0], fitMask, 0xff, x_tst, y_tst); - psVector *e1fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E1], fitMask, 0xff, x_tst, y_tst); - psVector *e2fit = pmTrend2DEvalVector (psf->params->data[PM_PAR_E2], fitMask, 0xff, x_tst, y_tst); - - psVector *e0res = (psVector *) psBinaryOp (NULL, (void *) e0obs_tst, "-", (void *) e0fit); - psVector *e1res = (psVector *) psBinaryOp (NULL, (void *) e1obs_tst, "-", (void *) e1fit); - psVector *e2res = (psVector *) psBinaryOp (NULL, (void *) e2obs_tst, "-", (void *) e2fit); - - // measure scatter for the unfitted points - // psTraceSetLevel ("psLib.math.vectorSampleStdev", 10); - // psTraceSetLevel ("psLib.math.vectorClippedStats", 10); - pmPSFShapeParamsScatter (scatterTotal, e0res, e1res, e2res, fitMask, 0xff, psStatsStdevOption(psf->psfTrendStats->options)); - // psTraceSetLevel ("psLib.math.vectorSampleStdev", 0); - // psTraceSetLevel ("psLib.math.vectorClippedStats", 0); - - psLogMsg ("psphot.psftry", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny); - psLogMsg ("psphot.psftry", PS_LOG_INFO, "systematic scatter: %f\n", *scatterTotal); - - psFree (x_fit); - psFree (y_fit); - psFree (x_tst); - psFree (y_tst); - - psFree (e0obs_fit); - psFree (e1obs_fit); - psFree (e2obs_fit); - psFree (e0obs_tst); - psFree (e1obs_tst); - psFree (e2obs_tst); - - psFree (e0fit); - psFree (e1fit); - psFree (e2fit); - - psFree (e0res); - psFree (e1res); - psFree (e2res); - - // XXX copy fitMask values back to mask - for (int i = 0; i < fitMask->n; i++) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = fitMask->data.PS_TYPE_VECTOR_MASK_DATA[i]; - } - psFree (fitMask); - - return true; -} - -// calculate the scatter of the parameters -bool pmPSFShapeParamsScatter(float *scatterTotal, psVector *e0res, psVector *e1res, psVector *e2res, psVector *mask, psVectorMaskType maskValue, psStatsOptions stdevOpt) -{ - - // psStats *stats = psStatsAlloc(stdevOpt); - psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_STDEV); - - // calculate the root-mean-square of E0, E1, E2 - float dEsquare = 0.0; - psStatsInit (stats); - if (!psVectorStats (stats, e0res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); - - psStatsInit (stats); - if (!psVectorStats (stats, e1res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); - - psStatsInit (stats); - if (!psVectorStats (stats, e2res, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(stats, stdevOpt)); - - *scatterTotal = sqrtf(dEsquare); - - psFree(stats); - return true; -} - -// calculate the minimum scatter of the parameters -bool pmPSFShapeParamsErrors(float *errorFloor, psVector *mag, psVector *e0res, psVector *e1res, - psVector *e2res, psVector *mask, int nGroup, psStatsOptions stdevOpt) -{ - - psStats *statsS = psStatsAlloc(stdevOpt); - - // measure the trend in bins with 10 values each; if < 10 total, use them all - int nBin = PS_MAX (mag->n / nGroup, 1); - - // use mag to group parameters in sequence - psVector *index = psVectorSortIndex (NULL, mag); - - // subset vectors for mag and dap values within the given range - psVector *dE0subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *dE1subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *dE2subset = psVectorAllocEmpty (nGroup, PS_TYPE_F32); - psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK); - - int n = 0; - float min = INFINITY; // Minimum error - for (int i = 0; i < nBin; i++) { - int j; - int nValid = 0; - for (j = 0; (j < nGroup) && (n < mag->n); j++, n++) { - int N = index->data.U32[n]; - dE0subset->data.F32[j] = e0res->data.F32[N]; - dE1subset->data.F32[j] = e1res->data.F32[N]; - dE2subset->data.F32[j] = e2res->data.F32[N]; - - mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N]; - if (!mask->data.PS_TYPE_VECTOR_MASK_DATA[N]) nValid ++; - } - if (nValid < 3) continue; - - dE0subset->n = j; - dE1subset->n = j; - dE2subset->n = j; - mkSubset->n = j; - - // calculate the root-mean-square of E0, E1, E2 - float dEsquare = 0.0; - psStatsInit (statsS); - if (!psVectorStats (statsS, dE0subset, NULL, mkSubset, 0xff)) { - } - dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); - - psStatsInit (statsS); - if (!psVectorStats (statsS, dE1subset, NULL, mkSubset, 0xff)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); - - psStatsInit (statsS); - if (!psVectorStats (statsS, dE2subset, NULL, mkSubset, 0xff)) { - psError(PS_ERR_UNKNOWN, false, "failure to measure stats"); - return false; - } - dEsquare += PS_SQR(psStatsGetValue(statsS, stdevOpt)); - - if (isfinite(dEsquare)) { - float err = sqrtf(dEsquare); - if (err < min) { - min = err; - } - } - } - *errorFloor = min; - - psFree (dE0subset); - psFree (dE1subset); - psFree (dE2subset); - psFree (mkSubset); - - psFree(index); - - psFree(statsS); - - return true; -} - -float psVectorSystematicError (psVector *residuals, psVector *errors, float clipFraction) { - - psAssert(residuals, "residuals cannot be NULL"); - psAssert(errors, "errors cannot be NULL"); - psAssert(residuals->n == errors->n, "residuals and errors must be the same length"); - - // given a vector of residuals and their formal errors, calculated the necessary systematic - // error needed to yield a reduced chisq of 1.0, after first tossing out the clipFraction - // highest chi-square contributors (allowed outliers) - - psVector *mask = psVectorAlloc(residuals->n, PS_TYPE_VECTOR_MASK); - psVector *chisq = psVectorAlloc(residuals->n, PS_TYPE_F32); - - // calculate the chisq vector: - int Ngood = 0; - for (int i = 0; i < residuals->n; i++) { - chisq->data.F32[i] = PS_MAX_F32; - if (!isfinite(residuals->data.F32[i])) continue; - if (!isfinite(errors->data.F32[i])) continue; - if (errors->data.F32[i] <= 0.0) continue; - chisq->data.F32[i] = PS_SQR(residuals->data.F32[i] / errors->data.F32[i]); - Ngood ++; - } - - psVector *index = psVectorSortIndex(NULL, chisq); - - // toss out the clipFraction highest chisq values - for (int i = 0; i < residuals->n; i++) { - int n = index->data.S32[i]; - if (i < (1.0 - clipFraction)*Ngood) { - mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 0; - } else { - mask->data.PS_TYPE_VECTOR_MASK_DATA[n] = 1; - } - } - - // Ndof ~= Ngood - // Chisq_Ndof = sum(residuals_i^2 / error_i^2) / Ndof - // choose S2 such than Chisq^sys_Ndof = sum(residuals_i^2 / (error_i^2 + S2)) / Ndof = 1.0 - - // use Newton-Raphson to solve for S2: - - // use median sigma to calculate the initial guess for S2: - psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); - psVectorStats (stats, errors, NULL, mask, 1); - float errorMedian = stats->sampleMedian; - - float nPts = 0.0; - float res2mean = 0.0; - float ChiSq = 0.0; - for (int i = 0; i < residuals->n; i++) { - int n = index->data.S32[i]; - if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue; - res2mean += PS_SQR(residuals->data.F32[n]); - ChiSq += PS_SQR(residuals->data.F32[n] / errors->data.F32[n]); - nPts += 1.0; - } - res2mean /= nPts; - ChiSq /= nPts; - - float S2guess = res2mean - PS_SQR(errorMedian); - - psLogMsg ("psModules", 3, "ChiSquare: %f, Ntotal: %ld, Ngood: %d, Nkeep: %.0f, S2 guess: %f\n", - ChiSq, residuals->n, Ngood, nPts, S2guess); - - for (int iter = 0; iter < 10; iter++) { - - ChiSq = 0.0; - float dRdS = 0.0; - for (int i = 0; i < residuals->n; i++) { - int n = index->data.S32[i]; - if (mask->data.PS_TYPE_VECTOR_MASK_DATA[n]) continue; - float error2 = PS_SQR(errors->data.F32[n]) + S2guess; - ChiSq += PS_SQR(residuals->data.F32[n]) / error2; - dRdS += PS_SQR(residuals->data.F32[n]) / PS_SQR(error2); - } - ChiSq /= nPts; - dRdS /= nPts; - - // Note the sign on dS: dRdS above is -1 * dR/dS formally - float dS = (ChiSq - 1.0) / dRdS; - S2guess += dS; - - psLogMsg ("psModules", 3, "ChiSquare: %f, dS: %f, S2 guess: %f\n", ChiSq, dS, S2guess); - } - - // free local allocations - psFree (mask); - psFree (chisq); - psFree (stats); - psFree (index); - - return (sqrt(S2guess)); -} -