Changeset 36325

Timestamp:

Nov 28, 2013, 1:25:15 PM (13 years ago)

Author:

eugene

Message:

fix up grid search for sersic index

Location:

branches/eam_branches/ipp-20130904/psphot

Files:

• src/psphotSourceFits.c (modified) (9 diffs)
• test/tap_psphot_galaxygrid.pro (modified) (13 diffs)

branches/eam_branches/ipp-20130904/psphot/src/psphotSourceFits.c

@@ -653 +653 @@
     if (TIMING) { t4 = psTimerMark ("psphotFitPCM"); }
 
-    // psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
+    if (options.isInteractive) psphotPCMfitCheckSize (pcm, source, maskVal, psfSize);
     // if (pcm->modelConv->nIter == fitOptions->nIter) {
     //  psphotPCMfitRetry (pcm, source, &options, maskVal, markVal, psfSize);
@@ -791 +791 @@
 
 // 0.5 / n for (0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0)
-float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 1.666, 0.125, 0.10, 0.0833};
+float indexGuessInv[] = {1.00, 0.50, 0.333, 0.25, 0.166, 0.125, 0.10, 0.0833};
 float indexGuessR1q[] = {1.06, 1.19, 1.335, 1.48, 1.840, 2.290, 2.84, 3.5300};
 # define N_INDEX_GUESS_INV 8
@@ -902 +902 @@
                 rMin = reffGuess[j] / indexGuessR1q[i];
             }
-            // fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
-            // fprintf (stderr, "\n");
+            fprintf (stderr, "%d | %f %f %f %f | %f %f %f %f", i, indexGuessInv[i], reffGuess[j], Io, Chisq, sMin, rMin, iMin, xMin);
+            fprintf (stderr, "\n");
         }
     }
@@ -1006 +1006 @@
     psVector *Sidx = psVectorAllocEmpty (16, PS_TYPE_F32);
 
-    PAR[PM_PAR_7] = indexGuessInv[nStart];
+    float Sm = NAN, Sp = NAN, So = NAN;
+    if (nStart == 0) {
+        Sm = indexGuessInv[nStart];
+        So = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+        Sp = indexGuessInv[nStart + 1];
+    } else if (nStart == N_INDEX_GUESS_INV - 1) {
+        Sp = indexGuessInv[nStart];
+        So = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+        Sm = indexGuessInv[nStart - 1];
+    } else {
+        Sm = 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]);
+        So = indexGuessInv[nStart];
+        Sp = 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]);
+    }
+    
+    PAR[PM_PAR_7] = Sm;
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1012 +1027 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart < N_INDEX_GUESS_INV - 1) ? 0.5*(indexGuessInv[nStart + 1] + indexGuessInv[nStart]) : indexGuessInv[nStart - 1];
+    PAR[PM_PAR_7] = So;
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1018 +1033 @@
     psVectorAppend (chi2, model->chisqNorm);
 
-    PAR[PM_PAR_7] = (nStart > 0) ? 0.5*(indexGuessInv[nStart - 1] + indexGuessInv[nStart]) : indexGuessInv[nStart + 1];
+    PAR[PM_PAR_7] = Sp;
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);
     if (EXTRA_VERBOSE) fprintf (stderr, "%d >>> %d %f : %f - %f %f - %f %f %f - %f\n", source->id, model->nIter, model->chisqNorm, PAR[7], PAR[2], PAR[3], PAR[4], PAR[5], PAR[6], PAR[1]);
@@ -1036 +1051 @@
     float Smin = -0.5 * poly->coeff[1] / poly->coeff[2] / 100.0;
 
-    // constrain Smin to be in a valid range (1.0 - 0.1, corresponding to 0.5 (Gauss) to 5.0 (slightly peakier than Dev)
-    Smin = PS_MAX(PS_MIN(Smin, 1.0), 0.1);
+    // constrain Smin to be in a valid range: allow the fitted range to go a bit beyond the 3 trial points, but no further
+    float Smx = Sm - 0.25*(So - Sm);
+    float Spx = Sp + 0.25*(Sp - So);
+    Smin = PS_MAX(PS_MIN(Smin, Smx), Spx);
     PAR[PM_PAR_7] = Smin;
 
