Changeset 36085 for trunk/psModules/src/objects/pmPCM_MinimizeChisq.c

Timestamp:

Aug 31, 2013, 5:55:16 AM (13 years ago)

Author:

eugene

Message:

merge changes from EAM branch (pmModel_CentralPixel functions for sersic-like model support; add interactive mode to PCM LMM fitting; harder PCM LMM damping (nu scaled by 3 not 2); add EXT_AND SKY and SHAPE fitting modes, include sky in INDEX-only fits; ifdefs for trace-only values; allow positions to go to 100,000; major re-work of central pixel & flux for sersic-like models; add pmPCMMakeModel function to generate the flux for an arbitrary model

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/pmPCM_MinimizeChisq.c (modified) (8 diffs)

trunk/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"
 
@@ -91 +96 @@
     psF32 lambda = 0.001;
     psF32 dLinear = 0.0;
-    psF32 nu = 2.0;
+    psF32 nu = 3.0;
 
 # if (USE_FFT && PRE_CONVOLVE)
@@ -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;
+                    int Nscan = fscanf (stdin, "%d %f", &Npar, &value);
+                    if (Nscan != 2) {
+                      fprintf (stderr, "scan failure\n");
+                    }
+                    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 +267 @@
         // 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);
@@ -215 +280 @@
           case 0:
             if (rho >= -1e-6) {
-                lambda *= 0.25;
+                lambda *= 0.1;
             } else {
                 lambda *= 10.0;
@@ -234 +299 @@
             if (rho > 0.0) {
                 lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
-                nu = 2.0;
+                nu = 3.0;
             } else {
                 lambda *= nu;
-                nu *= 2.0;
+                nu *= 3.0;
             }
             break;
@@ -474 +539 @@
     // 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