1/* -*- C -*- */ 2/* 3 * A 1-d histogram is specified by a set of N grid points X_k and some set 4 * of values y_i to be grouped into the histogram. The bin size of the 5 * nth bin is given by x_{i+1} - x_i, except for the last bin, which is 6 * assumed to be of infinite width. 7 */ 8 9/* If reverse_indices is NON-NULL, it is assumed to point to an array of 10 * size npts. 11 * 12 * This routines assumes that the caller has initialized the histogram 13 * array to 0, and the reverse_indices array to -1. 14 */ 15 16#ifdef HISTOGRAM_1D 17static int HISTOGRAM_1D (PTS_TYPE *pts, SLuindex_Type npts, 18 BIN_EDGES_TYPE *bin_edges, SLuindex_Type nbins, 19 HistData_Type *histogram, 20 SLindex_Type *reverse_indices) 21{ 22 SLuindex_Type i, nbins_m1; 23 double xlo, xhi, dx; 24 25 if (nbins == 0) 26 return 0; 27 28 if (-1 == check_grid (bin_edges, nbins)) 29 return -1; 30 31 nbins_m1 = nbins - 1; 32 xlo = (double) bin_edges[0]; 33 xhi = (double) bin_edges[nbins_m1]; 34 dx = xhi - xlo; 35 36 if (dx < 0.0) 37 { 38 SLang_verror (SL_INVALID_PARM, "hist1d: bin edges array is not in increasing order"); 39 return -1; 40 } 41 42 for (i = 0; i < npts; i++) 43 { 44 PTS_TYPE val = pts[i]; 45 SLuindex_Type j; 46 47 if ( 48#if CHECK_NANS 49 isnan(val) || 50#endif 51 (val < xlo)) 52 continue; 53 54 if (val >= xhi) 55 j = nbins_m1; 56 else 57 { 58 /* Try linear interpolation since many grids will be linear */ 59 j = (SLuindex_Type) (((val - xlo)/dx)*nbins_m1); 60 if (j == nbins_m1) j--; 61 if ((bin_edges[j] > val) || (val >= bin_edges[j+1])) 62 j = BINARY_SEARCH (val, bin_edges, nbins); 63 } 64 histogram[j] += 1; 65 if (reverse_indices != NULL) 66 reverse_indices[i] = (SLindex_Type) j; 67 } 68 return 0; 69} 70#undef HISTOGRAM_1D 71#endif /* HISTOGRAM_1D */ 72 73#ifdef HISTOGRAM_2D 74static int HISTOGRAM_2D (PTS_TYPE *xpts, PTS_TYPE *ypts, SLuindex_Type npts, 75 BIN_EDGES_TYPE *xbin_edges, SLuindex_Type nxbins, 76 BIN_EDGES_TYPE *ybin_edges, SLuindex_Type nybins, 77 HistData_Type *histogram, 78 SLindex_Type *reverse_indices) 79{ 80 SLuindex_Type i, nxbins_m1, nybins_m1; 81 double xlo, xhi, dx; 82 double ylo, yhi, dy; 83 84 if ((nxbins == 0) || (nybins == 0)) 85 return 0; 86 87 if (-1 == check_grid (xbin_edges, nxbins)) 88 return -1; 89 90 if (-1 == check_grid (ybin_edges, nybins)) 91 return -1; 92 93 nxbins_m1 = nxbins - 1; 94 xlo = (double) xbin_edges[0]; 95 xhi = (double) xbin_edges[nxbins_m1]; 96 dx = xhi - xlo; 97 98 nybins_m1 = nybins - 1; 99 ylo = (double) ybin_edges[0]; 100 yhi = (double) ybin_edges[nybins_m1]; 101 dy = yhi - ylo; 102 103 if ((dx < 0.0) || (dy < 0.0)) 104 { 105 SLang_verror (SL_INVALID_PARM, "hist2d: bin edges array is not in increasing order"); 106 return -1; 107 } 108 109 for (i = 0; i < npts; i++) 110 { 111 PTS_TYPE xval = xpts[i]; 112 PTS_TYPE yval = ypts[i]; 113 SLuindex_Type jx, jy, j; 114 115 if ( 116#if CHECK_NANS 117 isnan(xval) || isnan(xval) || 118#endif 119 (xval < xlo) || (yval < ylo)) 120 continue; 121 122 if (xval >= xhi) 123 jx = nxbins_m1; 124 else 125 { 126 /* Try linear interpolation since many grids will be linear */ 127 jx = (SLuindex_Type) (((xval - xlo)/dx)*nxbins_m1); 128 if (jx == nxbins_m1) jx--; 129 if ((xbin_edges[jx] > xval) || (xval >= xbin_edges[jx+1])) 130 jx = BINARY_SEARCH (xval, xbin_edges, nxbins); 131 } 132 133 if (yval >= yhi) 134 jy = nybins_m1; 135 else 136 { 137 /* Try linear interpolation since many grids will be linear */ 138 jy = (SLuindex_Type) (((yval - ylo)/dy)*nybins_m1); 139 if (jy == nybins_m1) jy--; 140 if ((ybin_edges[jy] > yval) || (yval >= ybin_edges[jy+1])) 141 jy = BINARY_SEARCH (yval, ybin_edges, nybins); 142 } 143 144 j = jx*nybins + jy; 145 histogram[j] += 1; 146 if (reverse_indices != NULL) 147 reverse_indices[i] = (int) j; 148 } 149 return 0; 150} 151#undef HISTOGRAM_2D 152#endif 153 154#ifdef HISTOGRAM_REBIN 155static int HISTOGRAM_REBIN (double *new_grid, SLuindex_Type new_n, 156 double *old_grid, PTS_TYPE *old_h, SLuindex_Type old_n, 157 PTS_TYPE *new_h) 158{ 159 SLuindex_Type i, imax; 160 SLuindex_Type j, jmax; 161 double *old_left, *old_right, *new_left, *new_right; 162 163 if ((new_n == 0) || (old_n == 0)) 164 return 0; 165 166 for (i = 0; i < new_n; i++) 167 new_h[i] = 0; 168 169 old_left = old_grid; 170 old_right = old_grid + 1; 171 new_left = new_grid; 172 new_right = new_grid + 1; 173 174 imax = new_n - 1; 175 jmax = old_n - 1; 176 177 if (-1 == check_grid (old_grid, old_n)) 178 return -1; 179 180 if (-1 == check_grid (new_grid, new_n)) 181 return -1; 182 183 i = j = 0; 184 if (j < jmax) 185 { 186 double r_new, l_new, r_old, l_old; 187 double h_rho; 188 189 r_old = old_right[0]; 190 l_old = old_left [0]; 191 l_new = new_left [0]; 192 if (i == imax) 193 r_new = old_left[jmax]; 194 else 195 r_new = new_right[0]; 196 197 if (r_old > l_old) 198 h_rho = old_h[j] / (r_old - l_old); 199 else 200 h_rho = 0.0; 201 202 while (1) 203 { 204 if (r_new < r_old) 205 { 206 if (l_new >= l_old) 207 new_h[i] += h_rho * (r_new - l_new); 208 else if (r_new > l_old) 209 new_h[i] += h_rho * (r_new - l_old); 210 211 if (i != imax) 212 { 213 i++; 214 l_new = r_new; 215 if (i == imax) 216 r_new = old_left[jmax]; 217 else 218 r_new = new_right[i]; 219 } 220 } 221 else 222 { 223 if (l_new < l_old) 224 new_h[i] += old_h [j]; 225 else if (l_new < r_old) 226 new_h[i] += h_rho * (r_old - l_new); 227 228 j++; 229 if (j == jmax) 230 break; 231 l_old = r_old; 232 r_old = old_right[j]; 233 if (r_old > l_old) 234 h_rho = old_h[j] / (r_old - l_old); 235 else 236 h_rho = 0.0; 237 } 238 } 239 } 240 241 /* Now take care of the last bin, which is of infinite length. Because of 242 * its infinite length, it gets all of the stuff in the input histograms 243 * last bin. 244 */ 245 new_h [imax] += old_h [jmax]; 246 return 0; 247} 248#undef HISTOGRAM_REBIN 249#endif 250 251#undef PTS_TYPE 252#undef BIN_EDGES_TYPE 253#undef BINARY_SEARCH 254#undef CHECK_NANS 255