Changeset 36410

Timestamp:

Dec 24, 2013, 1:28:04 PM (13 years ago)

Author:

bills

Message:

progress on fulLforceSummary implementation. in psphotFullForce set
modelEXT to the nominal model so that those paramters are available in fullForceSummary

Location:

branches/bills_branches/bills_201312/psphot/src

Files:

• psphotFullForceSummary.c (modified) (6 diffs)
• psphotFullForceSummaryReadout.c (modified) (3 diffs)
• psphotGalaxyShape.c (modified) (3 diffs)

branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummary.c

@@ -140 +140 @@
     psFree(outputFPA);                        // Drop reference
     if (!output) {
-        psError(psErrorCodeLast(), false, _("Unable to generate output file from PPSTACK.OUTPUT"));
+        psError(psErrorCodeLast(), false, _("Unable to generate output file from PSPHOT.OUTPUT"));
         return false;
     }
@@ -183 +183 @@
         return NULL;
     }
-#ifdef notdef
-    pmHDU *lastHDU = NULL;              // Last HDU updated
-#endif
 
     // files associated with the science image
@@ -196 +193 @@
 
     while ((chip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
-        psLogMsg ("psphotFullForceSummary", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
+//        psLogMsg ("psphotFullForceSummary", 4, "Chip %d: %x %x\n", view->chip, chip->file_exists, chip->process);
         if (! chip->process || ! chip->file_exists) { continue; }
 
@@ -233 +230 @@
         // there is now only a single chip (multiple readouts?). loop over it and process
         while ((cell = pmFPAviewNextCell (view, input->fpa, 1)) != NULL) {
-            psLogMsg ("psphotFullForceSummary", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+ //           psLogMsg ("psphotFullForceSummary", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
 
             // process each of the readouts
             while ((readout = pmFPAviewNextReadout (view, input->fpa, 1)) != NULL) {
-                psLogMsg ("psphotFullForceSummary", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
+//                psLogMsg ("psphotFullForceSummary", 6, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process);
                 if (! readout->data_exists) { continue; }
 
-#ifdef notdef
-                pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
-                if (hdu && hdu != lastHDU) {
-                    // XXX: need to do this
-                    // psphotVersionHeaderFull(hdu->header);
-                    lastHDU = hdu;
-                }
-#endif
                 if (!psphotFullForceSummaryReadout(config, view)) {
                     ESCAPE ("failure in psphotFullForceSummaryReadout");
@@ -257 +246 @@
 
 
-    // If these keywords are not set we get a warning message on output.
+    // If these keywords are not set we get a warning message on output. Since we are not copying the input header
+    // and do not have an image pmFPA stuff doesn't have values for these
+    // XXX: What other keywords and concepts should be set (or copied)?
     int numCols = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS1");
     int numRows = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS2");
     psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS1", PS_META_REPLACE, "", numCols);
     psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS2", PS_META_REPLACE, "", numRows);
-    // XXX: What other keywords and concepts should be set (or copied)?
+
+    // XXX: Also add psphot version information
 
     if (!pmFPAfileIOChecks (config, view, PM_FPA_AFTER)) ESCAPE ("failed pmFPAfileIOChecks for FPA in psphot.");
@@ -270 +262 @@
 
     psFree (view);
+
     return true;
 }

branches/bills_branches/bills_201312/psphot/src/psphotFullForceSummaryReadout.c

@@ -1 +1 @@
 #include "psphotInternal.h"
+
+
+typedef struct {
+    int     numTrials;
+    psF32   fRmajorMin;
+    psF32   fRmajorMax;
+    psF32   fRmajorDel;
+    psF32   fRminorMin;
+    psF32   fRminorMax;
+    psF32   fRminorDel;
+    psVector    *rMajor;
+    psVector    *rMinor;
+    psArray *zeroPt;
+    psArray *exptime;
+} galaxyShapeOptions;
+
+static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *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);
+
 
 bool psphotFullForceSummaryReadout (pmConfig *config, const pmFPAview *view) {
@@ -15 +37 @@
         outputReadout = pmReadoutAlloc(outputCell);
     }
-        
+
+    // Get the exposure parameters for the output from recipe and set them on the output
+    psF32 outputZeroPoint = psMetadataLookupF32(&status, recipe, "PSPHOT.FULLFORCE.ZERO_PT");
+    psF32 outputExptime = psMetadataLookupF32(&status, recipe, "PSPHOT.FULLFORCE.EXPTIME");
+
+    psMetadataAddF32(output->fpa->concepts, PS_LIST_TAIL, "FPA.ZP", PS_META_REPLACE, "Magnitude zero point",
+        outputZeroPoint);
+    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
-    pmReadout *savedReadout = NULL;
+//    psArray *readouts = psArrayAlloc(num);
+    psVector *fluxScaleFactor = psVectorAlloc(num, PS_TYPE_F32);
+    galaxyShapeOptions options;
+    psArray *objects = NULL;
     for (int index = 0; index < num; index++) {
         // find the currently selected readout
@@ -27 +63 @@
         pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
         psAssert (readout, "missing readout?");
-        if (savedReadout == NULL) {
-            // save this readout for help in setting some concepts below
-            savedReadout = readout;
-        }
+
+        if (index == 0) {
+            // Get the galaxy shape recipe values, from the analysis if present
+            // or from the recipe if not
+            if (!setOptions(&options, readout, recipe, true)) {
+                psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options.");
+                return false;
+            }
+        } else { 
+            // Make sure that this input was created with the same galaxy shapes recipe
+            if (!checkOptions(&options, readout, recipe)) {
+                psError (PS_ERR_UNKNOWN, false, "galaxy shape options do not match for input %d", index);
+                return false;
+            }
+        }
+
+        // look up zero point
+        psF32 zero_point = psMetadataLookupF32(&status, readout->parent->parent->parent->concepts, "FPA.ZP");
+        psF32 exptime = psMetadataLookupF32(&status, readout->parent->parent->parent->concepts, "FPA.EXPOSURE");
+
+        psF32 scaleFactor = pow(10, 0.4 * (outputZeroPoint - zero_point)) * outputExptime / exptime;
+        fluxScaleFactor->data.F32[index] = scaleFactor;
+
+//        readouts->data[index] = readout;
 
         // find detections
         pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-
-        if (index == 0) {
-            // As a test, just copy readout zero's detections to the output
-            if (!psMetadataAddPtr (outputReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", detections)) {
-                psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
-                return false;
-            }
-
-            outputReadout->data_exists = true;
-            outputReadout->parent->data_exists = true;
-            outputReadout->parent->parent->data_exists = true;
-        }
-    }
+        psAssert (detections, "missing detections?");
+        psArray *sources = detections->allSources;
+        psAssert (sources, "missing sources?");
+        sources = psArraySort (sources, pmSourceSortBySeq);
+
+        if (objects == NULL) {
+            pmSource *lastSource = sources->data[sources->n - 1];
+            psAssert(lastSource, "last source is null!");
+            objects = psArrayAlloc(lastSource->seq + 1);
+        }
+
+        for (int i=0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            source->imageID = index;
+            findObjectForSource(&objects, source);
+        }
+    }
+
+    pmDetections *outputDetections = pmDetectionsAlloc();
+    if (!psMetadataAddPtr (outputReadout->analysis, PS_LIST_TAIL, "PSPHOT.DETECTIONS", 
+            PS_META_REPLACE | PS_DATA_UNKNOWN, "psphot detections", outputDetections)) {
+        psError (PSPHOT_ERR_CONFIG, false, "problem saving detections on readout");
+        return false;
+    }
+
+    psArray *outputSources = psArrayAllocEmpty (objects->n);
 
     // Loop over objects and compute the summaries
-    // save the summaries on the output readout->analysis
-
+    long nObjects = 0;
+    long nSources = 0;
+    for (int i=0 ; i<objects->n; i++) {
+        pmPhotObj *obj = objects->data[i];
+        if (!obj) continue;
+
+        ++nObjects;
+        pmSource *source = psphotFullForceSummarizeObject(config, obj, fluxScaleFactor, &options);
+        if (source) {
+            psArrayAdd (outputSources, 100, source);
+            ++nSources;
+        }
+    }
+
+    psLogMsg("psphot", PS_LOG_INFO, "constructed %ld output sources, from %ld objects.", nSources, nObjects);
+
+    psFree(fluxScaleFactor);
+
+    if (nSources) {
+        // We have data
+        outputDetections->allSources = outputSources;
+        outputReadout->data_exists = true;
+        outputReadout->parent->data_exists = true;
+        outputReadout->parent->parent->data_exists = true;
+    } else {
+        // XXX: tooo set a quality or fault code
+        return false;
+    }
 
     return true;
 }
+
+static pmPhotObj *findObjectForSource(psArray **pObjects, pmSource *source) {
+    int seq = source->seq;
+
+    psArray *objects = *pObjects;
+    if (seq >= objects->n) {
+        // We need to expand the object array. Kind of suprising.
+        objects = *pObjects = psArrayRealloc(objects, seq+1);
+    }
+
+    // Look up object for this seq
+    pmPhotObj *obj = objects->data[seq];
+    if (!obj) {
+        // not found allocate one
+        obj = pmPhotObjAlloc();
+        objects->data[seq] = obj;
+    }
+    pmPhotObjAddSource(obj, source);
+
+    return obj;
+}
+
+static pmSource *psphotFullForceSummarizeObject(pmConfig *config, pmPhotObj *obj, psVector *fluxScaleFactor, galaxyShapeOptions *options) {
+
+    pmSource *outSrc = NULL;
+
+    psVector *sumWeightedFlux = NULL;
+    psVector *sumInvSig2 = NULL;
+    psVector *numerator   = NULL;
+//    psVector *tmp = NULL;
+    psF32   totalNPix = 0;
+    long    vectorLength = 0;
+
+    for (int i=0; i < obj->sources->n; i++) {
+        pmSource *source = obj->sources->data[i];
+
+        // XXX: get cut from recipe
+        if (source->pixWeightNotPoor < .9) continue;
+
+        // For now just start the output source as a copy of the first input source that makes cuts
+        if (!outSrc) {
+            outSrc = pmSourceCopy(source);
+            // This copies 
+            //  the peak 
+            //  the moments which are mostly nan except for Mrf Mx, My, and some of the kron parameters
+            //
+            // type, mode, flags
+            // magnitudes
+            outSrc->imageID = 0;
+            outSrc->seq = source->seq;
+
+            if (source->modelPSF) {
+                outSrc->modelPSF = psMemIncrRefCounter(source->modelPSF);
+            }
+            if (source->extpars) {
+                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];
+
+            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
+            }
+
+            psFree(numerator);
+            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) {
+                // copy the best fit params to the model
+                psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
+                axes.major *= options->rMajor->data.F32[min_j];
+                axes.minor *= 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
+            }
+        }
+    }
+
+    return outSrc;
+}
+
+#define GETVAL(member, key) \
+    opt->member = psMetadataLookupF32(&status, md, key); \
+    if (!status) { \
+        psError (PSPHOT_ERR_CONFIG, true, "failed to looup value for %s in %s", key, \
+            useAnalysis ? "readout->analysis" : "recipe"); \
+    }
+
+static bool setOptions(galaxyShapeOptions *opt, pmReadout *readout, psMetadata *recipe, bool makeVectors) {
+    bool status;
+    bool useAnalysis;
+
+    psMetadataLookupF32(&useAnalysis, readout->analysis, "GALAXY_SHAPES_FR_MAJOR_MIN");
+    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) ;
+
+    if (makeVectors) {
+        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;
+                fRmajor += opt->fRmajorDel) {
+            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;
+                i++;
+            }
+        }
+        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;
+    }
+        
+    return true;
+}
+
+#define CHECKVAL(left, right, val, message) \
+    if (left->val != right.val) { \
+        psError (PSPHOT_ERR_CONFIG, true, message); \
+        return false; \
+    }
+
+static bool checkOptions(galaxyShapeOptions *options, pmReadout *readout, psMetadata *recipe) {
+    galaxyShapeOptions thisReadoutsOptions;
+    
+    if (!setOptions(&thisReadoutsOptions, readout, recipe, false)) {
+        psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options for readout");
+        return false;
+    }
+    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")
+
+    return true;
+}

branches/bills_branches/bills_201312/psphot/src/psphotGalaxyShape.c

@@ -283 +283 @@
     }
 
+#ifdef SAVE_BEST_MODEL
+    // Save model with smallest chisq
     if (isfinite(chisqBest)) {
+#else 
+    // Save model with nominal parameters
+    {
+        fRmajorBest = 1;
+        fRminorBest = 1;
+#endif
         // now save the best fitting model as the source's extended model
         psEllipseAxes testAxes = guessAxes;
@@ -293 +301 @@
         psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
 
+        // Replace modelEXT with this model
+        // XXX: only do this if the model is good
         psFree (source->modelEXT);
 
@@ -299 +309 @@
         source->mode |= PM_SOURCE_MODE_EXTMODEL;
         source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-        // adjust the window so the subtraction covers the faint wings
-        // psphotSetRadiusMoments(&fitRadius, &windowRadius, readout, source, markVal);
 
         // cache the model flux