Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile	(revision 40751)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/Makefile	(revision 40752)
@@ -90,4 +90,5 @@
 $(SRC)/lookup.$(ARCH).o	\
 $(SRC)/matrix.$(ARCH).o	\
+$(SRC)/match1d.$(ARCH).o	\
 $(SRC)/match2d.$(ARCH).o	\
 $(SRC)/mkrgb.$(ARCH).o	\
@@ -170,4 +171,5 @@
 $(SRC)/vgauss.$(ARCH).o            \
 $(SRC)/vlorentz.$(ARCH).o          \
+$(SRC)/vsigmoid.$(ARCH).o          \
 $(SRC)/vellipse.$(ARCH).o          \
 $(SRC)/vmaxwell.$(ARCH).o          \
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c	(revision 40751)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/init.c	(revision 40752)
@@ -80,4 +80,5 @@
 int lookup           PROTO((int, char **));
 int matrix           PROTO((int, char **));
+int match1d          PROTO((int, char **));
 int match2d          PROTO((int, char **));
 int mkrgb            PROTO((int, char **));
@@ -157,4 +158,5 @@
 int vgauss           PROTO((int, char **));
 int vlorentz         PROTO((int, char **));
+int vsigmoid         PROTO((int, char **));
 int vellipse         PROTO((int, char **));
 int vmaxwell         PROTO((int, char **));
@@ -279,4 +281,5 @@
   {1, "medimage",     medimage_command, "median image manipulation"},
   {1, "matrix",       matrix,           "matrix math operations"},
