Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/Makefile.am =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/Makefile.am (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/Makefile.am (revision 37403) @@ -97,6 +97,6 @@ # generalized forced photometry (including Kron, Petro, and Models) of specified positions given a specified psf psphotFullForce_SOURCES = \ - psphotFullForce.c \ - psphotFullForceArguments.c \ + psphotFullForce.c \ + psphotFullForceArguments.c \ psphotParseCamera.c \ psphotImageLoop.c \ @@ -259,5 +259,6 @@ psphotEfficiency.c \ psphotSetNFrames.c \ - psphotSourceMemory.c + psphotSourceMemory.c \ + psphotChipParams.c # not currently used Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphot.h =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphot.h (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphot.h (revision 37403) @@ -432,10 +432,8 @@ typedef struct { - float fRmajorMin; - float fRmajorMax; - float fRmajorDel; - float fRminorMin; - float fRminorMax; - float fRminorDel; + float Q; + float NSigma; + float clampSN; + int extModelType; } psphotGalaxyShapeOptions; @@ -537,6 +535,7 @@ bool psphotGalaxyShape_Threaded (psThreadJob *job); bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize); -bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults); +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults); psphotGalaxyShapeOptions *psphotGalaxyShapeOptionsAlloc(); +bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *defaultOptions); bool psphotRadialProfileWings (pmConfig *config, const pmFPAview *view, const char *filerule); @@ -572,6 +571,10 @@ bool psphotSourceMemoryReadout(pmConfig *config, const pmFPAview *view, const char *filerule, int index); +bool psphotChipParams(pmConfig *config, const pmFPAview *view, const char *filerule); +bool psphotChipParams_Threaded(psThreadJob *job); + const char * psphotGetFilerule(const char *baseRule); extern bool psphotINpsphotStack; + #endif Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChipParams.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChipParams.c (revision 37403) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChipParams.c (revision 37403) @@ -0,0 +1,293 @@ +# include "psphotInternal.h" + +bool psphotChipParamsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources); +bool psphotChipParams_Threaded (psThreadJob *job); +bool psphotChipParamsSource (pmSource *source, psArray *backmaps); + +typedef struct { + int index; + short chipNum; + psString name; + float C0x; + float C1x; + float C2x; + float C0y; + float C1y; + float C2y; + int xMin; + int xMax; + int yMin; + int yMax; +} psphotChipBackmap ; + +psArray * psphotChipParamsParseBackmaps(pmReadout *readout); + +bool psphotChipParams (pmConfig *config, const pmFPAview *view, const char *filerule) +{ + bool status = true; + + fprintf (stdout, "\n"); + psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Chip Parameters ---"); + + // select the appropriate recipe information + psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); + psAssert (recipe, "missing recipe?"); + + int num = psphotFileruleCount(config, filerule); + + // skip the chisq image (optionally?) + int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM"); + if (!status) chisqNum = -1; + + // loop over the available readouts + for (int i = 0; i < num; i++) { + if (i == chisqNum) continue; // skip chisq image + + // 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->newSources ? detections->newSources : detections->allSources; + psAssert (sources, "missing sources?"); + + if (!psphotChipParamsReadout (config, recipe, view, filerule, readout, sources)) { + psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i); + return false; + } + } + + return true; +} + +bool psphotChipParamsReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char *filerule, pmReadout *readout, psArray *sources) { + + bool status = false; + + if (!sources->n) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping chipParams"); + return true; + } + + psTimerStart ("psphot.chipparams"); + + // determine the number of allowed threads + int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads + if (!status) { + nThreads = 0; + } + psArray *backmaps = psphotChipParamsParseBackmaps(readout); + if (!backmaps) { + psError (PSPHOT_ERR_UNKNOWN, false, "failed to parse chip backmaps for %s", filerule); + return false; + } + + // threaded measurement of the source magnitudes + // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts) + int Cx = 1, Cy = 1; + psphotChooseCellSizes (&Cx, &Cy, readout, nThreads); + + psArray *cellGroups = psphotAssignSources (Cx, Cy, sources); + + for (int i = 0; i < cellGroups->n; i++) { + + psArray *cells = cellGroups->data[i]; + + for (int j = 0; j < cells->n; j++) { + + // allocate a job -- if threads are not defined, this just runs the job + psThreadJob *job = psThreadJobAlloc ("PSPHOT_CHIP_PARAMS"); + + psArrayAdd(job->args, 1, cells->data[j]); // sources + psArrayAdd(job->args, 1, backmaps); // sources + +// set this to 0 to run without threading +# if (1) + if (!psThreadJobAddPending(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(backmaps); + return false; + } +# else + if (!psphotChipParams_Threaded(job)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + // psFree(AnalysisRegion); + psFree(backmaps); + return false; + } + psFree(job); +# endif + } + + // wait for the threads to finish and manage results + if (!psThreadPoolWait (false, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to guess model."); + psFree(backmaps); + return false; + } + + // we have only supplied one type of job, so we can assume the types here + psThreadJob *job = NULL; + while ((job = psThreadJobGetDone()) != NULL) { + if (job->args->n < 1) { + fprintf (stderr, "error with job\n"); + } + psFree(job); + } + } + psFree (cellGroups); + psFree(backmaps); + + psLogMsg ("psphot.chip.params", PS_LOG_WARN, "measure chip paramters : %f sec for %ld objects\n", psTimerMark ("psphot.chipparams"), sources->n); + return true; +} + +bool psphotChipParams_Threaded (psThreadJob *job) { + + psArray *sources = job->args->data[0]; + psArray *backmaps = job->args->data[1]; + + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + if (!source->peak) continue; // XXX how can we have a peak-less source? + + psphotChipParamsSource(source, backmaps); + } + + return true; +} + + +bool psphotChipParamsSource (pmSource *source, psArray *backmaps) { + + PS_ASSERT_PTR_NON_NULL(source, false); + PS_ASSERT_PTR_NON_NULL(source->peak, false); + PS_ASSERT_PTR_NON_NULL(source->pixels, false); + + float warp_x = source->peak->xf; + float warp_y = source->peak->yf; + + int chipForSource = -1; + for (int i=0; i < backmaps->n; i++) { + psphotChipBackmap *thisChip = backmaps->data[i]; + + float chip_x = thisChip->C0x + thisChip->C1x * warp_x + thisChip->C2x * warp_y; + float chip_y = thisChip->C0y + thisChip->C1y * warp_x + thisChip->C2y * warp_y; + if (chip_x >= thisChip->xMin && chip_x <= thisChip->xMax && + chip_y >= thisChip->yMin && chip_y <= thisChip->yMax) { + + chipForSource = thisChip->chipNum; + if (0) { + printf("source %8d on chip: %s %2d at %6.1f, %6.1f\n", source->chipNum, thisChip->name, + chipForSource, chip_x, chip_y); + } + source->chipNum = chipForSource; + source->chipX = chip_x; + source->chipY = chip_y; + // Since it is on this chip we are done for this source + break; + } + } + if (chipForSource < 0) return false; + + return true; +} + + +psArray * psphotChipParamsParseBackmaps(pmReadout *readout) { + + // XXX: I'm assuming the structure of the FPA here. Find a more usual way + // to find header + psMetadata *header = readout->parent->parent->parent->hdu->header; + if (!header) { + psError(PS_ERR_UNKNOWN, true, "Unable to find header for readout."); + return NULL; + } + + psArray *backmaps = psArrayAllocEmpty(6); + + // find the keywords that define the warp to chip transformations stored by + // pswarpUpdateMetadata() + + for (int i = 0; ; i++) { + bool mdok; + char keyword[16]; + + sprintf(keyword, "SRC_%04d", i); + psString chipName = psMetadataLookupStr(&mdok, header, keyword); + if (!chipName) { + break; + } + psphotChipBackmap *thisChip = psAlloc(sizeof(psphotChipBackmap)); + psArrayAdd(backmaps, 1, thisChip); + psFree(thisChip); + + thisChip->index = i; + thisChip->name = chipName; // not incrementing the reference this gets freed before header + + // XXX: If chipName begins with XY assume gpc1 style class_id and use number as ID + if (!strncmp(chipName, "XY", 2)) { + thisChip->chipNum = atoi(chipName+2); + } else { + // otherwise just use the index + thisChip->chipNum = thisChip->index; + } + + sprintf(keyword, "SEC_%04d", i); + psString chipRegionStr = psMetadataLookupStr(&mdok, header, keyword); + if (!chipRegionStr) { + psError(PS_ERR_UNKNOWN, true, "failed to find required keyword %s", keyword); + psFree(backmaps); + return NULL; + } + int nvals = sscanf(chipRegionStr, "[%d:%d,%d:%d]", &thisChip->xMin, &thisChip->xMax, &thisChip->yMin, &thisChip->yMax); + if (nvals != 4) { + psError(PS_ERR_UNKNOWN, true, "failed to parse %s: %s", keyword, chipRegionStr); + psFree(backmaps); + return NULL; + } + // convert from fits extent convention to zero offset array + #define COORD_SLACK 2 + thisChip->xMin -= COORD_SLACK; thisChip->xMax += COORD_SLACK; + thisChip->yMin -= COORD_SLACK; thisChip->yMax += COORD_SLACK; + + sprintf(keyword, "MPX_%04d", i); + psString mpx = psMetadataLookupStr(&mdok, header, keyword); + if (!mpx) { + psError(PS_ERR_UNKNOWN, true, "failed to find %s", keyword); + psFree(backmaps); + return NULL; + } + nvals = sscanf(mpx, "[%f,%f,%f]", &thisChip->C0x, &thisChip->C1x, &thisChip->C2x); + if (nvals != 3) { + psError(PS_ERR_UNKNOWN, true, "failed to parse %s: %s", keyword, mpx); + psFree(backmaps); + return NULL; + } + sprintf(keyword, "MPY_%04d", i); + psString mpy = psMetadataLookupStr(&mdok, header, keyword); + if (!mpy) { + psError(PS_ERR_UNKNOWN, true, "failed to find %s", keyword); + psFree(backmaps); + return NULL; + } + nvals = sscanf(mpy, "[%f,%f,%f]", &thisChip->C0y, &thisChip->C1y, &thisChip->C2y); + if (nvals != 3) { + psError(PS_ERR_UNKNOWN, true, "failed to parse %s: %s", keyword, mpx); + psFree(backmaps); + return NULL; + } + } + + if (!backmaps->n) { + psFree(backmaps); + return NULL; + } + + return backmaps; +} Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChoosePSF.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChoosePSF.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotChoosePSF.c (revision 37403) @@ -448,5 +448,5 @@ // get the model full-width at half-max - float FWHM_MAJOR = 2*modelPSF->modelRadius (modelPSF->params, 0.5); + float FWHM_MAJOR = 2*modelPSF->class->modelRadius (modelPSF->params, 0.5); // XXX make sure this is consistent with the re-definition of PM_PAR_SXX Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEfficiency.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEfficiency.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEfficiency.c (revision 37403) @@ -71,5 +71,5 @@ continue; } - float flux = normModel->modelFlux(normModel->params); // Total flux for peak of 1.0 + float flux = normModel->class->modelFlux(normModel->params); // Total flux for peak of 1.0 psFree(normModel); if (!isfinite(flux)) { @@ -424,5 +424,5 @@ continue; } - float sourceRadius = PS_MAX(radius, model->modelRadius(model->params, minFlux)); // Radius for source + float sourceRadius = PS_MAX(radius, model->class->modelRadius(model->params, minFlux)); // Radius for source psFree(model); if (!isfinite(sourceRadius)) { Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEvalFLT.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEvalFLT.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotEvalFLT.c (revision 37403) @@ -53,5 +53,5 @@ } - keep = model->modelFitStatus(model); + keep = model->class->modelFitStatus(model); if (keep) return true; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceReadout.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceReadout.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceReadout.c (revision 37403) @@ -106,4 +106,6 @@ psphotKronFlux (config, view, filerule); + psphotChipParams (config, view, filerule); + // measure petro fluxes psphotPetroFlux (config, view, filerule); Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceSummary.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceSummary.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceSummary.c (revision 37403) @@ -186,13 +186,27 @@ psFitsClose(fits); - // Convert to an array indexed by ID - psArray *sortedTable = psArrayAlloc(4*inTable->n); + // Convert to a set of arrays indexed by model type + 1 + // which contain pointers to arrays indexed by ID +#define MAX_MODEL_TYPE 10 + psArray *sortedTables = psArrayAlloc(MAX_MODEL_TYPE+1); for (int i=0; in; i++) { psMetadata *row = inTable->data[i]; psS32 ID = psMetadataLookupS32(&status, row, "ID"); + psS32 modelType = psMetadataLookupS32(&status, row, "MODEL_TYPE"); + // XXX: need to use the lookup table functions to be ready for changes in the model type numbers + if (modelType+1 >= MAX_MODEL_TYPE) { + psError(PSPHOT_ERR_CONFIG, false, "found modelType %d max allowed is %d", modelType, MAX_MODEL_TYPE); + return false; + } + psArray *sortedTable = sortedTables->data[modelType+1]; + if (!sortedTable) { + sortedTable = psArrayAlloc(4*inTable->n); + sortedTables->data[modelType+1] = sortedTable; + // dont' free sortedTable the array of tables gets our reference + } if (ID >= sortedTable->n) { sortedTable = psArrayRealloc(sortedTable, 2*ID); // Why doesn't psArrayRealloc do this????? - sortedTable->n = sortedTable->nalloc; + sortedTable->n = sortedTable->nalloc; } if (sortedTable->data[ID]) { @@ -202,7 +216,7 @@ sortedTable->data[ID] = psMemIncrRefCounter(row); } - psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLES", PS_META_REPLACE, "cff tables", sortedTables); psFree(inTable); - psFree(sortedTable); + psFree(sortedTables); return true; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceSummaryReadout.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotFullForceSummaryReadout.c (revision 37403) @@ -3,23 +3,31 @@ typedef struct { + psF32 Q; + psF32 NSigma; + psF32 clampSN; +#ifdef notdef + // these are per object now int numTrials; - psF32 fRmajorMin; - psF32 fRmajorMax; - psF32 fRmajorDel; - psF32 fRminorMin; - psF32 fRminorMax; - psF32 fRminorDel; - psVector *rMajor; - psVector *rMinor; - psArray *zeroPt; - psArray *exptime; - psArray *cffTable; -} galaxyShapeOptions; - -static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options); + psF64 fRmajorMin; + psF64 fRmajorMax; + psF64 fRmajorDel; + psF64 fRminorMin; + psF64 fRminorMax; + psF64 fRminorDel; + psVector *fRmajor; + psVector *fRminor; +#endif + psArray *zeroPt; // zero points for each input used with exptime to scale flux + psArray *exptime; // exposure times for each input + psArray *cffTables; // one for each model type index is (model_type + 1) entry 0 is no extended model (star) +} galaxyShapeSummaryOptions; + + + +static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeSummaryOptions *options); static pmPhotObj *findObjectForSource(psArray **pObjects, pmSource *source); -static bool setOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors); -static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe); +static bool setOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe, bool saveVectors); +static bool checkOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe); @@ -47,13 +55,10 @@ psMetadataAddF32(outputCell->concepts, PS_LIST_TAIL, "CELL.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", outputExptime); - // don't think this one matters -// psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.EXPOSURE", PS_META_REPLACE, "Exposure time (sec)", - // outputExptime); // Create objects from the various input's sources // loop over the available readouts -// psArray *readouts = psArrayAlloc(num); + psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32); - galaxyShapeOptions options; + galaxyShapeSummaryOptions options; psArray *objects = NULL; for (int index = 0; index < num; index++) { @@ -72,6 +77,6 @@ return false; } - options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE"); - if (!options.cffTable) { + options.cffTables = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLES"); + if (!options.cffTables) { psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments."); return false; @@ -91,6 +96,4 @@ psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime; fluxScaleFactor->data.F32[index] = scaleFactor; - -// readouts->data[index] = readout; // find detections @@ -177,15 +180,27 @@ } -static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options) { +static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeSummaryOptions *options) { pmSource *outSrc = NULL; - psVector *sumWeightedFlux = NULL; - psVector *sumInvSig2 = NULL; - psVector *numerator = NULL; -// psVector *tmp = NULL; - psF32 totalNPix = 0; - long vectorLength = 0; - + + pmSource *source0 = obj->sources->data[0]; + int ID = source0->seq; + + // Find the parameters from the cff tables. nModels is the number of extended model types that + // this source has an entry for + int nModels = 0; + // skip table 0 which is modelType -1 not extended + for (int i=1; icffTables->n; i++) { + psArray *table = options->cffTables->data[i]; + if (!table) continue; + if (ID >= table->n) continue; + if (table->data[ID]) { + nModels++; + } + } + + // Loop over sources for this object. Start output source for first input that has + // good pixel frac for (int i=0; i < obj->sources->n; i++) { pmSource *source = obj->sources->data[i]; @@ -212,87 +227,193 @@ outSrc->extpars = psMemIncrRefCounter(source->extpars); } - if (source->modelEXT) { - outSrc->modelEXT = psMemIncrRefCounter(source->modelEXT); - } - if (source->modelFits) { - outSrc->modelFits = psMemIncrRefCounter(source->modelFits); - } - } - - if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) { - if (numerator == NULL) { - vectorLength = source->galaxyFits->chisq->n; - // tmp = psVectorAlloc(vectorLength, PS_TYPE_F32); - sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32); - psVectorInit(sumWeightedFlux, 0.0); - sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32); - psVectorInit(sumInvSig2, 0.0); - numerator = psVectorAlloc(vectorLength, PS_TYPE_F32); - psVectorInit(numerator, 0.0); - } - - psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d", - vectorLength, options->numTrials); - psAssert(source->galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld", - source->galaxyFits->chisq->n, vectorLength); - - - psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; + } + if (nModels && source->modelFits && source->modelFits->n) { + if (!outSrc->modelFits) { + outSrc->modelFits = psArrayAllocEmpty(nModels); + } + for (int iModel = 0; iModel < source->modelFits->n; iModel++) { + pmModel *newModel = source->modelFits->data[iModel]; + pmModel *outModel = NULL; + for (int jModel = 0; jModel < outSrc->modelFits->n; jModel++) { + outModel = outSrc->modelFits->data[jModel]; + if (newModel->type == outModel->type) { + // already have an entry for this type + break; + } + outModel = NULL; + } + if (!outModel) { + // no model of this type in output source yet. Add this one. + // Note: the parameters that we need (position, angle, index, if applicable) + // are the same for all sources so copying any of them is fine + psArrayAdd(outSrc->modelFits, 1, newModel); + } + } + } + } + + if (!outSrc) { + // no good measurements for this source (how?) + return NULL; + } + + // now loop over the model fits and summarize the galaxy shape measurements for each + for (int iModel=0; outSrc->modelFits && iModel < outSrc->modelFits->n; iModel++) { + pmModel *model = outSrc->modelFits->data[iModel]; + int modelType = model->type; + psAssert(modelType >= 0 && modelType+1 < options->cffTables->n, "model type out of range"); + + psArray *cffTable = options->cffTables->data[modelType+1]; + psAssert(cffTable, "missing cff table?"); + + psMetadata *cffRow = cffTable->data[outSrc->seq]; + if (!cffRow) { + psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d model %d", outSrc->seq, modelType); + psFree(outSrc); + return NULL; + } + bool mdok; + psF32 R_MAJOR = psMetadataLookupF32(&mdok, cffRow, "R_MAJOR"); + if (!mdok) { + psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d type %d", outSrc->seq, modelType); + psFree(outSrc); + return NULL; + } + psF32 R_MINOR = psMetadataLookupF32(&mdok, cffRow, "R_MINOR"); + if (!mdok) { + psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d type %d", outSrc->seq, modelType); + psFree(outSrc); + return NULL; + } + + psVector *sumWeightedFlux = NULL; + psVector *sumInvSig2 = NULL; + psVector *numerator = NULL; + psVector *fRmajorVec = NULL; + psVector *fRminorVec = NULL; + psF32 totalNPix = 0; + long vectorLength = 0; + psF32 fRmajorMin = 0; + psF32 fRmajorMax = 0; + psF32 fRmajorDel = 0; + psF32 fRminorMin = 0; + psF32 fRminorMax = 0; + psF32 fRminorDel = 0; + for (int i=0; i < obj->sources->n; i++) { + pmSource *source = obj->sources->data[i]; + + if (source->pixWeightNotPoor < .9) continue; + + // The only parameters that we are summarizing currently are galaxy shapes. + if (!source->modelFits) continue; + if (!source->galaxyFits) continue; + + // XXX: put this into a function + pmSourceGalaxyFits *galaxyFits = NULL; + for (int j=0; jgalaxyFits->n; j++) { + galaxyFits = source->galaxyFits->data[j]; + if (galaxyFits && galaxyFits->modelType == modelType) break; + galaxyFits = NULL; + } + + if (galaxyFits && isfinite(galaxyFits->nPix) && galaxyFits->chisq->n) { + if (numerator == NULL) { + // first source with galaxyFits allocate accumulators + vectorLength = galaxyFits->chisq->n; + sumWeightedFlux = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(sumWeightedFlux, 0.0); + sumInvSig2 = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(sumInvSig2, 0.0); + numerator = psVectorAlloc(vectorLength, PS_TYPE_F32); + psVectorInit(numerator, 0.0); + + // Initialize the fractional radius vectors + // save these to initialize the output galaxy fits + fRmajorMin = galaxyFits->fRmajorMin; + fRmajorMax = galaxyFits->fRmajorMax; + fRmajorDel = galaxyFits->fRmajorDel; + fRminorMin = galaxyFits->fRminorMin; + fRminorMax = galaxyFits->fRminorMax; + fRminorDel = galaxyFits->fRminorDel; + fRmajorVec = psVectorAlloc(vectorLength, PS_TYPE_F32); + fRminorVec = psVectorAlloc(vectorLength, PS_TYPE_F32); + int k = 0; + for (float fRmajor = fRmajorMin; fRmajor < fRmajorMax + 0.5*fRmajorDel; fRmajor += fRmajorDel) { + for (float fRminor = fRminorMin; fRminor < fRminorMax + 0.5*fRminorDel; fRminor += fRminorDel) { + fRminorVec->data.F32[k] = fRminor; + fRmajorVec->data.F32[k] = fRmajor; + k++; + } + } + psAssert(k == vectorLength, "mismatched vectors"); + } + + // Die if the lengths of the vectors is not the same in all sources +#ifdef notdef + psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d", + vectorLength, options->numTrials); +#endif + psAssert(galaxyFits->chisq->n == vectorLength, "length of chisq vectors do not match %ld %ld", + galaxyFits->chisq->n, vectorLength); + + psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; + + totalNPix += galaxyFits->nPix; + + for (int k = 0; k < vectorLength; k++) { + psF32 chisq = galaxyFits->chisq->data.F32[k]; + psF32 flux = galaxyFits->Flux->data.F32[k] * scaleFactor; + psF32 dFlux = galaxyFits->dFlux->data.F32[k] * scaleFactor; + + numerator->data.F32[k] += chisq * galaxyFits->nPix; + + psF32 invSig2 = 1.0 / (dFlux * dFlux); + sumInvSig2->data.F32[k] += invSig2; + + sumWeightedFlux->data.F32[k] += flux * invSig2; + } + } + } + + if (vectorLength) { + // allocate galaxyFits for the output source + if (!outSrc->galaxyFits) { + outSrc->galaxyFits = psArrayAllocEmpty(1); + } + pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc(); + psArrayAdd(outSrc->galaxyFits, 1, galaxyFits); + psFree(galaxyFits); + + galaxyFits->nPix = totalNPix; + galaxyFits->modelType = modelType; + psVector *fluxVec = + galaxyFits->Flux = psVectorRecycle(galaxyFits->Flux, vectorLength, PS_TYPE_F32); + psVector *dFluxVec = + galaxyFits->dFlux = psVectorRecycle(galaxyFits->dFlux, vectorLength, PS_TYPE_F32); + psVector *chisqVec = + galaxyFits->chisq = psVectorRecycle(galaxyFits->chisq, vectorLength, PS_TYPE_F32); + + galaxyFits->fRmajorMin = fRmajorMin; + galaxyFits->fRmajorMax = fRmajorMax; + galaxyFits->fRmajorDel = fRmajorDel; + galaxyFits->fRminorMin = fRminorMin; + galaxyFits->fRminorMax = fRminorMax; + galaxyFits->fRminorDel = fRminorDel; + + // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq + + int min_k = -1; + psF32 minChisq = NAN; for (int k = 0; k < vectorLength; k++) { - - psF32 chisq = source->galaxyFits->chisq->data.F32[k]; - psF32 flux = source->galaxyFits->Flux->data.F32[k] * scaleFactor; - psF32 dFlux = source->galaxyFits->dFlux->data.F32[k] * scaleFactor; - // for chisq - // Numerator = sum( nPix[i] * chisq[i] ) - // denominator will be sum( nPix ) - - -// tmp = (psVector *) psBinaryOp(tmp, - // source->galaxyFits->chisq, "*", psScalarAlloc(source->galaxyFits->nPix, PS_TYPE_F32)); -// numerator = (psVector *) psBinaryOp(numerator, numerator, "+", tmp); - - numerator->data.F32[k] += chisq * source->galaxyFits->nPix; - totalNPix += source->galaxyFits->nPix; - - // sumInvSig2 = sum( 1 / dFlux**2 ) -// tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->dFlux, "*", source->galaxyFits->dFlux); -// tmp = (psVector *) psBinaryOp(tmp, psScalarAlloc(1.0, PS_TYPE_F32), "/", tmp); -// sumInvSig2 = (psVector *) psBinaryOp(sumInvSig2, sumInvSig2, "+", tmp); -// - psF32 invSig2 = 1.0 / (dFlux * dFlux); - sumInvSig2->data.F32[k] += invSig2; - - // sumWeightedFlux = sum ( Flux / dFlux**2 ) -// tmp = (psVector *) psBinaryOp(tmp, source->galaxyFits->Flux, "*", tmp); -// sumWeightedFlux = (psVector *) psBinaryOp(sumWeightedFlux, sumWeightedFlux, "+", tmp); - sumWeightedFlux->data.F32[k] += flux * invSig2; - } - } - } - - if (outSrc) { - if (vectorLength) { - outSrc->galaxyFits = pmSourceGalaxyFitsAlloc(); - outSrc->galaxyFits->nPix = totalNPix; - outSrc->galaxyFits->Flux = psVectorRecycle(outSrc->galaxyFits->Flux, vectorLength, PS_TYPE_F32); - outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32); - outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32); - for (int k = 0; k < vectorLength; k++) { - outSrc->galaxyFits->Flux->data.F32[k] = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k]; - outSrc->galaxyFits->dFlux->data.F32[k] = 1.0 / sumInvSig2->data.F32[k]; - outSrc->galaxyFits->chisq->data.F32[k] = numerator->data.F32[k] / totalNPix; - -#ifdef nodef - outSrc->galaxyFits->Flux = (psVector *) psBinaryOp(outSrc->galaxyFits->Flux, - sumWeightedFlux, "/", sumInvSig2); - - outSrc->galaxyFits->dFlux = (psVector *) psBinaryOp(outSrc->galaxyFits->dFlux, - psScalarAlloc(1.0, PS_TYPE_F32), "/", sumInvSig2); - - outSrc->galaxyFits->chisq = (psVector *) psBinaryOp(outSrc->galaxyFits->chisq, - numerator, "/", psScalarAlloc(totalNPix, PS_TYPE_F32)); -#endif + fluxVec->data.F32[k] = sumWeightedFlux->data.F32[k] / sumInvSig2->data.F32[k]; + dFluxVec->data.F32[k] = 1.0 / sumInvSig2->data.F32[k]; + + psF32 thischisq = chisqVec->data.F32[k] = numerator->data.F32[k] / totalNPix; + + if (isfinite(thischisq) && (!isfinite(minChisq) || thischisq < minChisq)) { + min_k = k; + minChisq = thischisq; + } } @@ -300,52 +421,230 @@ psFree(sumInvSig2); psFree(sumWeightedFlux); -// psFree(tmp); - - int min_j = -1; - psF32 minChisq = NAN; - psVector *chisq = outSrc->galaxyFits->chisq; - for (int j=0; j < chisq->n; j++) { - psF32 thischisq = chisq->data.F32[j]; - if (isfinite(thischisq) && (!isfinite(minChisq) || thischisq < minChisq)) { - min_j = j; - minChisq = thischisq; + + if (min_k >= 0 && isfinite(minChisq)) { + // copy the best fit params to the model + // fractional radii with the lowest chisq + psEllipseAxes axes = pmPSF_ModelToAxes(model->params->data.F32, model->class->useReff); + + // examine the params for the trial with minimum chisq. + bool fitMajor = true; + bool fitMinor = true; + psF64 fRmajorBest = fRmajorVec->data.F32[min_k]; + if ((fabs(fRmajorBest - galaxyFits->fRmajorMin) < galaxyFits->fRmajorDel) || + (fabs(fRmajorBest - galaxyFits->fRmajorMax) < galaxyFits->fRmajorDel)) { + fitMajor = false; } - } - if (min_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) { - // copy the best fit params to the model - psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type); - if (outSrc->seq >= options->cffTable->n) { - psError(PS_ERR_PROGRAMMING, true, "object sequence number %d is larger than cff table length %ld", - outSrc->seq, options->cffTable->n); - return NULL; + psF64 fRminorBest = fRminorVec->data.F32[min_k]; + if ((fabs(fRminorBest - galaxyFits->fRminorMin) < galaxyFits->fRminorDel) || + (fabs(fRminorBest - galaxyFits->fRminorMax) < galaxyFits->fRminorDel)) { + fitMinor = false; } - psMetadata *row = options->cffTable->data[outSrc->seq]; - if (!row) { - psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq); - psFree(outSrc); - return NULL; + // If either major or minor is at one of the limits do not fit report the minimum value + bool useFit = fitMajor && fitMinor; + + // save the flux and dFlux values from entry with lowest chisq + psF64 fluxBest = fluxVec->data.F32[min_k]; + psF64 dFluxBest = dFluxVec->data.F32[min_k]; + psF64 dFlux0 = NAN; + psF64 flux0 = NAN; + + if (useFit) { + #define NUM_TRIALS_INIT 9 + psVector *major = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *minor = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *chisq = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *flux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); + + // XXX: use a recipe parameter instead of 2.2 + psF64 maxDeltaMaj = 2.2 * galaxyFits->fRmajorDel; + psF64 maxDeltaMin = 2.2 * galaxyFits->fRminorDel; + + psF64 majorMin = NAN; + psF64 majorMax = NAN; + psF64 minorMin = NAN; + psF64 minorMax = NAN; + psF64 chisqMin = NAN; + psF64 chisqMax = NAN; + for (int k = 0; k < vectorLength; k++) { + + if (fabs(fRmajorVec->data.F32[k] - fRmajorBest) < maxDeltaMaj && + fabs(fRminorVec->data.F32[k] - fRminorBest) < maxDeltaMin) { + + if (isfinite(chisqVec->data.F32[k]) && + isfinite(fluxVec->data.F32[k]) && + isfinite(dFluxVec->data.F32[k])) { + + // compute major and minor radius vectors from nominal and trial fractions + // also find the ranges in the vectors + psF64 thisMajor = R_MAJOR * fRmajorVec->data.F32[k]; + if (!isfinite(majorMin) || thisMajor < majorMin) { + majorMin = thisMajor; + } + if (!isfinite(majorMax) || thisMajor > majorMax) { + majorMax = thisMajor; + } + psVectorAppend(major, thisMajor); + + psF64 thisMinor = R_MINOR * fRminorVec->data.F32[k]; + if (!isfinite(minorMin) || thisMinor < minorMin) { + minorMin = thisMinor; + } + if (!isfinite(minorMax) || thisMinor > minorMax) { + minorMax = thisMinor; + } + psVectorAppend(minor, thisMinor); + + psF64 thisChisq = chisqVec->data.F32[k]; + if (!isfinite(chisqMin) || thisChisq < chisqMin) { + chisqMin = thisChisq; + } + if (!isfinite(chisqMax) || thisChisq > chisqMax) { + chisqMax = thisChisq; + } + psVectorAppend(chisq, thisChisq); + + psVectorAppend(flux, fluxVec->data.F32[k]); + psVectorAppend(dFlux, dFluxVec->data.F32[k]); + } + } + } + + + // see if we ever get too few good values (haven't seen this happen) + if (major->n < NUM_TRIALS_INIT) { + fprintf(stderr, "only found %ld good points wanted %d for seq: %d\n", + major->n, NUM_TRIALS_INIT, outSrc->seq); + } + + // rescale data - this helps avoid precision errors. + psF64 majorRange = majorMax - majorMin; + if (majorRange == 0) { + majorRange = 1; + } + psF64 minorRange = minorMax - minorMin; + if (minorRange == 0) { + minorRange = 1; + } + psF64 chisqRange = chisqMax - chisqMin; + if (chisqRange == 0) { + chisqRange = 1; + } + for (int k = 0; k < major->n; k++) { + major->data.F64[k] = (major->data.F64[k] - majorMin) / majorRange; + minor->data.F64[k] = (minor->data.F64[k] - minorMin) / minorRange; + chisq->data.F64[k] = (chisq->data.F64[k] - chisqMin) / chisqRange; + } + + // Fit chisq versus rMajor and rMinor + psPolynomial2D *chisqFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); + psPolynomial2D *fluxFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); + psPolynomial2D *fluxErrorFit = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); + chisqFit->coeffMask[2][2] = PS_POLY_MASK_SET; + chisqFit->coeffMask[2][1] = PS_POLY_MASK_SET; + chisqFit->coeffMask[1][2] = PS_POLY_MASK_SET; + fluxFit->coeffMask[2][2] = PS_POLY_MASK_SET; + fluxFit->coeffMask[2][1] = PS_POLY_MASK_SET; + fluxFit->coeffMask[1][2] = PS_POLY_MASK_SET; + fluxErrorFit->coeffMask[2][2] = PS_POLY_MASK_SET; + fluxErrorFit->coeffMask[2][1] = PS_POLY_MASK_SET; + fluxErrorFit->coeffMask[1][2] = PS_POLY_MASK_SET; + + bool goodFit = psVectorFitPolynomial2D (chisqFit, NULL, 0xff, chisq, NULL, major, minor); + + if (goodFit) { + // find the axes for the minimum chisq from the formula + psF64 **coeff = chisqFit->coeff; + psF64 D = 4*coeff[2][0]*coeff[0][2] - coeff[1][1]*coeff[1][1]; + psF64 major0 = ((coeff[1][1]*coeff[0][1] - 2*coeff[0][2]*coeff[1][0]) / D) * majorRange + majorMin; + psF64 minor0 = ((coeff[1][1]*coeff[1][0] - 2*coeff[2][0]*coeff[0][1]) / D) * minorRange + minorMin; + + // estimated chisq at the minimum + psF64 chisq0 = psPolynomial2DEval(chisqFit, major0, minor0) * chisqRange + chisqMin; + + // now fit the flux ... + bool goodFluxFit = psVectorFitPolynomial2D(fluxFit, NULL, 0xFF, flux, NULL, major, minor); + + // .. and compute flux at the minimum chisq position + if (goodFluxFit) { + flux0 = psPolynomial2DEval(fluxFit, major0, minor0); + } else { + flux0 = fluxBest; + } + + // .. and compute dFlux at the minimum chisq position + bool goodFluxErrorFit = psVectorFitPolynomial2D(fluxErrorFit, NULL, 0xFF, dFlux, NULL, major, minor); + if (goodFluxErrorFit) { + dFlux0 = psPolynomial2DEval(fluxErrorFit, major0, minor0); + } else { + dFlux0 = dFluxBest; + } + +#ifdef PRINTVALS + fprintf (stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %4.2f | %6.3f %6.3f %6.3f %4.2f | %7.4f %7.4f %7.1f %7.1f\n", + outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest*R_MAJOR, major0, major0/R_MAJOR, R_MINOR, fRminorBest*R_MINOR, minor0, minor0/R_MINOR, minChisq, chisq0, fluxBest, flux0); +#endif + + + axes.major = major0; + axes.minor = minor0; + model->chisq = chisq0; +#ifdef DUMPVECTORS + char fn[80]; + + sprintf(fn, "vectors/%s.%05d.txt", goodFit ? "g" : "b", outSrc->seq); + + FILE *f = fopen(fn, "w"); + fprintf(f, "#major minor chisq i\n"); + for (int k = 0; k < chisq->n; k++) { + fprintf(f, "%10.6f %10.6f %10.6f %5d\n", + major->data.F64[k]/R_MAJOR, minor->data.F64[k]/R_MINOR, chisq->data.F64[k], k); + } + fclose(f); +#endif + } else { + model->flags |= PM_MODEL_STATUS_NONCONVERGE; + +#ifdef PRINTVALS + fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f bad fit\n", + outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest); +#endif + psErrorClear(); + // psFree(outSrc->galaxyFits); + + } + psFree(chisqFit); + psFree(fluxFit); + psFree(fluxErrorFit); + psFree(major); + psFree(minor); + psFree(chisq); + psFree(flux); + psFree(dFlux); + } else { + // No fit + // Set the values based on the trial with the best chisq + axes.major = R_MAJOR * fRmajorBest; + axes.minor = R_MINOR * fRminorBest; + model->chisq = minChisq; + flux0 = fluxBest; + dFlux0 = dFluxBest; +#ifdef PRINTVALS + fprintf(stderr, "%4d %3d %3ld | %6.3f %6.3f %6.3f %6.3f skip fit\n", + outSrc->seq, min_k, 0L, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest); +#endif } - bool mdok; - psF32 rMajor = psMetadataLookupF32(&mdok, row, "R_MAJOR"); - if (!mdok) { - psError (PS_ERR_PROGRAMMING, true, "can't find R_MAJOR for object %d", outSrc->seq); - psFree(outSrc); - return NULL; - } - psF32 rMinor = psMetadataLookupF32(&mdok, row, "R_MINOR"); - if (!mdok) { - psError (PS_ERR_PROGRAMMING, true, "can't find R_MINOR for object %d", outSrc->seq); - psFree(outSrc); - return NULL; - } - axes.major = rMajor * options->rMajor->data.F32[min_j]; - axes.minor = rMinor * options->rMinor->data.F32[min_j]; - pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type); - outSrc->modelEXT->chisq = minChisq; - outSrc->modelEXT->mag = -2.5 * log10(outSrc->galaxyFits->Flux->data.F32[min_j]); - outSrc->modelEXT->magErr = 1.0 / - (outSrc->galaxyFits->Flux->data.F32[min_j]/outSrc->galaxyFits->dFlux->data.F32[min_j]); // 1 / SN - } - } + // now save the model parameters in the model structure + + pmPSF_AxesToModel (model->params->data.F32, axes, model->class->useReff); + + model->mag = -2.5 * log10(flux0); + model->magErr = dFlux0 / flux0; // 1 / SN + // XXX: should there be a different flag if we didn't do a fit of the chisq + model->flags |= PM_MODEL_STATUS_FITTED; + } + } + psFree(fRminorVec); + psFree(fRmajorVec); } @@ -360,11 +659,20 @@ } -static bool setOptions(galaxyShapeOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) { +static bool setOptions(galaxyShapeSummaryOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) { bool status; - bool useAnalysis; - + bool useAnalysis; // fall back to recipe if we dont' find values in analysis. Probably should no longer do this + +#ifdef notdef psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_FR_MAJOR_MIN"); +#endif + psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_Q"); psMetadata *md = useAnalysis ? readout->analysis : recipe; + GETVAL(Q, "GALAXY_SHAPES_Q"); + GETVAL(NSigma, "GALAXY_SHAPES_NSIGMA"); + GETVAL(clampSN, "GALAXY_SHAPES_CLAMP_SN"); + +#ifdef notdef + // these are per object now GETVAL(fRmajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN"); GETVAL(fRmajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX"); @@ -378,6 +686,6 @@ if (makeVectors) { - opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); - opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->fRminor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->fRmajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); int i = 0; for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; @@ -385,6 +693,6 @@ for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) { - opt->rMinor->data.F32[i] = fRminor; - opt->rMajor->data.F32[i] = fRmajor; + opt->fRminor->data.F32[i] = fRminor; + opt->fRmajor->data.F32[i] = fRmajor; i++; } @@ -392,7 +700,8 @@ psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials); } else { - opt->rMinor = NULL; - opt->rMajor = NULL; - } + opt->fRminor = NULL; + opt->fRmajor = NULL; + } +#endif return true; @@ -405,6 +714,6 @@ } -static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe) { - galaxyShapeOptions thisReadoutsOptions; +static bool checkOptions(galaxyShapeSummaryOptions *options, pmReadout *readout, psMetadata *recipe) { + galaxyShapeSummaryOptions thisReadoutsOptions; if (!setOptions(&thisReadoutsOptions, readout, recipe, false)) { @@ -412,4 +721,10 @@ return false; } + CHECKVAL(options, thisReadoutsOptions, Q, "mismatched Q"); + CHECKVAL(options, thisReadoutsOptions, NSigma, "mismatched NSIGMA"); + CHECKVAL(options, thisReadoutsOptions, clampSN, "mismatched cleampSN"); + +#ifdef notdef + // these are per object now CHECKVAL(options, thisReadoutsOptions, numTrials, "mismatched number of trials") CHECKVAL(options, thisReadoutsOptions, fRmajorMin, "mismatched fRmajorMin") @@ -419,4 +734,5 @@ CHECKVAL(options, thisReadoutsOptions, fRminorMax, "mismatched fRminorMax") CHECKVAL(options, thisReadoutsOptions, fRminorDel, "mismatched fRminorDel") +#endif return true; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotGalaxyShape.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotGalaxyShape.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotGalaxyShape.c (revision 37403) @@ -80,9 +80,21 @@ assert (markVal); - // maskVal is used to test for rejected pixels, and must include markVal maskVal |= markVal; - // bit-mask to mark pixels not used in analysis psphotGalaxyShapeOptions *opt = psphotGalaxyShapeOptionsAlloc(); + opt->Q = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_Q"); + psAssert (status, "missing GALAXY_SHAPES_Q"); + opt->NSigma = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_NSIGMA"); + psAssert (status, "missing GALAXY_SHAPES_NSIGMA"); + opt->clampSN = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_CLAMP_SN"); + psAssert (status, "missing GALAXY_SHAPES_NSIGMA"); + psString modelTypeToSave = psMetadataLookupStr(&status, recipe, "EXT_MODEL_TYPE_FORCE"); + if (modelTypeToSave && strcmp(modelTypeToSave, "ALL") && strcmp(modelTypeToSave, "BEST")) { + opt->extModelType = pmModelClassGetType(modelTypeToSave); + } else { + opt->extModelType = pmModelClassGetType("PS_MODEL_SERSIC"); + } + +#ifdef notdef opt->fRmajorMin = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MIN"); if (!status) opt->fRmajorMin = 0.5; opt->fRmajorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MAJOR_MAX"); if (!status) opt->fRmajorMax = 2.0; @@ -91,4 +103,5 @@ opt->fRminorMax = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_MAX"); if (!status) opt->fRminorMax = 2.0; opt->fRminorDel = psMetadataLookupF32(&status, recipe, "GALAXY_SHAPES_FR_MINOR_DEL"); if (!status) opt->fRminorDel = 0.1; +#endif // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box) @@ -238,89 +251,97 @@ bool psphotGalaxyShapeGrid (pmSource *source, pmSourceFitOptions *fitOptions, psphotGalaxyShapeOptions *opt, psImageMaskType maskVal, int psfSize) { - // we use the 0th model to define the initial guess shape - pmModel *model = source->modelFits->data[0]; - if (!model) return false; - - pmModelType modelType = model->type; - - // we are using fitOptions->mode : be sure this makes sense - pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); - if (!pcm) return false; - - // we are fitting only PM_PAR_I0; the shape elements are generated from a grid - psF32 *PAR = pcm->modelConv->params->data.F32; - - // double check that the guess is carried along... - - // I have some source guess (e0, e1, e2) - psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, modelType); - - if (!source->galaxyFits) { - source->galaxyFits = pmSourceGalaxyFitsAlloc(); - } - - float fRmajorBest = NAN; - float fRminorBest = NAN; - float chisqBest = NAN; - for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; fRmajor += opt->fRmajorDel) { - for (float fRminor = opt->fRminorMin; fRminor < opt->fRminorMax + 0.5*opt->fRminorDel; fRminor += opt->fRminorDel) { - - psEllipseAxes testAxes = guessAxes; - testAxes.major = guessAxes.major * fRmajor; - testAxes.minor = guessAxes.minor * fRminor; - - pmPSF_AxesToModel (PAR, testAxes, modelType); - - psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, true); - - int i = source->galaxyFits->chisq->n - 1; - float flux = source->galaxyFits->Flux->data.F32[i]; - if (isfinite(flux)) { - float thisChisq = source->galaxyFits->chisq->data.F32[i]; - if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) { - chisqBest = thisChisq; - fRmajorBest = fRmajor; - fRminorBest = fRminor; + + for (int iModel = 0 ; iModel < source->modelFits->n; iModel++) { + pmModel *model = source->modelFits->data[iModel]; + if (!model) return false; + + pmModelType modelType = model->type; + model->flags = PM_MODEL_STATUS_NONE; + + // we are using fitOptions->mode : be sure this makes sense + pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize); + if (!pcm) return false; + + // we are fitting only PM_PAR_I0; the shape elements are generated from a grid + psF32 *PAR = pcm->modelConv->params->data.F32; + + pmSourceGalaxyFits *galaxyFits = pmSourceGalaxyFitsAlloc(); + + // set the options for this source + if (!psphotGalaxyShapeOptionsSet(source, galaxyFits, opt)) { + // this source won't work + psFree(pcm); + psFree(galaxyFits); + continue; + } + + if (!source->galaxyFits) { + source->galaxyFits = psArrayAllocEmpty(1); + } + psArrayAdd(source->galaxyFits, 1, galaxyFits); + psFree(galaxyFits); + + // I have some source guess (e0, e1, e2) + psEllipseAxes guessAxes = pmPSF_ModelToAxes (PAR, model->class->useReff); + + float fRmajorBest = NAN; + float fRminorBest = NAN; + float chisqBest = NAN; + for (float fRmajor = galaxyFits->fRmajorMin; fRmajor < galaxyFits->fRmajorMax + 0.5*galaxyFits->fRmajorDel; fRmajor += galaxyFits->fRmajorDel) { + for (float fRminor = galaxyFits->fRminorMin; fRminor < galaxyFits->fRminorMax + 0.5*galaxyFits->fRminorDel; fRminor += galaxyFits->fRminorDel) { + + psEllipseAxes testAxes = guessAxes; + testAxes.major = guessAxes.major * fRmajor; + testAxes.minor = guessAxes.minor * fRminor; + + pmPSF_AxesToModel (PAR, testAxes, model->class->useReff); + + psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, true); + + int i = galaxyFits->chisq->n - 1; + float flux = galaxyFits->Flux->data.F32[i]; + if (isfinite(flux)) { + float thisChisq = galaxyFits->chisq->data.F32[i]; + if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) { + chisqBest = thisChisq; + fRmajorBest = fRmajor; + fRminorBest = fRminor; + } } + // reset I0 to avoid potential problems on the next iteration + PAR[PM_PAR_I0] = 1.0; } - // reset I0 to avoid potential problems on the next iteration - PAR[PM_PAR_I0] = 1.0; - } - } - -#define SAVE_BEST_MODEL -#ifdef SAVE_BEST_MODEL - // Save model with smallest chisq - if (isfinite(chisqBest)) { -#else - // Save model with nominal parameters - { - fRmajorBest = 1; - fRminorBest = 1; + } + + if (isfinite(chisqBest)) { + // now save the best fitting model as the source's extended model ... + psEllipseAxes testAxes = guessAxes; + + // ... unless this macro is defined +#ifndef SAVE_NOMINAL_MODEL + testAxes.major = guessAxes.major * fRmajorBest; + testAxes.minor = guessAxes.minor * fRminorBest; #endif - // now save the best fitting model as the source's extended model - psEllipseAxes testAxes = guessAxes; - testAxes.major = guessAxes.major * fRmajorBest; - testAxes.minor = guessAxes.minor * fRminorBest; - - pmPSF_AxesToModel (PAR, testAxes, modelType); - - psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false); - - // Replace modelEXT with this model, if the model is good - if (isfinite(PAR[PM_PAR_I0])) { - psFree (source->modelEXT); - - source->modelEXT = psMemIncrRefCounter (pcm->modelConv); - source->type = PM_SOURCE_TYPE_EXTENDED; - source->mode |= PM_SOURCE_MODE_EXTMODEL; - source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; - - // cache the model flux - pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma); - } - } - - psFree (pcm); + + pmPSF_AxesToModel (PAR, testAxes, model->class->useReff); + + psphotGalaxyShapeSource (pcm, source, galaxyFits, maskVal, psfSize, false); + + // Replace modelEXT with the best model from the first of the model fits if one of them is good + if (isfinite(PAR[PM_PAR_I0]) && modelType == opt->extModelType) { + psFree (source->modelEXT); + + source->modelEXT = psMemIncrRefCounter (pcm->modelConv); + source->type = PM_SOURCE_TYPE_EXTENDED; + source->mode |= PM_SOURCE_MODE_EXTMODEL; + source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT; + + // cache the model flux + pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma); + } + } + + psFree (pcm); + } return true; } @@ -328,5 +349,5 @@ // fit the given model to the source and find chisq & normalization // XXX is this a single-component model? sersic with a supplied index, Reff, axis ratio, and theta? -bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, psImageMaskType maskVal, int psfSize, bool saveResults) { +bool psphotGalaxyShapeSource (pmPCMdata *pcm, pmSource *source, pmSourceGalaxyFits *galaxyFits, psImageMaskType maskVal, int psfSize, bool saveResults) { PS_ASSERT_PTR_NON_NULL(source, false); @@ -379,13 +400,14 @@ pcm->modelConv->params->data.F32[PM_PAR_I0] = Io; - float flux = pcm->modelConv->modelFlux (pcm->modelConv->params); + float flux = pcm->modelConv->class->modelFlux (pcm->modelConv->params); float dflux = flux * (dIo / Io); if (saveResults) { - psVectorAppend (source->galaxyFits->Flux, flux); - psVectorAppend (source->galaxyFits->dFlux, dflux); - psVectorAppend (source->galaxyFits->chisq, Chisq); - source->galaxyFits->nPix = nPix; + psVectorAppend (galaxyFits->Flux, flux); + psVectorAppend (galaxyFits->dFlux, dflux); + psVectorAppend (galaxyFits->chisq, Chisq); + galaxyFits->nPix = nPix; + galaxyFits->modelType = pcm->modelConv->type; } @@ -404,12 +426,28 @@ psMemSetDeallocator(opt, (psFreeFunc) psphotGalaxyShapeOptionsFree); - opt->fRmajorMin = 0.5; - opt->fRmajorMax = 2.0; - opt->fRmajorDel = 0.1; - opt->fRminorMin = 0.5; - opt->fRminorMax = 2.0; - opt->fRminorDel = 0.1; - return opt; } +bool psphotGalaxyShapeOptionsSet(pmSource *source, pmSourceGalaxyFits *galaxyFits, psphotGalaxyShapeOptions *options) +{ + // XXX: put these in recipe + // doesn't make sense to use too small of a spacing + float clampedSN = source->extSN < options->clampSN ? source->extSN : options->clampSN; + + float f_del = options->Q / clampedSN; + float f_min = 1 - options->NSigma * f_del; + float f_max = 1 + options->NSigma * f_del; + + // if f_min goes negative skip this object + if (f_min < 0) return false; + + galaxyFits->fRmajorMin = f_min; + galaxyFits->fRmajorMax = f_max; + galaxyFits->fRmajorDel = f_del; + galaxyFits->fRminorMin = f_min; + galaxyFits->fRminorMax = f_max; + galaxyFits->fRminorDel = f_del; + + return true; +} + Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotKronIterate.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotKronIterate.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotKronIterate.c (revision 37403) @@ -532,4 +532,6 @@ source->moments->KronFlux = NAN; source->moments->KronFluxErr = NAN; + source->moments->KronFinner = NAN; + source->moments->KronFouter = NAN; return false; } Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLensing.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLensing.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLensing.c (revision 37403) @@ -12,4 +12,10 @@ psMetadata *recipe = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE); psAssert (recipe, "missing recipe?"); + + // perform full extended source non-linear fits? + if (!psMetadataLookupBool (&status, recipe, "LENSING_PARAMETERS")) { + psLogMsg ("psphot", PS_LOG_INFO, "skipping extended source fits\n"); + return true; + } int num = psphotFileruleCount(config, filerule); @@ -87,4 +93,5 @@ } + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "LENS_OBJ", PS_META_REPLACE, "per-object lensing stats measured", true); psLogMsg ("psphot.lensing", PS_LOG_DETAIL, "calculate lensing parameters for %d objects\n", (int) sources->n); @@ -317,4 +324,5 @@ psFree (stats); + psMetadataAddBool(readout->analysis, PS_LIST_TAIL, "LENS_PSF", PS_META_REPLACE, "psf-trend lensing stats measured", true); psLogMsg ("psphot.lensing", PS_LOG_DETAIL, "calculate lensing parameters for %d objects: %f sec\n", (int) sources->n, psTimerMark ("psphot.lensing")); Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLoadSRCTEXT.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLoadSRCTEXT.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotLoadSRCTEXT.c (revision 37403) @@ -73,5 +73,5 @@ dPAR[PM_PAR_I0] = 0.0; - pmPSF_AxesToModel (PAR, axes, modelType); + pmPSF_AxesToModel (PAR, axes, model->class->useReff); float peakFlux = 1.0; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotMakeFluxScale.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotMakeFluxScale.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotMakeFluxScale.c (revision 37403) @@ -38,5 +38,5 @@ } else { // measure the fitMag for this model - fitSum = model->modelFlux (model->params); + fitSum = model->class->modelFlux (model->params); } if (fitSum < 1.e-6) continue; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRadiusChecks.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRadiusChecks.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotRadiusChecks.c (revision 37403) @@ -60,7 +60,7 @@ if (radiusFit <= 0) { // use fixed radius if (moments == NULL) { - radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + radiusFit = model->class->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); } else { - radiusFit = model->modelRadius(model->params, 1.0); + radiusFit = model->class->modelRadius(model->params, 1.0); } model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); @@ -95,7 +95,7 @@ if (radiusFit <= 0) { // use fixed radius if (moments == NULL) { - radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); + radiusFit = model->class->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky); } else { - radiusFit = model->modelRadius(model->params, 1.0); + radiusFit = model->class->modelRadius(model->params, 1.0); } model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING); @@ -294,5 +294,5 @@ float flux = deep ? EXT_FIT_NSIGMA*EXT_FIT_SKY_SIG : 0.1 * model->params->data.F32[PM_PAR_I0]; - float rawRadius = model->modelRadius (model->params, flux); + float rawRadius = model->class->modelRadius (model->params, flux); if (isnan(rawRadius)) return false; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSersicModelClass.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSersicModelClass.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSersicModelClass.c (revision 37403) @@ -586,5 +586,5 @@ coord->data.F32[1] = (psF32) (i + source->pixels->row0); - pcm->modelConvFlux->data.F32[i][j] = pcm->modelConv->modelFunc (NULL, params, coord); + pcm->modelConvFlux->data.F32[i][j] = pcm->modelConv->class->modelFunc (NULL, params, coord); } } Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSetThreads.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSetThreads.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSetThreads.c (revision 37403) @@ -87,4 +87,9 @@ psFree(task); + task = psThreadTaskAlloc("PSPHOT_CHIP_PARAMS", 2); + task->function = &psphotChipParams_Threaded; + psThreadTaskAdd(task); + psFree(task); + return true; } Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSourceFits.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSourceFits.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotSourceFits.c (revision 37403) @@ -700,5 +700,5 @@ model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i]; - if (!model->modelGuess(model, source, maskVal, markVal)) { + if (!model->class->modelGuess(model, source, maskVal, markVal)) { model->flags |= PM_MODEL_STATUS_BADARGS; return false; @@ -723,5 +723,5 @@ model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin]; - model->modelGuess(model, source, maskVal, markVal); + model->class->modelGuess(model, source, maskVal, markVal); return true; @@ -752,5 +752,5 @@ model->params->data.F32[PM_PAR_7] = indexGuess[i]; - if (!model->modelGuess(model, source, maskVal, markVal)) { + if (!model->class->modelGuess(model, source, maskVal, markVal)) { model->flags |= PM_MODEL_STATUS_BADARGS; return false; Index: branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotVisual.c =================================================================== --- branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotVisual.c (revision 36680) +++ branches/eam_branches/ps2-tc3-20130727/psphot/src/psphotVisual.c (revision 37403) @@ -1444,14 +1444,14 @@ coord->data.F32[1] = r; coord->data.F32[0] = 0.0; - fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + fmaj->data.F32[i] = log10(source->modelPSF->class->modelFunc (NULL, params, coord)); coord->data.F32[0] = r; coord->data.F32[1] = 0.0; - fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord)); + fmin->data.F32[i] = log10(source->modelPSF->class->modelFunc (NULL, params, coord)); } psFree (coord); psFree (params); - float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); + float FWHM_MAJOR = 2.0*source->modelPSF->class->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]); float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major); if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR);