Changeset 36588

Timestamp:

Mar 14, 2014, 1:59:13 PM (12 years ago)

Author:

bills

Message:

in psphotGalaxyShape save best model as modelEXT instead of nominal.
Pass cff to psphotFullForceSummary and use it to obtain the nominal parameter values

Location:

tags/ipp-20140114/psphot/src

Files:

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

tags/ipp-20140114/psphot/src/psphotFullForceSummary.c

@@ -78 +78 @@
         usage();
     }
+    if ((N = psArgumentGet(argc, argv, "-cff"))) {
+        if (argc <= N+1) {
+          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
+          usage();
+        }
+        psArgumentRemove(N, &argc, argv);
+
+        // We read CFF file as table not as sources. All we want are the nominal parameters input to psphotFullForce
+        psMetadataAddStr(config->arguments, PS_LIST_TAIL, "CFF_FILE", 0, "CFF file with nominal parameters", argv[N]);
+        psArgumentRemove(N, &argc, argv);
+    } else {
+        psErrorStackPrint(stderr, "-cff must be supplied.");
+        usage();
+    }
  
     if (argc < 2) {
@@ -148 +162 @@
     }
     output->save = true;
+
+    // We read cff file as a table not as sources
+    psString cffName = psMetadataLookupStr(NULL, config->arguments, "CFF_FILE");
+    if (!cffName) {
+        psError(PSPHOT_ERR_CONFIG, true, "CFF_FILE is missing from arguments");
+        return false;
+    }
+    psString resolvedName = pmConfigConvertFilename(cffName, config, false, false);
+    if (!resolvedName) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to resolve CFF_FILE %s", cffName);
+        return false;
+    }
+    psFits *fits = psFitsOpen(resolvedName, "r");
+    if (!fits) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to open fits CFF_FILE %s", resolvedName);
+        return false;
+    }
+    psArray *inTable = psFitsReadTable(fits);
+    if (!inTable) {
+        psError(PSPHOT_ERR_CONFIG, false, "failed to read cff fits table from CFF_FILE %s", resolvedName);
+        return false;
+    }
+    psFitsClose(fits);
+
+    // Convert to an array indexed by ID
+    psArray *sortedTable = psArrayAlloc(4*inTable->n);
+    for (int i=0; i<inTable->n; i++) {
+        psMetadata *row = inTable->data[i];
+        psS32 ID = psMetadataLookupS32(&status, row, "ID");
+        if (ID >= sortedTable->n) {
+            sortedTable = psArrayRealloc(sortedTable, 2*ID);
+        }
+        if (sortedTable->data[ID]) {
+            psError(PSPHOT_ERR_CONFIG, true, "Duplicate row with ID %d", ID);
+            return false;
+        }
+        sortedTable->data[ID] = psMemIncrRefCounter(row);
+    }
+    psMetadataAddArray(config->arguments, PS_LIST_TAIL, "CFF_TABLE", PS_META_REPLACE, "cff table", sortedTable);
+    psFree(inTable);
+    psFree(sortedTable);
 
     return true;

tags/ipp-20140114/psphot/src/psphotFullForceSummaryReadout.c

@@ -14 +14 @@
     psArray *zeroPt;
     psArray *exptime;
+    psArray *cffTable;
 } galaxyShapeOptions;
 
@@ -69 +70 @@
             if (!setOptions(&options, readout, recipe, true)) {
                 psError (PS_ERR_UNKNOWN, false, "problem determining galaxy shape options.");
+                return false;
+            }
+            options.cffTable = psMetadataLookupPtr(NULL, config->arguments, "CFF_TABLE");
+            if (!options.cffTable) {
+                psError (PS_ERR_UNKNOWN, true, "Cannot find cff table in arguments.");
                 return false;
             }
@@ -309 +315 @@
                 // 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];
+                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; 

tags/ipp-20140114/psphot/src/psphotGalaxyShape.c

@@ -274 +274 @@
 
             int i = source->galaxyFits->chisq->n - 1;
-            float thisChisq = source->galaxyFits->chisq->data.F32[i];
-            if (isfinite(thisChisq) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
-                chisqBest = thisChisq;
-                fRmajorBest = fRmajor;
-                fRminorBest = fRminor;
+            float flux = source->galaxyFits->Flux->data.F32[i];
+            if (isfinite(flux)) {
+                float thisChisq = source->galaxyFits->chisq->data.F32[i];
+                if (isfinite(thisChisq) && isfinite(flux) && (!isfinite(chisqBest) || thisChisq < chisqBest)) {
+                    chisqBest = thisChisq;
+                    fRmajorBest = fRmajor;
+                    fRminorBest = fRminor;
+                }
             }
+            // reset I0 to avoid potential problems on the next iteration
+            PAR[PM_PAR_I0] = 1.0;
         }
     }
 
+#define SAVE_BEST_MODEL
 #ifdef SAVE_BEST_MODEL
     // Save model with smallest chisq
@@ -301 +307 @@
         psphotGalaxyShapeSource (pcm, source, maskVal, psfSize, false);
 
-        // Replace modelEXT with this model
-        // XXX: only do this if the model is good
-        psFree (source->modelEXT);
-
-        source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
-        source->type = PM_SOURCE_TYPE_EXTENDED;
-        source->mode |= PM_SOURCE_MODE_EXTMODEL;
-        source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
-
-        // cache the model flux
-        pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+        // Replace modelEXT with this model, if the model is good
+        if (isfinite(PAR[PM_PAR_I0])) {
+            psFree (source->modelEXT);
+
+            source->modelEXT = psMemIncrRefCounter (pcm->modelConv);
+            source->type = PM_SOURCE_TYPE_EXTENDED;
+            source->mode |= PM_SOURCE_MODE_EXTMODEL;
+            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+
+            // cache the model flux
+            pmPCMCacheModel (source, maskVal, psfSize, fitOptions->nsigma);
+        }
     }