Index: trunk/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36663) +++ trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 36664) @@ -10,12 +10,10 @@ psF32 fRminorMax; psF32 fRminorDel; - psVector *fRMajor; - psVector *fRMinor; + psVector *rMajor; + psVector *rMinor; psArray *zeroPt; psArray *exptime; psArray *cffTable; } galaxyShapeOptions; - - static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options); @@ -49,8 +47,11 @@ 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; @@ -90,4 +91,6 @@ psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime; fluxScaleFactor->data.F32[index] = scaleFactor; + +// readouts->data[index] = readout; // find detections @@ -181,4 +184,5 @@ psVector *sumInvSig2 = NULL; psVector *numerator = NULL; +// psVector *tmp = NULL; psF32 totalNPix = 0; long vectorLength = 0; @@ -216,11 +220,8 @@ } - // 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); @@ -236,7 +237,6 @@ source->galaxyFits->chisq->n, vectorLength); + psF32 scaleFactor = fluxScaleFactor->data.F32[source->imageID]; - - totalNPix += source->galaxyFits->nPix; for (int k = 0; k < vectorLength; k++) { @@ -245,10 +245,27 @@ 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; } @@ -259,53 +276,23 @@ if (vectorLength) { outSrc->galaxyFits = pmSourceGalaxyFitsAlloc(); - 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; - + 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++) { - 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; - } + 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 } @@ -313,194 +300,51 @@ psFree(sumInvSig2); psFree(sumWeightedFlux); - - if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) { +// 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) { // copy the best fit params to the model - -#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]); - } - } + 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; } - - // 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); + 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; } - - // rescale data - psF64 majorRange = majorMax - majorMin; - if (majorRange == 0) { - majorRange = 1; + 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; } - psF64 minorRange = minorMax - minorMin; - if (minorRange == 0) { - minorRange = 1; + 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; } - 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 + 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 + } } } @@ -534,6 +378,6 @@ if (makeVectors) { - opt->fRMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); - opt->fRMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->rMinor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); + opt->rMajor = psVectorAlloc(opt->numTrials, PS_TYPE_F32); int i = 0; for (float fRmajor = opt->fRmajorMin; fRmajor < opt->fRmajorMax + 0.5*opt->fRmajorDel; @@ -541,6 +385,6 @@ 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; + opt->rMinor->data.F32[i] = fRminor; + opt->rMajor->data.F32[i] = fRmajor; i++; } @@ -548,6 +392,6 @@ 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->rMinor = NULL; + opt->rMajor = NULL; }