Index: trunk/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36633) +++ trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36662) @@ -10,10 +10,12 @@ psF32 fRminorMax; psF32 fRminorDel; - psVector *rMajor; - psVector *rMinor; + psVector *fRMajor; + psVector *fRMinor; psArray *zeroPt; psArray *exptime; psArray *cffTable; } galaxyShapeOptions; + + static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options); @@ -47,11 +49,8 @@ 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; @@ -91,6 +90,4 @@ psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime; fluxScaleFactor->data.F32[index] = scaleFactor; - -// readouts->data[index] = readout; // find detections @@ -184,5 +181,4 @@ psVector *sumInvSig2 = NULL; psVector *numerator = NULL; -// psVector *tmp = NULL; psF32 totalNPix = 0; long vectorLength = 0; @@ -220,8 +216,11 @@ } + // The only parameters that we are summarizing currently are galaxy shapes. + // Continue for inputs that do not have modelEXT + if (!source->modelEXT) continue; + 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); @@ -237,6 +236,7 @@ source->galaxyFits->chisq->n, vectorLength); - psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; + + totalNPix += source->galaxyFits->nPix; for (int k = 0; k < vectorLength; k++) { @@ -245,27 +245,10 @@ 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; } @@ -276,23 +259,53 @@ 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); + outSrc->galaxyFits->nPix = totalNPix; + psVector *fluxVec = + outSrc->galaxyFits->Flux = psVectorRecycle(outSrc->galaxyFits->Flux, vectorLength, PS_TYPE_F32); + psVector *dFluxVec = + outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32); + psVector *chisqVec = + outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32); + + // Get nominal radii paramters from the cff table + 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; + } + psMetadata *cffRow = options->cffTable->data[outSrc->seq]; + if (!cffRow) { + psError (PS_ERR_PROGRAMMING, true, "NO cff data for object %d", outSrc->seq); + 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", outSrc->seq); + 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", outSrc->seq); + psFree(outSrc); + return NULL; + } + + // fill the summary galaxyFits vectors + + int min_k = -1; + psF32 minChisq = NAN; + 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,51 +313,194 @@ 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_j >= 0 && isfinite(minChisq) && outSrc->modelEXT) { + + if (min_k >= 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; - } - 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; - } - 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 - } + +#define NUM_TRIALS_TO_FIT 9 + psVector *major = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64); + psVector *minor = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64); + psVector *chisq = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64); + psVector *flux = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64); + psVector *dFlux = psVectorAllocEmpty(NUM_TRIALS_TO_FIT, PS_TYPE_F64); + + // fractional radii with the lowest chisq + psF32 fRMajorBest = options->fRMajor->data.F32[min_k]; + psF32 fRMinorBest = options->fRMinor->data.F32[min_k]; + // XXX: use a recipe parameter + psF32 maxDelta = 0.11; + 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(options->fRMajor->data.F32[k] - fRMajorBest) < maxDelta && + fabs(options->fRMinor->data.F32[k] - fRMinorBest) < maxDelta) { + + if (isfinite(chisqVec->data.F32[k]) && + isfinite(fluxVec->data.F32[k]) && + isfinite(dFluxVec->data.F32[k])) { + + // compute major and minor radius vectors are in pixel space from nominal and trial fractions + psF64 thisMajor = R_MAJOR * options->fRMajor->data.F32[k]; + if (!isfinite(majorMin) || thisMajor < majorMin) { + majorMin = thisMajor; + } + if (!isfinite(majorMax) || thisMajor > majorMax) { + majorMax = thisMajor; + } + psVectorAppend(major, thisMajor); + + psF64 thisMinor = R_MINOR * options->fRMinor->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_TO_FIT) { + fprintf(stderr, "only found %ld good points wanted %d for seq: %d\n", + major->n, NUM_TRIALS_TO_FIT, outSrc->seq); + } + + // rescale data + 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); + 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; + + 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; + + psF64 chisq0 = psPolynomial2DEval(chisqFit, major0, minor0) * chisqRange + chisqMin; + + // flux and dFlux values from entry with lowest chisq + psF64 fluxBest = fluxVec->data.F32[min_k]; + psF64 dFlux0 = dFluxVec->data.F32[min_k]; + + bool goodFluxFit = psVectorFitPolynomial2D(fluxFit, NULL, 0xFF, flux, NULL, major, minor); + // flux from fit + psF64 flux0; + if (goodFluxFit) { + flux0 = psPolynomial2DEval(fluxFit, major0, minor0); + } else { + flux0 = fluxBest; + } + + // compare to the minimum value in the vector + + 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); + + psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type); + +#ifdef USE_MINIMUM_CHISQ + axes.major = R_MAJOR * fRMajorBest; + axes.minor = R_MINOR * fRMinorBest; + outSrc->modelEXT->chisq = minChisq; + flux0 = fluxBest; +#else + axes.major = major0; + axes.minor = minor0; + outSrc->modelEXT->chisq = chisq0; +#endif + pmPSF_AxesToModel (outSrc->modelEXT->params->data.F32, axes, outSrc->modelEXT->type); + + outSrc->modelEXT->mag = -2.5 * log10(flux0); + outSrc->modelEXT->magErr = dFlux0 / flux0; // 1 / SN + outSrc->modelEXT->flags |= PM_MODEL_STATUS_FITTED; + } else { +#ifdef notdef + psErrorStackPrint(stderr, NULL); +#endif + outSrc->modelEXT->flags |= PM_MODEL_STATUS_NONCONVERGE; + + 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); + psErrorClear(); + // psFree(outSrc->galaxyFits); + + // XXX drop the model since it is bogus? Perhaps we do this or use the nominal values + // and add a flag +#ifdef notnow + psFree(outSrc->modelEXT); + outSrc->modelEXT = 0; + psFree(outSrc->modelFits); + outSrc->modelFits = 0; +#endif + } + 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); + psFree(chisqFit); + psFree(major); + psFree(minor); + psFree(chisq); + psFree(flux); + psFree(dFlux); + } + } else { +#ifdef notnow + // drop the extended source information so that we can tell which objects that we fit + psFree(outSrc->modelEXT); + psFree(outSrc->modelFits); +#endif } } @@ -378,6 +534,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 +541,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,6 +548,6 @@ 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; } Index: trunk/psphot/src/psphotLensing.c =================================================================== --- trunk/psphot/src/psphotLensing.c (revision 36633) +++ trunk/psphot/src/psphotLensing.c (revision 36662) @@ -193,4 +193,8 @@ // set the max order (0 = constant) which the number of psf stars can support: // we require only 3 stars for n = 0, increase stars / cell for higher order + if (Npsf < 3) { + psLogMsg ("psphot", PS_LOG_INFO, "no sources passed selection criteria, skipping lensing psf trends"); + return true; + } int LensingTrendOrder = 0; if (Npsf >= 16) LensingTrendOrder = 1; // 4 cells