@@ -1249 +1266 @@
 
     // loop over Reff, keeping the ARatio and Theta constant
-    for (int j = -20; j < 21; j++) {
-
-        float dref = j * 0.02;
+    for (int j = -4; j <= 4; j++) {
+
+        float dref = j * 0.01;
 
         psEllipseAxes guessAxes;
@@ -1306 +1323 @@
             rMin = dref;
         }
-        fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
+        // fprintf (stderr, "%d | %f %f %f | %f %f %f\n", j, dref, Io, Chisq, rMin, iMin, xMin);
     }
 

branches/eam_branches/ipp-20130904/psphot/test/tap_psphot_galaxygrid.pro

@@ -325 +325 @@
   label -x sequence -y "M_out| - M_in|"
   resize 700 320
+
+  # check on magnitude
+  set dI = Iot_s - Iin_s
+  lim -n 4$1 n -1.0 1.0; clear; box; plot n dI
+  label -x sequence -y "I_out| - I_in|"
+  resize 700 320
 end
 
@@ -624 +630 @@
 
 macro cmf.load.reset
-  $fields = X Y M T R r MT
+  $fields = X Y M T R r MT I
   foreach field $fields
     foreach set in ot
@@ -636 +642 @@
 macro cmf.load.concat
   if ($0 != 4)
-    echo "USAGE: cmf.load.concat (dat) (cmf) (type)"
+    echo "USAGE: cmf.load.concat (dat) (cmf) (inType)"
     break
   end
@@ -643 +649 @@
   read Xin_all 1 Yin_all 2 Fin_all 3 Type 4 Min_all 5 RmajIn_all 7 RminIn_all 8 ThetaIn_all 9 IndexIn_all 10
 
+  $TYPE_S = 83
+  $TYPE_D = 68
+  $TYPE_E = 69
+  if ("$3" == "SERSIC")
+    $InType = $TYPE_S
+  end
+  if ("$3" == "DEV")
+    $InType = $TYPE_D
+    set IndexIn_all = 0.125 + zero(Xin_all)
+  end
+  if ("$3" == "EXP")
+    $InType = $TYPE_E
+    set IndexIn_all = 0.5 + zero(Xin_all)
+  end
+
+  # select only the galaxies 
   subset Xin = Xin_all if (Type == 1)
   subset Yin = Yin_all if (Type == 1)
@@ -657 +679 @@
   subset IndexIn = IndexIn_all if (Type == 1)
 
-  $TYPE_S = 83
-  $TYPE_D = 68
-  $TYPE_E = 69
-  if ("$3" == "SERSIC")
-    $InType = $TYPE_S
-  end
-  if ("$3" == "DEV")
-    $InType = $TYPE_D
-  end
-  if ("$3" == "EXP")
-    $InType = $TYPE_E
-  end
-
   data $2
 
@@ -675 +684 @@
   output -err /dev/null
   read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE EXT_PAR_07
+  $reread = not($STATUS)
   output -err stderr
-  $reread = not($STATUS)
   break -auto on
   if ($reread)
     read -fits Chip.xfit X_EXT Y_EXT EXT_INST_MAG EXT_WIDTH_MAJ EXT_WIDTH_MIN EXT_THETA MODEL_TYPE 
-    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*4 + (MODEL_TYPE:9 == 69)
+    set EXT_PAR_07 = (MODEL_TYPE:9 == 68)*0.125 + (MODEL_TYPE:9 == 69)*0.5
   end
 
@@ -709 +718 @@
   reindex rin_m = RminIn using index2
   
