Changeset 35961 for branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c

Timestamp:

Aug 15, 2013, 5:56:56 PM (13 years ago)

Author:

eugene

Message:

more work on the central pixel optimizations -- perhaps not needed (not so expensive?); add some interactive support for PCM chisq fitting; EXP and DEV are for the moment using subdivided central pixels, but this is perhaps too slow?; turn on sky fitting for the PCM model fitting

File:

• branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c (modified) (5 diffs)

branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c

@@ -42 +42 @@
 #include "pmPCMdata.h"
 
+# define SAVE_IMAGES 0
+# if (SAVE_IMAGES) 
+int psphotSaveImage (psMetadata *header, psImage *image, char *filename);
+# endif
+
 # define FACILITY "psModules.objects"
 
@@ -130 +135 @@
         }
 
+        char key[10]; // used for interactive responses
+        bool testValue = false;
+
         // set a new guess for Alpha, Beta, Params
         if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
+            if (min->isInteractive) {
+                fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] ");
+                if (!fgets(key, 8, stdin)) {
+                    psWarning("Unable to read option");
+                }
+                switch (key[0]) {
+                  case 'n':
+                  case 'N':
+                    done = true;
+                    break;
+                  case 'y':
+                  case 'Y':
+                  case '\n':
+                    lambda *= 10.0;
+                    continue;
+                  default:
+                    lambda *= 10.0;
+                    continue;
+                }
+                if (done) break;
+            }
             min->iter ++;
             if (min->iter >=  min->maxIter) break;
@@ -138 +167 @@
         }
 
+        if (min->isInteractive) {
+            p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: ");
+            fprintf (stdout, "last chisq : %f\n", min->value);
+            bool getOptions = true;
+            while (getOptions) {
+                fprintf (stdout, "options: (m)odify, (g)o, (q)uit: ");
+                if (!fgets(key, 8, stdin)) {
+                    psWarning("Unable to read option");
+                }
+                switch (key[0]) {
+                  case 'm':
+                  case 'M':
+                    testValue = TRUE;
+                    fprintf (stdout, "enter (Npar) (value): ");
+                    int Npar = 0;
+                    float value= 0;
+                    fscanf (stdin, "%d %f", &Npar, &value);
+                    Params->data.F32[Npar] = value;
+                    break;
+                  case 'g':
+                  case 'G':
+                  case '\n':
+                    getOptions = false;
+                    break;
+                  default:
+                    done = true;
+                    break;
+                }
+                fprintf (stderr, "foo\n");
+            }
+            if (done) break;
+        }
+            
         // dump some useful info if trace is defined
         if (psTraceGetLevel(FACILITY) >= 6) {
@@ -202 +264 @@
         // XXX : Madsen gives suggestion for better use of rho
         // rho is positive if the new chisq is smaller
-        if (rho >= -1e-6) {
+        if (testValue || (rho >= -1e-6)) {
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
@@ -474 +536 @@
     // XXX TEST : SAVE IMAGES
 # if (SAVE_IMAGES)
-    psphotSaveImage (NULL, pcm->psf->image, "psf.fits");
-    psphotSaveImage (NULL, pcm->modelFlux, "model.fits");
-    psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits");
-    psphotSaveImage (NULL, source->pixels, "obj.fits");
-    psphotSaveImage (NULL, source->maskObj, "mask.fits");
-    psphotSaveImage (NULL, source->variance, "variance.fits");
+    static int Npass = 0;
+    char name[128]; 
+    snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name);
+    snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name);
+    snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name);
+    snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name);
+    snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name);
+    snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name);
+    for (int n = 0; n < pcm->dmodelsFlux->n; n++) {
+        psImage *dmodelConv = pcm->dmodelsConvFlux->data[n];
+        if (!dmodelConv) continue;
+        snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name);
+    }
+    Npass ++;
 # endif