Index: trunk/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36731) +++ trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36747) @@ -3,25 +3,31 @@ typedef struct { + psF32 Q; + psF32 NSigma; + psF32 clampSN; +#ifdef notdef + // these are per object now int numTrials; - psF64 fRMajorMin; - psF64 fRMajorMax; - psF64 fRMajorDel; - psF64 fRMinorMin; - psF64 fRMinorMax; - psF64 fRMinorDel; - psVector *fRMajor; - psVector *fRMinor; + 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) -} galaxyShapeOptions; - - - -static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options); +} 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); @@ -54,5 +60,5 @@ psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32); - galaxyShapeOptions options; + galaxyShapeSummaryOptions options; psArray *objects = NULL; for (int index = 0; index < num; index++) { @@ -174,5 +180,5 @@ } -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; @@ -220,5 +226,4 @@ outSrc->extpars = psMemIncrRefCounter(source->extpars); } - outSrc->modelEXT = psFree(outSrc->modelEXT); } if (nModels && !outSrc->modelFits && source->modelFits && source->modelFits->n == nModels) { @@ -264,6 +269,14 @@ 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]; @@ -282,5 +295,4 @@ galaxyFits = NULL; } - // end of proposed function if (galaxyFits && isfinite(galaxyFits->nPix) && galaxyFits->chisq->n) { @@ -294,9 +306,31 @@ 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); @@ -339,4 +373,10 @@ 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 @@ -369,12 +409,12 @@ bool fitMajor = true; bool fitMinor = true; - psF64 fRMajorBest = options->fRMajor->data.F32[min_k]; - if ((fabs(fRMajorBest - options->fRMajorMin) < options->fRMajorDel) || - (fabs(fRMajorBest - options->fRMajorMax) < options->fRMajorDel)) { + psF64 fRmajorBest = fRmajorVec->data.F32[min_k]; + if ((fabs(fRmajorBest - galaxyFits->fRmajorMin) < galaxyFits->fRmajorDel) || + (fabs(fRmajorBest - galaxyFits->fRmajorMax) < galaxyFits->fRmajorDel)) { fitMajor = false; } - psF64 fRMinorBest = options->fRMinor->data.F32[min_k]; - if ((fabs(fRMinorBest - options->fRMinorMin) < options->fRMinorDel) || - (fabs(fRMinorBest - options->fRMinorMax) < options->fRMinorDel)) { + psF64 fRminorBest = fRminorVec->data.F32[min_k]; + if ((fabs(fRminorBest - galaxyFits->fRminorMin) < galaxyFits->fRminorDel) || + (fabs(fRminorBest - galaxyFits->fRminorMax) < galaxyFits->fRminorDel)) { fitMinor = false; } @@ -396,7 +436,7 @@ psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_INIT, PS_TYPE_F64); - // XXX: use a recipe parameter - psF64 maxDeltaMaj = 2.2 * options->fRMajorDel; - psF64 maxDeltaMin = 2.2 * options->fRMinorDel; + // XXX: use a recipe parameter instead of 2.2 + psF64 maxDeltaMaj = 2.2 * galaxyFits->fRmajorDel; + psF64 maxDeltaMin = 2.2 * galaxyFits->fRminorDel; psF64 majorMin = NAN; @@ -408,6 +448,6 @@ for (int k = 0; k < vectorLength; k++) { - if (fabs(options->fRMajor->data.F32[k] - fRMajorBest) < maxDeltaMaj && - fabs(options->fRMinor->data.F32[k] - fRMinorBest) < maxDeltaMin) { + if (fabs(fRmajorVec->data.F32[k] - fRmajorBest) < maxDeltaMaj && + fabs(fRminorVec->data.F32[k] - fRminorBest) < maxDeltaMin) { if (isfinite(chisqVec->data.F32[k]) && @@ -417,5 +457,5 @@ // compute major and minor radius vectors from nominal and trial fractions // also find the ranges in the vectors - psF64 thisMajor = R_MAJOR * options->fRMajor->data.F32[k]; + psF64 thisMajor = R_MAJOR * fRmajorVec->data.F32[k]; if (!isfinite(majorMin) || thisMajor < majorMin) { majorMin = thisMajor; @@ -426,5 +466,5 @@ psVectorAppend(major, thisMajor); - psF64 thisMinor = R_MINOR * options->fRMinor->data.F32[k]; + psF64 thisMinor = R_MINOR * fRminorVec->data.F32[k]; if (!isfinite(minorMin) || thisMinor < minorMin) { minorMin = thisMinor; @@ -449,4 +489,5 @@ } } + // see if we ever get too few good values (haven't seen this happen) @@ -521,5 +562,5 @@ #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); + 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 @@ -546,5 +587,5 @@ #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); + outSrc->seq, min_k, major->n, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest); #endif psErrorClear(); @@ -563,6 +604,6 @@ // No fit // Set the values based on the trial with the best chisq - axes.major = R_MAJOR * fRMajorBest; - axes.minor = R_MINOR * fRMinorBest; + axes.major = R_MAJOR * fRmajorBest; + axes.minor = R_MINOR * fRminorBest; model->chisq = minChisq; flux0 = fluxBest; @@ -570,5 +611,5 @@ #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); + outSrc->seq, min_k, 0L, R_MAJOR, fRmajorBest, R_MINOR, fRminorBest); #endif } @@ -583,4 +624,6 @@ } } + psFree(fRminorVec); + psFree(fRmajorVec); } @@ -595,31 +638,40 @@ } -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; // 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(fRMajorMin, "GALAXY_SHAPES_FR_MAJOR_MIN"); - GETVAL(fRMajorMax, "GALAXY_SHAPES_FR_MAJOR_MAX"); - GETVAL(fRMajorDel, "GALAXY_SHAPES_FR_MAJOR_DEL"); - GETVAL(fRMinorMin, "GALAXY_SHAPES_FR_MINOR_MIN"); - GETVAL(fRMinorMax, "GALAXY_SHAPES_FR_MINOR_MAX"); - GETVAL(fRMinorDel, "GALAXY_SHAPES_FR_MINOR_DEL"); - - opt->numTrials = ceil((opt->fRMajorMax - opt->fRMajorMin + 0.5*opt->fRMajorDel) / opt->fRMajorDel) * - ceil((opt->fRMinorMax - opt->fRMinorMin + 0.5*opt->fRMinorDel) / opt->fRMinorDel) ; + 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"); + GETVAL(fRmajorDel, "GALAXY_SHAPES_FR_MAJOR_DEL"); + GETVAL(fRminorMin, "GALAXY_SHAPES_FR_MINOR_MIN"); + GETVAL(fRminorMax, "GALAXY_SHAPES_FR_MINOR_MAX"); + GETVAL(fRminorDel, "GALAXY_SHAPES_FR_MINOR_DEL"); + + opt->numTrials = ceil((opt->fRmajorMax - opt->fRmajorMin + 0.5*opt->fRmajorDel) / opt->fRmajorDel) * + ceil((opt->fRminorMax - opt->fRminorMin + 0.5*opt->fRminorDel) / opt->fRminorDel) ; if (makeVectors) { - opt->fRMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); - opt->fRMajor = 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; - fRMajor += opt->fRMajorDel) { - for (float fRMinor = opt->fRMinorMin; fRMinor < opt->fRMinorMax + 0.5*opt->fRMinorDel; - fRMinor += opt->fRMinorDel) { - opt->fRMinor->data.F32[i] = fRMinor; - opt->fRMajor->data.F32[i] = fRMajor; + 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) { + opt->fRminor->data.F32[i] = fRminor; + opt->fRmajor->data.F32[i] = fRmajor; i++; } @@ -627,7 +679,8 @@ psAssert(i == opt->numTrials, "Something's wrong with my loop got %d entries expected %d", i, opt->numTrials); } else { - opt->fRMinor = NULL; - opt->fRMajor = NULL; - } + opt->fRminor = NULL; + opt->fRmajor = NULL; + } +#endif return true; @@ -640,6 +693,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)) { @@ -647,11 +700,18 @@ 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") - CHECKVAL(options, thisReadoutsOptions, fRMajorMax, "mismatched fRMajorMax") - CHECKVAL(options, thisReadoutsOptions, fRMajorDel, "mismatched fRMajorDel") - CHECKVAL(options, thisReadoutsOptions, fRMinorMin, "mismatched fRMinorMin") - CHECKVAL(options, thisReadoutsOptions, fRMinorMax, "mismatched fRMinorMax") - CHECKVAL(options, thisReadoutsOptions, fRMinorDel, "mismatched fRMinorDel") + CHECKVAL(options, thisReadoutsOptions, fRmajorMin, "mismatched fRmajorMin") + CHECKVAL(options, thisReadoutsOptions, fRmajorMax, "mismatched fRmajorMax") + CHECKVAL(options, thisReadoutsOptions, fRmajorDel, "mismatched fRmajorDel") + CHECKVAL(options, thisReadoutsOptions, fRminorMin, "mismatched fRminorMin") + CHECKVAL(options, thisReadoutsOptions, fRminorMax, "mismatched fRminorMax") + CHECKVAL(options, thisReadoutsOptions, fRminorDel, "mismatched fRminorDel") +#endif return true;