-  # if ("$3" == "SERSIC")
-    reindex Iot_m = EXT_PAR_07 using index1
-    reindex Iin_m = IndexIn using index2
-   $fields = X Y M T R r MT I
-  # else
-  # $fields = X Y M T R r MT
-  # end
+  reindex Iot_m = EXT_PAR_07 using index1
+  reindex Iin_m = IndexIn using index2
+  $fields = X Y M T R r MT I
   
   foreach field $fields
@@ -1286 +1291 @@
     foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
       mkexp.devexp.single tests.20131120/testrad.ps1v1.exp.$Nrun EXP $radius 1.0
-      fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4p EXP_CONV
+      fitexp tests.20131120/testrad.ps1v1.exp.$Nrun tests.20131120/testrad.ps1v1.exp.$Nrun.t4b EXP_CONV
       $Nrun ++
     end
@@ -1293 +1298 @@
   $Nrun = 0
   foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
-    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4p.cmf EXP
-    set dM_m = Mot_m - Min_m
-    vstat -q dM_m
-    $mag_off = $MEDIAN
-    set dR_m = Rot_m - Rin_m
-    vstat -q dR_m
-    $rad_off = $MEDIAN
-    echo $radius $mag_off $rad_off
+    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
+    echo -no-return $radius ""
+    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.t4b.cmf EXP
     $Nrun ++
   end
@@ -1319 +1319 @@
   $Nrun = 0
   foreach radius 1.0 1.5 2.0 2.5 3.0 4.0 6.0 8.0
-    cmf.load.concat tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP
-    set dM_m = Mot_m - Min_m
-    vstat -q dM_m
-    $mag_off = $MEDIAN
-    set dR_m = Rot_m - Rin_m
-    vstat -q dR_m
-    $rad_off = $MEDIAN
-    echo $radius $mag_off $rad_off
+    echo -no-return $radius ""
+    check.fit tests.20131120/testrad.ps1v1.exp.$Nrun.dat tests.20131120/testrad.ps1v1.exp.$Nrun.$1.cmf EXP
     $Nrun ++
   end
@@ -1388 +1382 @@
 end
 
-macro check.convolve.loop
-  if ($0 != 2)
-    echo "check.convolve.loop (intype)"
+# EXP : tests.20131120/test.nsig
+# DEV : tests.20131120/test.dev
+
+# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.nsig   EXP EXP_CONV
+# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.dev    DEV DEV_CONV
+# run.convolve.loop tests.20131120/test.nsig tests.20131120/test.serexp EXP SER_CONV
+# run.convolve.loop tests.20131120/test.dev  tests.20131120/test.serdev DEV SER_CONV
+
+macro run.convolve.loop
+  if ($0 != 5)
+    echo "run.convolve.loop (inName) (fitName) (inType) (fitType)"
     break
   end
 
   local radius
-
   $radius = 3.0
-
   cmf.load.reset
-
-  # foreach Cin 3 5 7 9 11
-  foreach Cin 11
-    sprintf nameIn tests.20131120/test.nsig.%02d $Cin
-    foreach Cot 3 5 7 9 11
-      sprintf nameOt tests.20131120/test.serdev.%02d.%02d $Cin $Cot
-      echo -no-return $Cin $Cot " "
-      check.fit $nameIn.dat $nameOt.cmf $1
-    end
-  end
-  grid.plots.devexp a
-end
-
-macro run.convolve.loop
-  local radius
-
-  $radius = 3.0
 
   # foreach Cin 3 5 7 9 11
   foreach Cin 11
     $CONVOLVE_NSIGMA = $Cin
-    sprintf nameIn tests.20131120/test.nsig.%02d $Cin
-    # mkexp.devexp.single $nameIn EXP $radius 1.0
+    sprintf nameIn %s.%02d $1 $Cin
+    mkexp.devexp.single $nameIn $3 $radius 1.0
     foreach Cot 3 5 7 9 11
       $NSIGMA_CONV = $Cot
-      sprintf nameOt tests.20131120/test.nsig.x.%02d.%02d $Cin $Cot
-      fitexp $nameIn $nameOt EXP_CONV
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      fitexp $nameIn $nameOt $4
     end
   end
