Merge pull request #76 from hamogu/debug_fef_search

Propagate changes to binary search

jdmath/src/dinterpo.c

@@ -44,32 +44,38 @@ unsigned int JDMbinary_search_d (double x, double *xp, unsigned int n)
 	  }
 	else n0 = n2;
      }
+   if (x >= xp[n0]) return n1;
    return n0;
 }
 
 double JDMinterpolate_d (double x, double *xp, double *yp, unsigned int n)
 {
-   unsigned int n0, n1;
-   double x0, x1;
-
-   n0 = JDMbinary_search_d (x, xp, n);
-
-   x0 = xp[n0];
-   n1 = n0 + 1;
-
-   if (x == x0)
-     return yp[n0];
-   if (n1 == n)
-     {
-	if (n0 == 0)
-	  return yp[n0];
-	n1 = n0 - 1;
-     }
-
-   x1 = xp[n1];
-   if (x1 == x0) return yp[n0];
-
-   return yp[n0] + (yp[n1] - yp[n0]) / (x1 - x0) * (x - x0);
+	unsigned int n0, n1;
+	double x0, x1;
+
+	if (n == 1)
+		return yp[0];
+	n1 = JDMbinary_search_f(x, xp, n);
+	n0 = n1 - 1;
+
+	if (x == xp[n1])
+		return yp[n1];
+	if (n1 == n)
+	{
+		n1--;
+		n0--;
+	}
+	if (n1 == 0)
+	{
+		n0 = 1;
+	}
+
+	x0 = xp[n0];
+	x1 = xp[n1];
+	if (x1 == x0)
+		return yp[n1];
+
+	return yp[n0] + (yp[n1] - yp[n0]) / (x1 - x0) * (x - x0);
 }
 
 int JDMinterpolate_dvector (double *xp, double *yp, unsigned int n,

jdmath/src/finterpo.c

@@ -27,6 +27,13 @@
 
 #include "jdmath.h"
 
+/*
+ * Binary Search
+ *
+ * This routine returns the index of the first element in the array xp
+ * that is greater than or equal to x.  If all elements are less than x,
+ * then n is returned.
+ */
 unsigned int JDMbinary_search_f (float x, float *xp, unsigned int n)
 {
    unsigned int n0, n1, n2;
@@ -42,10 +49,7 @@ unsigned int JDMbinary_search_f (float x, float *xp, unsigned int n)
 	     if (xp[n2] == x) return n2;
 	     n1 = n2;
 	  }
-	else
-	  {
-	     n0 = n2;
-	  }
+	else n0 = n2;
      }
 
    if (x >= xp[n0]) return n1;
@@ -57,22 +61,25 @@ float JDMinterpolate_f (float x, float *xp, float *yp, unsigned int n)
    unsigned int n0, n1;
    double x0, x1;
 
-   n0 = JDMbinary_search_f (x, xp, n);
-
-   x0 = xp[n0];
-   n1 = n0 + 1;
+   if (n == 1) return yp[0];
+   n1 = JDMbinary_search_f (x, xp, n);
+   n0 = n1 - 1;
 
-   if (x == x0)
-     return yp[n0];
+   if (x == xp[n1])
+	   return yp[n1];
    if (n1 == n)
      {
-	if (n0 == 0)
-	  return yp[n0];
-	n1 = n0 - 1;
+	    n1--;
+		n0--;
      }
+   if (n1 == 0)
+	 {
+		 n0 = 1;
+	 }
 
+   x0 = xp[n0];
    x1 = xp[n1];
-   if (x1 == x0) return yp[n0];
+   if (x1 == x0) return yp[n1];
 
    return yp[n0] + (yp[n1] - yp[n0]) / (x1 - x0) * (x - x0);
 }
@@ -192,14 +199,15 @@ float JDMlog_interpolate_f (float x, float *xp, float *yp, unsigned int n)
    unsigned int n0, n1;
    double x0, x1;
 
-   n0 = JDMbinary_search_f (x, xp, n);
+   n1 = JDMbinary_search_f (x, xp, n);
+   n0 = n1 - 1;
 
