Changeset 36718 for trunk/psphot/src/psphotFullForceSummaryReadout.c

Timestamp:

May 1, 2014, 4:36:26 PM (12 years ago)

Author:

bills

Message:

Modify psphotFullForce and psphotFullForceSummary to support forcing
multiple extended model fits per source.

File:

• trunk/psphot/src/psphotFullForceSummaryReadout.c (modified) (10 diffs)

trunk/psphot/src/psphotFullForceSummaryReadout.c

@@ -14 +14 @@
     psArray *zeroPt;
     psArray *exptime;
-    psArray *cffTable;
+    psArray *cffTables;
 } galaxyShapeOptions;
 
@@ -71 +71 @@
                 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;
@@ -178 +178 @@
     pmSource *outSrc = NULL;
 
-    psVector *sumWeightedFlux = NULL;
-    psVector *sumInvSig2 = NULL;
-    psVector *numerator   = NULL;
-    psF32   totalNPix = 0;
-    long    vectorLength = 0;
+
+    pmSource *source0 = obj->sources->data[0];
+    int ID = source0->seq;
+
+    // Find the parameters from the cff tables
+    int nModels = 0;
+    // skip table 0 which is modelType -1 not extended
+    for (int i=1; i<options->cffTables->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
@@ -210 +220 @@
                 outSrc->extpars =  psMemIncrRefCounter(source->extpars);
             }
-            if (source->modelEXT) {
-                outSrc->modelEXT =  psMemIncrRefCounter(source->modelEXT);
-            }
-            if (source->modelFits) {
-                outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
-            }
-        }
-
-        // The only parameters that we are summarizing currently are galaxy shapes. 
-        // Continue for inputs that do not have modelEXT which contains the model paramters
-        if (!source->modelEXT) continue;
-
-        if (source->galaxyFits && isfinite(source->galaxyFits->nPix) && source->galaxyFits->chisq->n) {
-            if (numerator == NULL) {
-                // first source with galaxyFits allocate accumulators
-                vectorLength = source->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);
-            }
-
-            // Die if the lengths of the vectors is not the same in all sources
-            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];
-
-            totalNPix += source->galaxyFits->nPix;
-
-            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;
-
-                numerator->data.F32[k] += chisq * source->galaxyFits->nPix;
-
-                psF32 invSig2 = 1.0 / (dFlux * dFlux);
-                sumInvSig2->data.F32[k] += invSig2;
-
-                sumWeightedFlux->data.F32[k] += flux * invSig2;
-            }
-        }
-    }
-
-    if (outSrc) {
+            outSrc->modelEXT = psFree(outSrc->modelEXT);
+        }
+        if (nModels && !outSrc->modelFits && source->modelFits && source->modelFits->n == nModels) {
+            outSrc->modelFits =  psMemIncrRefCounter(source->modelFits);
+        }
+    }
+
+    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;
+        psF32   totalNPix = 0;
+        long    vectorLength = 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; j<source->galaxyFits->n; j++) {
+                galaxyFits = source->galaxyFits->data[j];
+                if (galaxyFits && galaxyFits->modelType == modelType) break;
+                galaxyFits = NULL;
+            }
+            // end of proposed function
+
+            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);
+                }
+
+                // Die if the lengths of the vectors is not the same in all sources
+                psAssert(vectorLength == options->numTrials, "length of chisq vector %ld does not match options %d",
+                    vectorLength, options->numTrials);
+                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
-            outSrc->galaxyFits = pmSourceGalaxyFitsAlloc();
-            outSrc->galaxyFits->nPix  = totalNPix;
+            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 = 
-                outSrc->galaxyFits->Flux  = psVectorRecycle(outSrc->galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
+                galaxyFits->Flux  = psVectorRecycle(galaxyFits->Flux,  vectorLength, PS_TYPE_F32);
             psVector *dFluxVec = 
-                outSrc->galaxyFits->dFlux = psVectorRecycle(outSrc->galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
+                galaxyFits->dFlux = psVectorRecycle(galaxyFits->dFlux, vectorLength, PS_TYPE_F32);
             psVector *chisqVec = 
-                outSrc->galaxyFits->chisq = psVectorRecycle(outSrc->galaxyFits->chisq, vectorLength, PS_TYPE_F32);
-
-            // Find the nominal radii paramters in 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;
-            }
+                galaxyFits->chisq = psVectorRecycle(galaxyFits->chisq, vectorLength, PS_TYPE_F32);
+
 
             // fill the summary galaxyFits vectors and find the trial with the minimum value for chisq
@@ -319 +361 @@
             psFree(sumWeightedFlux);
 
-            if (min_k >= 0 && isfinite(minChisq) && outSrc->modelEXT) {
+            if (min_k >= 0 && isfinite(minChisq)) {
                 // copy the best fit params to the model
                 // fractional radii with the lowest chisq
-                psEllipseAxes axes = pmPSF_ModelToAxes(outSrc->modelEXT->params->data.F32, outSrc->modelEXT->type);
+                psEllipseAxes axes = pmPSF_ModelToAxes(model->params->data.F32, model->type);
 
                 // examine the params for the trial with minimum chisq.
@@ -472 +514 @@
                         axes.major = major0;
                         axes.minor = minor0;
-                        outSrc->modelEXT->chisq = chisq0;
+                        model->chisq = chisq0;
 #ifdef DUMPVECTORS
                         char fn[80];
@@ -487 +529 @@
 #endif
                     } else {
-                        outSrc->modelEXT->flags |= PM_MODEL_STATUS_NONCONVERGE;
+                        model->flags |= PM_MODEL_STATUS_NONCONVERGE;
 
 #ifdef PRINTVALS
@@ -496 +538 @@
                         // 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
                     }
                     psFree(chisqFit);
@@ -516 +550 @@
                     axes.major = R_MAJOR * fRMajorBest;
                     axes.minor = R_MINOR * fRMinorBest;
-                    outSrc->modelEXT->chisq = minChisq; 
+                    model->chisq = minChisq; 
                     flux0 = fluxBest;
 #ifdef PRINTVALS
@@ -525 +559 @@
                 // now save the model parameters in the model structure
 
-                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
+                pmPSF_AxesToModel (model->params->data.F32, axes, modelType);
+
+                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
-                outSrc->modelEXT->flags |= PM_MODEL_STATUS_FITTED;
+                model->flags |= PM_MODEL_STATUS_FITTED;
             }
         }