@@ -1431 +1414 @@
   # foreach Cin 3 5 7 9 11
   foreach Cin 11
-    sprintf nameIn tests.20131120/test.nsig.%02d $Cin
+    sprintf nameIn %s.%02d $1 $Cin
     foreach Cot 3 5 7 9 11
-      sprintf nameOt tests.20131120/test.nsig.x.%02d.%02d $Cin $Cot
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
       check.fit $nameIn.dat $nameOt.cmf 
     end
@@ -1440 +1423 @@
 end
 
-macro run.convolve.loop.dev
+macro check.convolve.loop
+  if ($0 != 4)
+    echo "run.convolve.loop (inName) (fitName) (inType)"
+    break
+  end
+
   local radius
-
   $radius = 3.0
-
-  foreach Cin 11
-    $CONVOLVE_NSIGMA = $Cin
-    sprintf nameIn tests.20131120/test.dev.%02d $Cin
-    mkexp.devexp.single $nameIn DEV $radius 1.0
-    foreach Cot 3 5 7 9 11
-      $NSIGMA_CONV = $Cot
-      sprintf nameOt tests.20131120/test.dev.%02d.%02d $Cin $Cot
-      fitexp $nameIn $nameOt DEV_CONV
-    end
-  end
+  cmf.load.reset
 
   # foreach Cin 3 5 7 9 11
   foreach Cin 11
-    sprintf nameIn tests.20131120/test.dev.%02d $Cin
+    sprintf nameIn %s.%02d $1 $Cin
     foreach Cot 3 5 7 9 11
-      sprintf nameOt tests.20131120/test.dev.%02d.%02d $Cin $Cot
-      check.fit $nameIn.dat $nameOt.cmf DEV
+      sprintf nameOt %s.%02d.%02d $2 $Cin $Cot
+      echo -no-return $Cin $Cot " "
+      check.fit $nameIn.dat $nameOt.cmf $3
     end
   end
   grid.plots.devexp a
 end
-
-macro run.convolve.loop.serexp
-  local radius
-
-  $radius = 3.0
-
-  foreach Cin 11
-    $CONVOLVE_NSIGMA = $Cin
-    sprintf nameIn tests.20131120/test.nsig.%02d $Cin
-    foreach Cot 3 5 7 9 11
-      $NSIGMA_CONV = $Cot
-      sprintf nameOt tests.20131120/test.serexp.%02d.%02d $Cin $Cot
-      fitexp $nameIn $nameOt SER_CONV
-    end
-  end
-
-  # foreach Cin 3 5 7 9 11
-  foreach Cin 11
-    sprintf nameIn tests.20131120/test.nsig.%02d $Cin
-    foreach Cot 3 5 7 9 11
-      sprintf nameOt tests.20131120/test.serexp.%02d.%02d $Cin $Cot
-      check.fit $nameIn.dat $nameOt.cmf EXP
-    end
-  end
-  grid.plots.devexp a
-end
-
-
-macro run.convolve.loop.serdev
-  local radius
-
-  $radius = 3.0
-
-  cmf.load.reset
-
-  foreach Cin 11
-    $CONVOLVE_NSIGMA = $Cin
-    sprintf nameIn tests.20131120/test.dev.%02d $Cin
-    foreach Cot 3 5 7 9 11
-      $NSIGMA_CONV = $Cot
-      sprintf nameOt tests.20131120/test.serdev.%02d.%02d $Cin $Cot
-      fitexp $nameIn $nameOt SER_CONV
-    end
-  end
-
-  # foreach Cin 3 5 7 9 11
-  foreach Cin 11
-    sprintf nameIn tests.20131120/test.dev.%02d $Cin
-    foreach Cot 3 5 7 9 11
-      sprintf nameOt tests.20131120/test.serdev.%02d.%02d $Cin $Cot
-      check.fit $nameIn.dat $nameOt.cmf DEV
-    end
-  end
-  grid.plots.devexp a
-end
-