Changeset 40645 for trunk

Timestamp:

Mar 19, 2019, 6:30:13 AM (7 years ago)

Author:

eugene

Message:

working on fitobj

File:

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

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

@@ -1 +1 @@
 # include "data.h"
 # include "deimos.h"
+
+int deimos_fitobj (int argc, char **argv) {
+  OHANA_UNUSED_PARAM(argc);
+  OHANA_UNUSED_PARAM(argv);
+  return FALSE;
+}
+
+# if (0)
+
+/*
+  this is starting to work OK.  some improvements to make
+  * save all fitted obj,sky,bck vectors
+  * save chisq for each fit set
+  * use Gaussdev to sample
+  * do not scale down range (or user-set scale-down)
+ */
 
 // internal functions to fitobj
@@ -144 +160 @@
   }
 
-  // Save the best 3 results for each of obj, sky, bck:
-  Vector *obj1, *obj2, *obj3;
-  if ((obj1 = SelectVector ("obj1", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (obj1, OPIHI_FLT, Ny);
-  if ((obj2 = SelectVector ("obj2", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (obj2, OPIHI_FLT, Ny);
-  if ((obj3 = SelectVector ("obj3", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (obj3, OPIHI_FLT, Ny);
-  Vector *sky1, *sky2, *sky3;
-  if ((sky1 = SelectVector ("sky1", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (sky1, OPIHI_FLT, Ny);
-  if ((sky2 = SelectVector ("sky2", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (sky2, OPIHI_FLT, Ny);
-  if ((sky3 = SelectVector ("sky3", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (sky3, OPIHI_FLT, Ny);
-  Vector *bck1, *bck2, *bck3;
-  if ((bck1 = SelectVector ("bck1", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bck1, OPIHI_FLT, Ny);
-  if ((bck2 = SelectVector ("bck2", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bck2, OPIHI_FLT, Ny);
-  if ((bck3 = SelectVector ("bck3", ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (bck3, OPIHI_FLT, Ny);
-
-  for (int iy = 0; iy < Ny; iy++) {
-    obj1->elements.Flt[iy] = 0.0; sky1->elements.Flt[iy] = 0.0; bck1->elements.Flt[iy] = 0.0;
-    obj2->elements.Flt[iy] = 0.0; sky2->elements.Flt[iy] = 0.0; bck2->elements.Flt[iy] = 0.0;
-    obj3->elements.Flt[iy] = 0.0; sky3->elements.Flt[iy] = 0.0; bck3->elements.Flt[iy] = 0.0;
-  }
-
   // for the cross-dispersion reference pixle, use the user value or set to Nx/2 if < 0
   deimos_set_cross_ref (atoi(argv[2]), Nx); 
@@ -175 +171 @@
   int NRESULT = Niter + Ntry*Niter;
   ALLOCATE_PTR (result, DeimosResult, NRESULT);
+
+  // Save all result vectors for each of obj, sky, bck:
+  ALLOCATE_PTR (objOut, Vector, NRESULT);
+  ALLOCATE_PTR (skyOut, Vector, NRESULT);
+  ALLOCATE_PTR (bckOut, Vector, NRESULT);
+  for (int i = 0; i < NRESULT; i++) {
+    char name[64];
+    snprintf (name, 64, "obj_%04d", i); if ((objOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (&objOut[i], OPIHI_FLT, Ny);
+    snprintf (name, 64, "sky_%04d", i); if ((skyOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (&skyOut[i], OPIHI_FLT, Ny);
+    snprintf (name, 64, "bck_%04d", i); if ((bckOut[i] = SelectVector (name, ANYVECTOR, TRUE)) == NULL) { return (FALSE); } ResetVector (&bckOut[i], OPIHI_FLT, Ny);
+    for (int iy = 0; iy < Ny; iy++) {
+      objOut[i].elements.Flt[iy] = 0.0;
+      skyOut[i].elements.Flt[iy] = 0.0;
+      bckOut[i].elements.Flt[iy] = 0.0;
+    }
+  }
+  Vector *chisqVect;
+  if ((chisqVect = SelectVector ("chisqOut", ANYVECTOR, TRUE)) == NULL) { return FALSE; }
+  ResetVector (chisqVect, OPIHI_FLT, NRESULT);
+  for (int i = 0; i < NRESULT; i++) { chisqVect->elements.Flt[i] = 0.0; }
 
   float noiseVar = SQ(noise);
@@ -302 +318 @@
     // int *index = sort_result (result, Nresult);
     // copy result[index[0]].obj,sky,bck to obj,sky,bck
-    int Nv = 0;
-    for (int iy = 0; iy < Nrow; iy++) {
-      Nv = IDX[0]; obj1->elements.Flt[iy + row] = result[Nv].obj[iy]; sky1->elements.Flt[iy + row] = result[Nv].sky[iy]; bck1->elements.Flt[iy + row] = result[Nv].bck[iy];
-      Nv = IDX[1]; obj2->elements.Flt[iy + row] = result[Nv].obj[iy]; sky2->elements.Flt[iy + row] = result[Nv].sky[iy]; bck2->elements.Flt[iy + row] = result[Nv].bck[iy];
-      Nv = IDX[2]; obj3->elements.Flt[iy + row] = result[Nv].obj[iy]; sky3->elements.Flt[iy + row] = result[Nv].sky[iy]; bck3->elements.Flt[iy + row] = result[Nv].bck[iy];
+    for (int i = 0; i < Nresult; i++) {
+      for (int iy = 0; iy < Nrow; iy++) {
+        objOut[i].elements.Flt[iy + row] = result[IDX[i]].obj[iy];
+        skyOut[i].elements.Flt[iy + row] = result[IDX[i]].sky[iy];
+        bckOut[i].elements.Flt[iy + row] = result[IDX[i]].bck[iy];
+      }
+      chisqVect->elements.Flt[i] += chisq;
     }
   }
@@ -402 +420 @@
   return bin;
 }
+# endif