Changeset 40643

Timestamp:

Mar 18, 2019, 5:15:30 PM (7 years ago)

Author:

eugene

Message:

cleanup opihi sort functions; add deimos fitobj function

File:

• trunk/Ohana/src/opihi/mana/deimos_fitobj.c (modified) (11 diffs)

trunk/Ohana/src/opihi/mana/deimos_fitobj.c

@@ -6 +6 @@
 int       *deimos_sort_result (DeimosResult *result, int Nresult);
 int        deimos_get_random_sample (int Nresult);
-opihi_flt *deimos_make_test (opihi_flt *guess, float range, int row, int Nrow, int fullInput);
+opihi_flt *deimos_make_test (opihi_flt *guess, float *sigma, int row, int Nrow, int fullInput);
 
 int deimos_fitobj (int argc, char **argv) {
@@ -36 +36 @@
 
   float stilt = 0.0; // angle of the slit
+
+  // a crude noise model valid for both the guess vectors and the
+  // data buffer:
+  float gain = 1.0;
+  if ((N = get_argument (argc, argv, "-gain"))) {
+    remove_argument (N, &argc, argv);
+    gain = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+  float noise = 0.0;
+  if ((N = get_argument (argc, argv, "-noise"))) {
+    remove_argument (N, &argc, argv);
+    noise = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
 
   int Niter = 25;
@@ -161 +176 @@
   ALLOCATE_PTR (result, DeimosResult, NRESULT);
 
+  float noiseVar = SQ(noise);
+  ALLOCATE_PTR (objNoise, float, Nrow);
+  ALLOCATE_PTR (skyNoise, float, Nrow);
+  ALLOCATE_PTR (bckNoise, float, Nrow);
+
   // work on segments of Nrow at a time
   // XXX worry about last segment (< Nrow rows)
@@ -174 +194 @@
     // likely have large errors in how the guess follows the truth. 
 
+    /** old concept:
     float objRange = 0, skyRange = 0, bckRange = 0;
     for (int iy = 0; iy < Nrow; iy++) {
@@ -180 +201 @@
       bckRange = MAX(bckVal[iy + row], bckRange);
     }
+    **/
+
+    // new concept: start with 2-sigma range, then ramp down below
+    for (int iy = 0; iy < Nrow; iy++) {
+      objNoise[iy] = 2*hypot(sqrt(fabs(objVal[iy + row]) / gain), noiseVar);
+      skyNoise[iy] = 2*hypot(sqrt(fabs(skyVal[iy + row]) / gain), noiseVar);
+      bckNoise[iy] = 2*hypot(sqrt(fabs(bckVal[iy + row]) / gain), noiseVar);
+    }
 
     // first step, generate Niter samples (chisq + parameters) 
@@ -187 +216 @@
 
       // generate a set of test values for obj, sky, backgnd based on the current guess
-      opihi_flt *objtest = deimos_make_test (objVal, objRange, row, Nrow, TRUE);
-      opihi_flt *skytest = deimos_make_test (skyVal, skyRange, row, Nrow, TRUE);
-      opihi_flt *bcktest = deimos_make_test (bckVal, bckRange, row, Nrow, TRUE);
+      opihi_flt *objtest = deimos_make_test (objVal, objNoise, row, Nrow, TRUE);
+      opihi_flt *skytest = deimos_make_test (skyVal, skyNoise, row, Nrow, TRUE);
+      opihi_flt *bcktest = deimos_make_test (bckVal, bckNoise, row, Nrow, TRUE);
       
       for (int itmp = 0; FALSE && (itmp < Nrow); itmp++) {
@@ -223 +252 @@
 
       // XXX not sure this scaling factor, or even concept is good
-      objRange *= 0.75;
-      skyRange *= 0.75;
-      bckRange *= 0.75;
+      // objRange *= 0.75;
+      // skyRange *= 0.75;
+      // bckRange *= 0.75;
+      for (int iy = 0; iy < Nrow; iy++) {
+        objNoise[iy] = 0.75*objNoise[iy];
+        skyNoise[iy] = 0.75*skyNoise[iy];
+        bckNoise[iy] = 0.75*bckNoise[iy];
+      }
 
       // run Niter iterations, selecting guess values from the result set
@@ -234 +268 @@
 
         // generate a set of test values for obj, sky, backgnd based on the current guess
-        opihi_flt *objtest = deimos_make_test (result[entry].obj, objRange, row, Nrow, FALSE);
-        opihi_flt *skytest = deimos_make_test (result[entry].sky, skyRange, row, Nrow, FALSE);
-        opihi_flt *bcktest = deimos_make_test (result[entry].bck, bckRange, row, Nrow, FALSE);
+        opihi_flt *objtest = deimos_make_test (result[entry].obj, objNoise, row, Nrow, FALSE);
+        opihi_flt *skytest = deimos_make_test (result[entry].sky, skyNoise, row, Nrow, FALSE);
+        opihi_flt *bcktest = deimos_make_test (result[entry].bck, bckNoise, row, Nrow, FALSE);
       
         for (int itmp = 0; FALSE && (itmp < Nrow); itmp++) {
@@ -285 +319 @@
 /****************** fitobj Support Functions *******************/
 
-opihi_flt *deimos_make_test (opihi_flt *guess, float range, int row, int Nrow, int fullInput) {
+opihi_flt *deimos_make_test (opihi_flt *guess, float *sigma, int row, int Nrow, int fullInput) {
 
   // we have a vector, guess, with some length >= row + Nrow.
@@ -293 +327 @@
   // if fullInput is TRUE,  the guess vector runs from 0 to Ny (full wavelength range)
   // if fullInput is FALSE, the guess vector runs from row to row + Nrow (subset range)
+  // sigma is always only Nrow long
 
   ALLOCATE_PTR (value, opihi_flt, Nrow);
@@ -301 +336 @@
   int iv = 0;
   for (int i = start_pix; i < stop_pix; i++, iv++) {
-    value[iv] = range*(drand48() - 0.5) + guess[i];
+    value[iv] = sigma[iv]*(drand48() - 0.5) + guess[i];
   }
   return value;