-   x0 = xp[n0];
-   n1 = n0 + 1;
-
-   if ((x == x0) || (n1 == n)) return yp[n0];
+   if (n1 == 0) return yp[0];
+   if (n1 == n) return yp[n - 1];
 
+   x0 = xp[n0];
    x1 = xp[n1];
+   if (x == x1) return yp[n1];
    if (x1 == x0) return yp[n0];
 
    return yp[n0] + (yp[n1] - yp[n0]) * (log(x/x0) / log (x1/x0));

jdmath/src/hist.c

@@ -40,34 +40,14 @@
 #define FLOAT_TYPE float
 #include "hist.c"
 #else
-#if 0
-static unsigned int binary_search_d (double x, double *xp, unsigned int n)
-{
-   unsigned int n0, n1, n2;
-
-   n0 = 0;
-   n1 = n;
-
-   while (n1 > n0 + 1)
-     {
-	n2 = (n0 + n1) / 2;
-	if (xp[n2] >= x)
-	  {
-	     if (xp[n2] == x) return n2;
-	     n1 = n2;
-	  }
-	else n0 = n2;
-     }
-   return n0;
-}
-#endif
 /*
  * histogram routines
  *
  * A 1-d histogram is specified by a set of N grid points X_k and some set
  * of values y_i to be grouped into the histogram.  The bin size of the
  * nth bin is given by x_{i+1} - x_i, except for the last bin, which is
  * assumed to be of infinite width.
+ * This means that values of pts that are less than x_0 are ignored.
  */
 
 /* If reverse_indices is NON-NULL, it is assumed to point to an array of
@@ -102,10 +82,10 @@ int JDMHISTOGRAM_FUNCTION (FLOAT_TYPE *pts, unsigned int npts,
 	  continue;
 
 	j = BINARY_SEARCH (val, bin_edges, nbins);
-	histogram[j] += 1;
+	histogram[j - 1] += 1;
 
 	if (reverse_indices != NULL)
-	  reverse_indices[i] = (int) j;
+	  reverse_indices[i] = (int) j - 1;
      }
    return 0;
 }

marx/doc/examples/user-source/eventlist.c

@@ -594,7 +594,9 @@ static double map_pi_to_line (short pi)
 	return b->energies[0];
      }
    i = JDMbinary_search_f (JDMrandom (), b->cum_strengths, b->num_energies);
-   return b->energies[i];
+   if (i == 0)
+     return b->energies[0];
+   return b->energies[i - 1];
 }
 
 static int read_data (double *tp, double *xp, double *yp, short *chanp)

marx/libsrc/acis_fef.c

@@ -553,19 +553,20 @@ static int normalize_gaussians (Gauss_Parm_Type *gaussians, unsigned int num,
 	/* return 0; */
      }
 
-   if (total_neg_area != 0.0)
-     flags |= HAS_NEG_AMP_GAUSSIANS;
-
-   g = gaussians;
-   if (total_pos_area > 0) while (g < gmax)
-     {
-	/* Since we interpolate over pairs of gaussians, do not normalize
-	 * each member of the pair separately.
-	 */
-	/* g->amp /= total_pos_area; */
-	g->cum_area /= total_pos_area;
-	g++;
-     }
+     if (total_neg_area != 0.0)
+       flags |= HAS_NEG_AMP_GAUSSIANS;
+
+     g = gaussians;
+     if (total_pos_area > 0)
+       while (g < gmax)
+       {
+         /* Since we interpolate over pairs of gaussians, do not normalize
+          * each member of the pair separately.
+          */
+         /* g->amp /= total_pos_area; */
+         g->cum_area /= total_pos_area;
+         g++;
+       }
 
    total_area = total_pos_area - total_neg_area;
    if (total_area != 0)
@@ -990,18 +991,22 @@ int marx_apply_acis_rmf (int ccd_id, float x, float y,
      }
 
    num_energies = f->num_energies;
-   i = JDMbinary_search_f (energy, f->energies, num_energies);
-   j = i + 1;
+   j = JDMbinary_search_f (energy, f->energies, num_energies);
+   if (j == 0)
+   {
+     /* extrapolate */
+     j++;
+   }
    if (j == num_energies)