+  {1, "match1d",      match1d,          "match 2 vectors and return matched indexes"},
   {1, "match2d",      match2d,          "match 2 pairs of X,Y vectors and return matched indexes"},
   {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
@@ -353,4 +356,5 @@
   {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
   {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
+  {1, "vsigmoid",     vsigmoid,         "fit a Sigmoid to a vector"},
   {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
   {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/match1d.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/match1d.c	(revision 40752)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/match1d.c	(revision 40752)
@@ -0,0 +1,248 @@
+# include "data.h"
+
+int find_matches1d (Vector *X1, Vector *X2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch);
+int find_matches1d_closest (Vector *X1, Vector *X2, double Radius, Vector *index);
+
+// match1d (X1) (X2) (Radius) [-index1 (index1)] [-index2 (index2)] [-nomatch1 nomatch1] [-nomatch2 nomatch2]
+int match1d (int argc, char **argv) {
+  
+  int N, CLOSEST;
+  double Radius;
+  char *endptr;
+  Vector *X1vec, *X2vec;
+  Vector *index1, *index2;
+
+  if ((N = get_argument (argc, argv, "-h"))) goto usage;
+  if ((N = get_argument (argc, argv, "--help"))) goto usage;
+
+  CLOSEST = FALSE;
+  if ((N = get_argument (argc, argv, "-closest"))) {
+    remove_argument (N, &argc, argv);
+    CLOSEST = TRUE;
+  }
+
+  if ((N = get_argument (argc, argv, "-index1"))) {
+    remove_argument (N, &argc, argv);
+    if ((index1 = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  } else {
+    if ((index1 = SelectVector ("index1", ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+  }
+
+  if ((N = get_argument (argc, argv, "-index2"))) {
+    remove_argument (N, &argc, argv);
+    if ((index2 = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  } else {
+    if ((index2 = SelectVector ("index2", ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+  }
+
+  Vector *radiusMatch = NULL;
+  if ((N = get_argument (argc, argv, "-radius"))) {
+    if (CLOSEST) {
+      gprint (GP_ERR, "error: -radius and -closest are currently incompatible\n");
+      return (FALSE);
+    }
+    remove_argument (N, &argc, argv);
+    if ((radiusMatch = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);    
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 4) {
+    gprint (GP_ERR, "USAGE: match1d X1 X2 Radius [-index1 (index1)] [-index2 (index2)] [-closest]\n");
+    gprint (GP_ERR, "  use -h or --help for more detail\n");
+    return (FALSE);
+  }
+
+  /*
+    we have two modes of operation:  
+
+    without -closest, we are finding all matched pairs within the match radius.  in this
+    case, the two index vectors have the same length, one entry per matched pair.
+    x1[index1],y1[index1] matches to x2[index2],y2[index2].
+
+    with -closest selected, we are finding the closest element of set 1 to each of set 2
+    and vice versa.  in this case, index1 is always the same length as x1,y1, while index2
+    is the same lengths as x2,y2.  x2[index1],y2[index1] matches x1,y1 while
+    x1[index2],y1[index2] matches x2,y2
+
+   */
+
+  if ((X1vec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
+  if ((X2vec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);    
+
+  REQUIRE_VECTOR_FLT (X1vec, FALSE); 
+  REQUIRE_VECTOR_FLT (X2vec, FALSE); 
+
+  Radius = strtod (argv[3], &endptr);
+  if (*endptr) {
+    gprint (GP_ERR, "Radius must be numerical (%s)\n", argv[3]);
+    return (FALSE);
+  }
+
+  if (CLOSEST) {
+    find_matches1d_closest (X1vec, X2vec, Radius, index1);
+    find_matches1d_closest (X2vec, X1vec, Radius, index2);
+  } else {
+    find_matches1d (X1vec, X2vec, Radius, index1, index2, radiusMatch);
+  }
+  return (TRUE);
+
+usage:
+  gprint (GP_ERR, "we have two modes of operation:\n\n");
+
+  gprint (GP_ERR, "without -closest, we are finding all matched pairs within the match radius.  in this\n");
+  gprint (GP_ERR, "case, the two index vectors have the same length, one entry per matched pair.\n");
+  gprint (GP_ERR, "x1[index1] matches to x2[index2].\n\n");
+
+  gprint (GP_ERR, "with -closest selected, we are finding the closest element of set 1 to each of set 2\n");
+  gprint (GP_ERR, "and vice versa.  in this case, index1 is always the same length as x1, while index2\n");
+  gprint (GP_ERR, "is the same lengths as x2.  x2[index1] matches x1 while\n");
+  gprint (GP_ERR, "x1[index2] matches x2\n\n");
+
+  gprint (GP_ERR, "if -index1 or -index2 is not supplied, the vectors are created with names index1 or index2\n");
+  gprint (GP_ERR, "use 'reindex' to generate new vectors based on these index vectors\n");
+
+  gprint (GP_ERR, "if -radius (vector) is supplied, the vector will be filled with the distance between the matched pairs\n");
+  gprint (GP_ERR, "  not valid with -closest\n");
+  return FALSE;
+}
+
+// Radius is a 1-D separation
+int find_matches1d (Vector *X1, Vector *X2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch) {
+  
+  off_t i, j, first_j, I, J, *N1, *N2, Nmatch, NMATCH, DMATCH;
+  double dX, dR;
+
+  NMATCH = MAX(MAX(0.01*X1->Nelements, 0.01*X2->Nelements), 100);
+  DMATCH = NMATCH;
+
+  ResetVector (index1, OPIHI_INT, NMATCH);
+  ResetVector (index2, OPIHI_INT, NMATCH);
+  if (radiusMatch) ResetVector (radiusMatch, OPIHI_FLT, NMATCH);
+
+  ALLOCATE (N1, off_t, X1->Nelements);
+  ALLOCATE (N2, off_t, X2->Nelements);
+
+  for (i = 0; i < X1->Nelements; i++) { N1[i] = i; }
+  for (i = 0; i < X2->Nelements; i++) { N2[i] = i; }
+
+  dsort_indexonly (X1->elements.Flt, N1, X1->Nelements);
+  dsort_indexonly (X2->elements.Flt, N2, X2->Nelements);
+
+  Nmatch = 0;
+  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
+    I = N1[i];
+    J = N2[j];
+
+    if (!isfinite(X1->elements.Flt[I])) { i++; continue; }
+    if (!isfinite(X2->elements.Flt[J])) { j++; continue; }
+
+    dX = X1->elements.Flt[I] - X2->elements.Flt[J];
+
+    if (dX <= -1.02*Radius) { i++; continue; }
+    if (dX >= +1.02*Radius) { j++; continue; }
+
+    // look for all matches of list2() to list1(i)
+    first_j = j;
+    for (j = first_j; (dX > -1.02*Radius) && (j < X2->Nelements); j++) {
+      J = N2[j];
+      dX = X1->elements.Flt[I] - X2->elements.Flt[J];
+      dR = fabs(dX);
+      if (dR < Radius) {
+	index1->elements.Int[Nmatch] = I;
+	index2->elements.Int[Nmatch] = J;
+	if (radiusMatch) radiusMatch->elements.Flt[Nmatch] = dR;
+
+	// XXX track matches 1 and 2 with internal vector, save new nomatch index vectors
+	// after this loop
+
+	Nmatch ++;
+	if (Nmatch >= NMATCH) {
+	  NMATCH += DMATCH;
+	  REALLOCATE (index1->elements.Int, opihi_int, NMATCH);
+	  REALLOCATE (index2->elements.Int, opihi_int, NMATCH);
+	  if (radiusMatch) { REALLOCATE (radiusMatch->elements.Flt, opihi_flt, NMATCH); }
+	}
+      }
+    }
+    j = first_j;
+    i++;
+  }
+  index1->Nelements = Nmatch;
+  index2->Nelements = Nmatch;
+  if (radiusMatch) radiusMatch->Nelements = Nmatch;
+
+  free (N1);
+  free (N2);
+
+  return (TRUE);
+}
+
+// find the elements of X2,Y2 which are closest to each element of X1,Y1 (-1 if no match)
+int find_matches1d_closest (Vector *X1, Vector *X2, double Radius, Vector *index) {
+  
+  off_t i, j, Jmin, Ji, I, J, *N1, *N2, NMATCH;
+  double dX, dR, Rmin;
+
+  NMATCH = X1->Nelements;
+  ResetVector (index, OPIHI_INT, NMATCH);
+
+  for (i = 0; i < index->Nelements; i++) { index->elements.Int[i] = -1; }
+
+  ALLOCATE (N1, off_t, X1->Nelements);
+  ALLOCATE (N2, off_t, X2->Nelements);
+
+  for (i = 0; i < X1->Nelements; i++) { N1[i] = i; }
+  for (i = 0; i < X2->Nelements; i++) { N2[i] = i; }
+
+  dsort_indexonly (X1->elements.Flt, N1, X1->Nelements);
+  dsort_indexonly (X2->elements.Flt, N2, X2->Nelements);
+
+  // find the closest entry in list 2 to the current entry in list 1:
+  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
+    I = N1[i];
+    J = N2[j];
+
+    dX = X1->elements.Flt[I] - X2->elements.Flt[J];
+
+    if (dX <= -1.02*Radius) { 
+      // no match in list 2 to this entry
+      index->elements.Int[I] = -1;
+      i++; 
+      continue; 
+    }
+    if (dX >= +1.02*Radius) { j++; continue; }
+
+    // look for closest matches of list2() to list1(i)
+    Jmin = -1;
+    for (Ji = j; (dX > -1.02*Radius) && (Ji < X2->Nelements); Ji++) {
+      J = N2[Ji];
+      dX = X1->elements.Flt[I] - X2->elements.Flt[J];
+      dR = fabs(dX);
+      if (dR > Radius) continue;
+      if (dR < Rmin) {
+	Rmin = dR;
+	Jmin  = J;
+      }
+    }
+
+    // no match in list 2 to this entry
+    if (Jmin == -1) {
+      index->elements.Int[I] = -1;
+      i++;
+      continue;
+    }
+    index->elements.Int[I] = Jmin;
+    i++;
+  }
+  index->Nelements = NMATCH;
+
+  free (N1);
+  free (N2);
+
+  return (TRUE);
+}
+
+
+
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/vsigmoid.sh
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/vsigmoid.sh	(revision 40752)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/test/vsigmoid.sh	(revision 40752)
@@ -0,0 +1,218 @@
+
+list tests
+ test1
+ test2
+ test3
+end
+
+macro test1
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 0 + 5 / (1 + exp((x - 3)/4))
+
+ $C0 = 1
+ $C1 = 2
+ $C2 = 10
+ $C3 = 1
+ set dy = sqrt(y)
+
+ vsigmoid -q x y dy yfit
+ lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 5.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.01)
+   $PASS = 0
+ end
+end
+
+# noise of 0.1
+macro test2
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 50 + 1000 / (1 + exp((x - 3)/2))
+ set dy = (rnd(y) - 0.5)/0.5
+ set y = y + dy
+
+ $C0 = 2
+ $C1 = 3
+ $C2 = 900
+ $C3 = 1
+
+ vsigmoid -q x y dy yfit
+ lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 1000.0) > 0.1)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.1)
+   $PASS = 0
+ end
+end
+
+# poisson-distributed noise
+macro test3
+ $PASS = 1
+ break -auto off
+
+ $C0o = 3
+ $C1o = 2
+ $C2o = 20
+ $C3o = -2
+
+ create x -10 10 0.1
+ set y = $C3o + $C2o / (1 + exp((x - $C0o)/$C1o))
+
+ gaussdev dY y[] 0.0 1.0
+ set dy = dY * sqrt(abs(y))
+ set y = y + dy
+
+ $C0 = $C0o + 2
+ $C1 = $C1o + 2
+ $C2 = $C2o + 50
+ $C3 = $C3o + 10
+
+ vsigmoid x y dy yfit
+ lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ $dS = 3.0 / sqrt(y[])
+
+ if (abs($C0 - $C0o) > $dS)
+   $PASS = 0
+ end
+ if (abs($C1 - $C1o)/$C1o > $dS)
+   $PASS = 0
+ end
+ if (abs($C2 - $C2o)/$C2o > $dS)
+   $PASS = 0
+ end
+ if (abs($C3 - $C3o) > $dS)
+   $PASS = 0
+ end
+end
+
+macro test1rev
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 5 / (1 + exp(-(x - 3)/4))
+
+ $C0 = 1
+ $C1 = 2
+ $C2 = 10
+ $C3 = 1
+ set dy = sqrt(abs(y))
+
+ vsigmoid -r -q x y dy yfit
+ lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 5.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.01)
+   $PASS = 0
+ end
+end
+
+# noise of 0.1
+macro test2rev
+ $PASS = 1
+ break -auto off
+
+ create x -10 10 0.1
+ set y = 50 + 1000 / (1 + exp(-(x - 3)/2))
+ set dy = (rnd(y) - 0.5)/0.5
+ set y = y + dy
+
+ $C0 = 2
+ $C1 = 3
+ $C2 = 900
+ $C3 = 1
+
+ vsigmoid -r -q x y dy yfit
+ lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ if (abs($C0 - 0.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C1 - 3.0) > 0.01)
+   $PASS = 0
+ end
+ if (abs($C2 - 1000.0) > 0.1)
+   $PASS = 0
+ end
+ if (abs($C3 - 0.0) > 0.1)
+   $PASS = 0
+ end
+end
+
+# poisson-distributed noise
+macro test3rev
+ $PASS = 1
+ break -auto off
+
+ $C0o = 3
+ $C1o = 2
+ $C2o = 20
+ $C3o = -2
+
+ create x -10 10 0.1
+ set y = $C3o + $C2o / (1 + exp(-(x - $C0o)/$C1o))
+
+ gaussdev dY y[] 0.0 1.0
+ set dy = dY * sqrt(abs(y))
+ set y = y + dy
+
+ $C0 = $C0o + 2
+ $C1 = $C1o + 2
+ $C2 = $C2o + 50
+ $C3 = $C3o + 10
+
+ vsigmoid -r -q x y dy yfit
+ lim x y; clear; box; plot -x 1 -c black x y; plot -c red x yfit
+ # cursor
+
+ $dS = 3.0 / sqrt(y[])
+
+ if (abs($C0 - $C0o) > $dS)
+   $PASS = 0
+ end
+ if (abs($C1 - $C1o)/$C1o > $dS)
+   $PASS = 0
+ end
+ if (abs($C2 - $C2o)/$C2o > $dS)
+   $PASS = 0
+ end
+ if (abs($C3 - $C3o) > $dS)
+   $PASS = 0
+ end
+end
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/threshold.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/threshold.c	(revision 40751)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/threshold.c	(revision 40752)
@@ -25,4 +25,14 @@
     remove_argument (N, &argc, argv);
     BinMax = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  float ValMin = NAN;
+  float ValMax = NAN;
+  if ((N = get_argument (argc, argv, "-vrange"))) {
+    remove_argument (N, &argc, argv);
+    ValMin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ValMax = atof (argv[N]);
     remove_argument (N, &argc, argv);
   }
@@ -62,4 +72,16 @@
     }
     if (BinMax < 0) BinMax += vecx[0].Nelements;
+  }
+
+  // ValMin, ValMax not compatible with BinMin, BinMax
+  if (isfinite(ValMin)) {
+    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMin);
+    BinMin = (vecx[0].elements.Flt[binTest+1] == ValMin) ? binTest : binTest + 1;
+    // XXX dangerous: check that this does not run off the end, etc
+  }
+  if (isfinite(ValMax)) {
+    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMax);
+    BinMax = (vecx[0].elements.Flt[binTest+1] == ValMax) ? binTest : binTest + 1;
+    // XXX dangerous: check that this does not run off the end, etc
   }
 
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vsigmoid.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vsigmoid.c	(revision 40752)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vsigmoid.c	(revision 40752)
@@ -0,0 +1,192 @@
+# include "data.h"
+
+typedef opihi_flt FitFunc (opihi_flt, opihi_flt *, int, opihi_flt *);
+
+/* local private functions */
+opihi_flt fsigmoidOD (opihi_flt, opihi_flt *, int, opihi_flt *);
+opihi_flt rsigmoidOD (opihi_flt, opihi_flt *, int, opihi_flt *);
+
+# define GET_VAR(V,A) {					\
+    char *c = get_variable (A);				\
+    if (c == NULL) {					\
+      gprint (GP_ERR, "missing fit parameter A\n");	\
+      return (FALSE);					\
+    }							\
+    V = atof (c);					\
+    free (c); }
+
+// XXX this fitting function works but fails if the input vectors have NANs
+int vsigmoid (int argc, char **argv) {
+
+  int N;
+  Vector *xvec, *yvec, *svec, *ovec;
+
+  FitFunc *FUNC = fsigmoidOD;
+  if ((N = get_argument (argc, argv, "-r"))) {
+    FUNC = rsigmoidOD;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-rev"))) {
+    FUNC = rsigmoidOD;
+    remove_argument (N, &argc, argv);
+  }
+
+  int Quiet = FALSE;
+  if ((N = get_argument (argc, argv, "-q"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+  if ((N = get_argument (argc, argv, "-quiet"))) {
+    Quiet = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  if (argc != 5) {
+    gprint (GP_ERR, "USAGE: vsigmoid [-q] [-r] [-rev] <x> <y> <dy> (out)\n");
+    gprint (GP_ERR, "  <dy> may be the words 'con[stant]' or 'poi[sson]', in which case the error is constructed'\n");
+    gprint (GP_ERR, " uses guesses: C0 (center), C1 (sigma), C2 (top), C3 (bottom)\n");
+    gprint (GP_ERR, " forward sigmoid : bottom + top / (1 + exp(-z)), z = (x - center) / sigma\n");
+    gprint (GP_ERR, " -r / -rev : reverse sigmoid : bottom + top / (1 + exp(z))\n");
+    return (FALSE);
+  }
+  
+  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((yvec = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  if ((ovec = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
+
+  CastVector (xvec, OPIHI_FLT);
+  CastVector (yvec, OPIHI_FLT);
+
+  int Nsvec = strlen(argv[3]);
+
+  if ((Nsvec >= 3) && (!strncasecmp ("poisson", argv[3], Nsvec) || !strncasecmp ("constant", argv[3], Nsvec))) {
+    svec = SelectVector ("vsigmoid_err", ANYVECTOR, TRUE);
+    if (svec == NULL) return (FALSE);
+    MatchVector (svec, xvec, OPIHI_FLT);
+    opihi_flt *v1 = svec[0].elements.Flt;
+    opihi_flt *v2 = yvec[0].elements.Flt;
+    if (!strncasecmp ("poisson", argv[3], Nsvec)) {
+      for (int i = 0; i < svec[0].Nelements; i++) {
+	v1[i] = (isfinite(v2[i]) && (v2[i] > 0)) ? sqrt(v2[i]) : 1.0;
+      }
+    } else {
+      for (int i = 0; i < svec[0].Nelements; i++) {
+	v1[i] = 1.0;
+      }
+    }
+  } else {
+    if ((svec = SelectVector (argv[3], OLDVECTOR, TRUE)) == NULL) return (FALSE);
+  }
+
+  CastVector (svec, OPIHI_FLT);
+  // XXX Cast is failing.
+
+  int Npts = xvec[0].Nelements;
+  ResetVector (ovec, OPIHI_FLT, Npts);
+
+  ALLOCATE_PTR (dy, opihi_flt, Npts);
+
+  opihi_flt par[4];
+  GET_VAR (par[0], "C0");
+  GET_VAR (par[1], "C1");
+  GET_VAR (par[2], "C2");
+  GET_VAR (par[3], "C3");
+  int Npar = 4;
+
+  // mrqmin takes the inverse variance (do not generate NANs)
+  opihi_flt *v1 = svec[0].elements.Flt;
+  opihi_flt *v2 = dy;
+  for (int i = 0; i < Npts; i++, v1++, v2++) {
+      *v2 = (*v1 == 0.0) ? 0.0 : 1.0 / (*v1 * *v1);
+  } 
+  
+  opihi_flt ochisq = mrqinit (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, FUNC, !Quiet);
+  opihi_flt dchisq = ochisq + 2*Npts;
+
+  int Niter;
+  for (Niter = 0; (Niter < 20) && ((dchisq > 0.1*(Npts - Npar)) || (dchisq <= 0.0)); Niter++) {
+    opihi_flt chisq = mrqmin (xvec[0].elements.Flt, yvec[0].elements.Flt, dy, Npts, par, Npar, FUNC, !Quiet);
+    dchisq = ochisq - chisq;
+    ochisq = chisq;
+    if (!Quiet) gprint (GP_ERR, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
+  }  
+  if (!Quiet) gprint (GP_ERR, "%d iterations\n", Niter); 
+
+  for (int i = 0; i < Npts; i++) {
+    ovec[0].elements.Flt[i] = FUNC (xvec[0].elements.Flt[i], par, Npar, dy);
+  }
+  ovec[0].Nelements = Npts;
+  /* set output *before* variable renomalization */
+
+  opihi_flt **covar = mrqcovar (Npar);
+  for (int i = 0; i < Npar; i++) {
+    char name[16];
+    sprintf (name, "C%d", i);
+    set_variable (name, par[i]);
+    if (!Quiet) gprint (GP_ERR, "%d  %f  %f\n", i, par[i], sqrt(covar[i][i]));
+    sprintf (name, "dC%d", i);
+    set_variable (name, sqrt(covar[i][i]));
+  }
+
+  free (dy);
+  mrqfree (Npar);
+  return (TRUE);
+}
+
+/* pars: x_o, sigma, top, bottom */
+opihi_flt rsigmoidOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  double z = (x - par[0])/par[1];
+  double r = exp (z);
+  double s = 1 + r;
+  double q = 1 / s;
+  double f = par[2]*q + par[3];
+
+  // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0
+  // df/dq  = par[2]
+  // dq/ds  = -s^-2 = -q/s
+  // ds/dr  = 1
+  // dr/dz  = r
+  // dz/dp0 = -1 / par[1]
+  // dz/dp1 = -z / par[1]
+
+  // df/dp0 = par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1])
+  // df/dp1 = par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1])
+
+  dpar[0] = par[2]*q*r/(s*par[1]); 
+  dpar[1] = dpar[0]*z;
+  dpar[2] = q;
+  dpar[3] = 1;
+  
+  return (f);
+}
+
+/* pars: x_o, sigma, top, bottom */
+opihi_flt fsigmoidOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
+  OHANA_UNUSED_PARAM(Npar);
+
+  double z = (x - par[0])/par[1];
+  double r = exp (-z);
+  double s = 1 + r;
+  double q = 1 / s;
+  double f = par[2]*q + par[3];
+
+  // df/dp0 = df/dq * dq/dr * dr/dz * dz/dp0
+  // df/dq  = par[2]
+  // dq/ds  = -s^-2 = -q/s
+  // ds/dr  = 1
+  // dr/dz  = -r
+  // dz/dp0 = -1 / par[1]
+  // dz/dp1 = -z / par[1]
+
+  // df/dp0 = -par[2]*(-q/s)*r*(-1/par[1]) = par[2]*q*r/(s*par[1])
+  // df/dp1 = -par[2]*(-q/s)*r*(-z/par[1]) = par[2]*q*r*z/(s*par[1])
+
+  dpar[0] = -par[2]*q*r/(s*par[1]); 
+  dpar[1] = dpar[0]*z;
+  dpar[2] = q;
+  dpar[3] = 1;
+  
+  return (f);
+}
Index: /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vtransitions.c
===================================================================
--- /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vtransitions.c	(revision 40751)
+++ /branches/eam_branches/ohana.20190329/src/opihi/cmd.data/vtransitions.c	(revision 40752)
@@ -137,5 +137,5 @@
   ResetVector (outpos, OPIHI_FLT,    Nout);
 
-  if (!QUIET) gprint (GP_LOG, "found %d transitions)\n", Nout);
+  if (!QUIET) gprint (GP_LOG, "found %d transitions\n", Nout);
 
   return (TRUE);