Changeset 40836

Timestamp:

Jul 19, 2019, 5:31:16 PM (7 years ago)

Author:

eugene

Message:

avoid NANs in peak detection; do not allow splines with < 3 pts; add vgauss -apply option

Location:

branches/eam_branches/ohana.20190329/src/opihi/cmd.data

Files:

• peak.c (modified) (2 diffs)
• spline_commands.c (modified) (6 diffs)
• vgauss.c (modified) (5 diffs)

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/peak.c

@@ -39 +39 @@
       if (*X < start) continue;
       if (*X > end) continue;
+      if (!isfinite(*Y)) continue;
       if (*Y < ymax) continue;
       xmax = *X;
@@ -63 +64 @@
       if (*X < start) continue;
       if (*X > end) continue;
+      if (!isfinite(*Y)) continue;
       if (*Y < ymax) continue;
       xmax = *X;

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/spline_commands.c

@@ -28 +28 @@
   if (xvec->Nelements != yvec->Nelements) {
     gprint (GP_ERR, "x and y vectors must be the same length\n");
+    return FALSE;
+  }
+
+  if (xvec->Nelements < 3) {
+    gprint (GP_ERR, "cannot make a spline with fewer than 3 knots\n");
     return FALSE;
   }
@@ -132 +137 @@
 
   int N, status;
-  Spline *spline;
 
   int QUIET = FALSE;
@@ -145 +149 @@
   }
 
-  spline = FindSpline (argv[1]);
+  Spline *spline = FindSpline (argv[1]);
   if (spline == NULL) {
     if (QUIET) return TRUE;
@@ -159 +163 @@
 int spline_rename (int argc, char **argv) {
 
-  Spline *spline;
-
   if (argc != 3) {
     gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n");
@@ -166 +168 @@
   }
 
-  spline = FindSpline (argv[1]);
+  Spline *spline = FindSpline (argv[1]);
   if (spline == NULL) {
     gprint (GP_ERR, "spline %s not found\n", argv[1]);
     return FALSE;
+  }
+
+  Spline *newname = FindSpline (argv[2]);
+  if (newname != NULL) {
+    DeleteSpline (newname);
   }
 
@@ -176 +183 @@
   return TRUE;
 }
-
-/* 
-int spline_listspline (int argc, char **argv) {
-
-  Spline *spline;
-
-  if (argc != 2) {
-    gprint (GP_ERR, "USAGE: spline listspline (spline)\n");
-    return FALSE;
-  }
-
-  spline = FindSpline (argv[1]);
-  if (spline == NULL) {
-    gprint (GP_ERR, "spline %s not found\n", argv[1]);
-    return FALSE;
-  }
-
-  ListPages (spline);
-  return TRUE;
-}
-*/
-

branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vgauss.c

@@ -3 +3 @@
 /* local private functions */
 opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *);
+int vgauss_apply (int argc, char **argv);
 
 # define GET_VAR(V,A) \
@@ -20 +21 @@
   Vector *xvec, *yvec, *svec, *ovec;
   char *c, name[16];
+
+  if ((N = get_argument (argc, argv, "-apply"))) {
+    remove_argument (N, &argc, argv);
+    int status = vgauss_apply (argc, argv);
+    return status;
+  }
 
   Quiet = FALSE;
@@ -104 +111 @@
   }  
   if (!Quiet) gprint (GP_ERR, "%d iterations\n", i); 
+  set_variable ("Chisq", chisq);
+  set_variable ("ChisqNu", chisq / (float) (Npts - Npar));
 
   for (i = 0; i < Npts; i++) {
@@ -125 +134 @@
 }
 
+int vgauss_apply (int argc, char **argv) {
+
+  opihi_flt par[4];
+  Vector *xvec, *ovec;
+  char *c;
+
+  // NOTE: -apply has already been stripped by vgauss call
+  if (argc != 3) {
+    gprint (GP_ERR, "USAGE: vgauss -apply <x> <out>\n");
+    return (FALSE);
+  }
+  
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((ovec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  CastVector (xvec, OPIHI_FLT);
+
+  GET_VAR (par[0], "C0");
+  GET_VAR (par[1], "C1");
+  GET_VAR (par[2], "C2");
+  GET_VAR (par[3], "C3");
+  int Npar = 4;
+
+  int Npts = xvec[0].Nelements;
+  ResetVector (ovec, OPIHI_FLT, Npts);
+
+  for (int i = 0; i < Npts; i++) {
+    ovec[0].elements.Flt[i] = fgaussOD (xvec[0].elements.Flt[i], par, Npar, NULL);
+  }
+  ovec[0].Nelements = Npts;
+
+  return TRUE;
+}
+
 /* pars: x_o, sigma, I, back */
 opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
@@ -135 +178 @@
   f = par[2]*r + par[3];
 
-  dpar[0] = par[2]*r*z/par[1];
-  dpar[1] = par[2]*r*z*z/par[1];
-  dpar[2] = r;
-  dpar[3] = 1;
+  if (dpar) {
+    dpar[0] = par[2]*r*z/par[1];
+    dpar[1] = par[2]*r*z*z/par[1];
+    dpar[2] = r;
+    dpar[3] = 1;
+  }
   
   return (f);