-     {
-	/* extrapolate */
-	i--;
-	j--;
-     }
+   {
+	   /* extrapolate */
+     j--;
+   }
 
-   t = (energy - f->energies[i])/(f->energies[j] - f->energies[i]);
+   t = (energy - f->energies[j-1])/(f->energies[j] - f->energies[j-1]);
    g0 = f->gaussians + i * num_gaussians;
    g1 = g0 + num_gaussians;
+   //marx_message("acis_fef: t = %f, energy= %f, f->energies[j-1] = %f, f->energies[j] = %f \n", t, energy, f->energies[j-1], f->energies[j]);
 
    for (i = 0; i < num_gaussians; i++)
      {
@@ -1031,10 +1036,10 @@ int marx_apply_acis_rmf (int ccd_id, float x, float y,
 
    if (flags & HAS_TOTAL_NEG_AREA)
      {
-	marx_message ("Warning occurred for ccdid=%d,chipx=%f,chipy=%f,energy=%f\n",
-		      ccd_id, x, y, energy);
-	return -1;
-	/* flags &= ~HAS_TOTAL_NEG_AREA; */
+       marx_message("Warning occurred for ccdid=%d,chipx=%f,chipy=%f,energy=%f,flags=%d,flagsandneg=%d\n",
+                    ccd_id, x, y, energy, flags, flags & HAS_TOTAL_NEG_AREA);
+       return -1;
+       /* flags &= ~HAS_TOTAL_NEG_AREA; */
      }
 
    if (flags == 0)

marx/libsrc/acis_subpix.c

@@ -284,7 +284,7 @@ int marx_compute_acis_subpix (Marx_Subpix_Table_Type *stt,
 			      int ccd, float energy, int fltgrade, float *dxp, float *dyp)
 {
    Subpix_Type *s;
-   unsigned int i, j, n;
+   unsigned int j, n;
    double w0, w1;
 
    if (stt == NULL)
@@ -301,18 +301,22 @@ int marx_compute_acis_subpix (Marx_Subpix_Table_Type *stt,
      }
 
    n = (unsigned int) s->num_energies;
-   i = JDMbinary_search_f (energy, s->energies, n);
-   j = i+1;
+   j = JDMbinary_search_f (energy, s->energies, n);
+   if (j == 0)
+   {
+     /* extrapolate */
+     j++;
+   }
    if (j == n)
-     {
-	j = i;
-	i--;
-     }
+   {
+     /* extrapolate */
+     j--;
+   }
    // Linear interpolation on the energy grid to get dxp, dyp
-   w1 = ((double)energy - s->energies[i])/(s->energies[j] - s->energies[i]);
+   w1 = ((double)energy - s->energies[j - 1])/(s->energies[j] - s->energies[j - 1]);
    w0 = 1.0 - w1;
-   *dxp = w0 * s->dxs[i] + w1 * s->dxs[j];
-   *dyp = w0 * s->dys[i] + w1 * s->dys[j];
+   *dxp = w0 * s->dxs[j - 1] + w1 * s->dxs[j];
+   *dyp = w0 * s->dys[j - 1] + w1 * s->dys[j];
 
    return 0;
 }

marx/libsrc/aciscontam.c

@@ -240,18 +240,26 @@ static void compute_tau_values (Single_Component_Contam_Type *contam,
    double t = _Marx_TStart_MJDsecs;
    double t0, t1;
 
-   n0 = JDMbinary_search_d (t, times, ntimes);
-   n1 = n0+1;
+   n1 = JDMbinary_search_d (t, times, ntimes);
+   n0 = n1 - 1;
 
-   if (n1 == ntimes)
-     {
-	if (n0 == 0)
+   if (ntimes == 1)
 	  {
 	     contam->tau_0s[layer] = tau0s[0];
 	     contam->tau_1s[layer] = tau1s[0];
 	     return;
 	  }
-	n1 = n0-1;
+
+  /* Out of range. Extrapolate. */
+  if (n1 == ntimes)
+     {
+      n1--;
+      n0--;
+     }
+  if (n1 == 0)
+     {
+      n1 = 1;
+      n0 = 0;
      }
 
    t0 = times[n0];