From f83db85f6c0cafef11531ff0dc213b93e6eab04f Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 14:53:44 -0800 Subject: [PATCH 01/14] add binning --- src/array_ops.cxx | 308 +++++++++++++++++++++++++++++++++++++++++ test/test_array_ops.py | 144 ++++++++++++++++++- 2 files changed, 451 insertions(+), 1 deletion(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index 21489acc..06170cf1 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1256,6 +1256,278 @@ void detrend(bp::object & tod, const std::string & method, const int linear_ncou } } +template +void _histogram(const T* data, const T* weights, int* histogram, + const T* bin_edges, const int nsamps, const int nbins, + const T lower, const T upper) +{ + for (int i = 0; i < nsamps; ++i) { + for (int j = 0; j < nbins; ++j) { + if (data[i] < lower || data[i] > upper) + continue; + + if ((data[i] >= bin_edges[j]) && (data[i] < bin_edges[j + 1])) { + histogram[j] += weights[i]; + break; + } + } + // Edge case to match np.histogram + if (data[i] == bin_edges[nbins]) { + histogram[nbins - 1] += weights[i]; + } + } +} + +template +int _find_bin_index(const T* bin_edges, T value, int nbins) { + int left = 0; + int right = nbins; + + // Mimic np.clip and assign out-of-bounds points to + // first and last bin + if (value < bin_edges[left]) { + return 0; + } + else if (value >= bin_edges[right]) { + return right - 1; + } + + while (left < right) { + int mid = left + (right - left) / 2; + + if (value >= bin_edges[mid] && value < bin_edges[mid + 1]) { + return mid; + } + else if (value >= bin_edges[mid + 1]) { + left = mid + 1; + } + else { + right = mid; + } + } + + return -1; +} + +template +void _bin_signal(const bp::object & bin_by, const bp::object & signal, + const bp::object & weight, bp::object & binned_sig, + bp::object & binned_sig_sigma, bp::object & bin_counts, + const bp::object & bin_edges, const T lower, const T upper, + const int* flags_data=nullptr, const bool is_flags_2d=false) +{ + // signal + BufferWrapper signal_buf ("signal", signal, false, std::vector{-1, -1}); + if (signal_buf->strides[1] != signal_buf->itemsize) + throw ValueError_exception("Argument 'signal' must be contiguous in last axis."); + const int ndets = signal_buf->shape[0]; + const int nsamps = signal_buf->shape[1]; + T* signal_data = (T*)signal_buf->buf; + + // bin_by + BufferWrapper bin_by_buf ("bin_by", bin_by, false, std::vector{nsamps}); + if (bin_by_buf->strides[0] != bin_by_buf->itemsize) + throw ValueError_exception("Argument 'bin_by' must be a C-contiguous 1d array"); + T* bin_by_data = (T*)bin_by_buf->buf; + + // weight + BufferWrapper weight_buf ("weight", weight, false); + if (weight_buf->ndim == 1 && weight_buf->strides[0] != weight_buf->itemsize) + throw ValueError_exception("Argument 'weight' must be a C-contiguous 1d array"); + else if (weight_buf->ndim == 2 && weight_buf->strides[1] != weight_buf->itemsize) + throw ValueError_exception("Argument 'weight' must be contiguous in last axis."); + + T* weight_data = (T*)weight_buf->buf; + + bool is_weight_2d = false; + if (weight_buf->ndim == 2) { + is_weight_2d = true; + } + + // bin_edges + BufferWrapper bin_edges_buf ("bin_edges", bin_edges, false, std::vector{-1}); + if (bin_edges_buf->strides[0] != bin_edges_buf->itemsize) + throw ValueError_exception("Argument 'bin_edges' must be a C-contiguous 1d array"); + const int nbins = bin_edges_buf->shape[0] - 1; + T* bin_edges_data = (T*)bin_edges_buf->buf; + + // binned_sig + BufferWrapper binned_sig_buf ("binned_sig", binned_sig, false, std::vector{ndets, nbins}); + if (binned_sig_buf->strides[1] != binned_sig_buf->itemsize) + throw ValueError_exception("Argument 'binned_sig' must be contiguous in last axis."); + T* binned_sig_data = (T*)binned_sig_buf->buf; + + // binned_sig_sigma + BufferWrapper binned_sig_sigma_buf ("binned_sig_sigma", binned_sig_sigma, false, std::vector{ndets, nbins}); + if (binned_sig_sigma_buf->strides[1] != binned_sig_sigma_buf->itemsize) + throw ValueError_exception("Argument 'binned_sig_sigma' must be contiguous in last axis."); + T* binned_sig_sigma_data = (T*)binned_sig_sigma_buf->buf; + + // bin_counts + BufferWrapper bin_counts_buf ("bin_counts", bin_counts, false, std::vector{ndets, nbins}); + if (bin_counts_buf->strides[1] != bin_counts_buf->itemsize) + throw ValueError_exception("Argument 'bin_counts' must be contiguous in last axis."); + int* bin_counts_data = (int*)bin_counts_buf->buf; + + // Map from data column to bin index + T* bin_indices = (T*) malloc(nsamps * sizeof(T)); + for (int i = 0; i < nsamps; ++i) { + bin_indices[i] = _find_bin_index(bin_edges_data, bin_by_data[i], nbins); + } + + // Make the histogram up front if no flag array given + if (!flags_data) { + for (int i = 0; i < nbins; ++i) { + bin_counts_data[i] = 0; + } + _histogram(bin_by_data, weight_data, bin_counts_data, bin_edges_data, + nsamps, nbins, lower, upper); + + // Set all other detectors to first det bins if no flags + #pragma omp parallel for + for (int i = 1; i < ndets; ++i) { + int binned_ioff = i * nbins; + int* bin_counts_row = bin_counts_data + binned_ioff; + for (int j = 0; j < nbins; ++j) { + bin_counts_row[j] = bin_counts_data[j]; + } + } + } + + T* binned_sig_sq_mean = (T*) malloc(nbins * ndets * sizeof(T)); + + #pragma omp parallel for + for (int i = 0; i < ndets; ++i) { + int ioff = i * nsamps; + int binned_ioff = i * nbins; + int weight_ioff = 0; + + // Weights may be 1D or 2D + if (is_weight_2d) { + weight_ioff = ioff; + } + + T* signal_row = signal_data + ioff; + T* binned_sig_row = binned_sig_data + binned_ioff; + T* binned_sig_sq_mean_row = binned_sig_sq_mean + binned_ioff; + T* binned_sig_sigma_row = binned_sig_sigma_data + binned_ioff; + T* weight_row = weight_data + weight_ioff; + int* bin_counts_row = bin_counts_data + binned_ioff; + + // Zero out binned data + for (int j = 0; j < nbins; ++j) { + binned_sig_row[j] = 0; + binned_sig_sq_mean_row[j] = 0; + + if (flags_data) { + bin_counts_row[j] = 0; + } + } + + // Populate binned data + for (int j = 0; j < nsamps; ++j) { + bool samp_flagged = false; + if (flags_data) { + // Flags may be 1D or 2D + int flags_ioff = 0; + + if (is_flags_2d) { + flags_ioff = ioff; + } + + samp_flagged = flags_data[flags_ioff + j]; + } + if (!samp_flagged) { + int bin = bin_indices[j]; + if (flags_data) { + bin_counts_row[bin] += weight_row[j]; + } + + if (bin_counts_row[bin] > 0) { + T ws = weight_row[j] * signal_row[j]; + binned_sig_row[bin] += ws; + binned_sig_sq_mean_row[bin] += ws * ws; + } + } + } + // Normalize + for (int j = 0; j < nbins; ++j) { + if (bin_counts_row[j] > 0) { + binned_sig_row[j] /= bin_counts_row[j]; + binned_sig_sq_mean_row[j] /= bin_counts_row[j]; + binned_sig_sigma_row[j] = + std::sqrt(std::abs(binned_sig_sq_mean_row[j] - + binned_sig_row[j] * binned_sig_row[j])) / + std::sqrt(bin_counts_row[j]); + } + else { + binned_sig_row[j] = std::numeric_limits::quiet_NaN(); + binned_sig_sigma_row[j] = std::numeric_limits::quiet_NaN(); + } + } + } +} + +void bin_signal(const bp::object & bin_by, const bp::object & signal, + const bp::object & weight, bp::object & binned_sig, + bp::object & binned_sig_sigma, bp::object & bin_counts, + const bp::object & bins, const double lower, + const double upper) +{ + // Get data type + int dtype = get_dtype(signal); + + if (dtype == NPY_FLOAT) { + _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, + bin_counts, bins, (float)lower, (float)upper); + } + else if (dtype == NPY_DOUBLE) { + _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, + bin_counts, bins, (double)lower, (double)upper); + } + else { + throw TypeError_exception("Only float32 or float64 arrays are supported."); + } +} + +void bin_flagged_signal(const bp::object & bin_by, const bp::object & signal, + const bp::object & weight, bp::object & binned_sig, + bp::object & binned_sig_sigma, bp::object & bin_counts, + const bp::object & bins, const double lower, + const double upper, const bp::object & flags) +{ + // Get data type + int dtype = get_dtype(signal); + + // flags + BufferWrapper flags_buf ("flags", flags, false); + if (flags_buf->ndim == 1 && flags_buf->strides[0] != flags_buf->itemsize) + throw ValueError_exception("Argument 'flags' must be a C-contiguous 1d array"); + else if (flags_buf->ndim == 2 && flags_buf->strides[1] != flags_buf->itemsize) + throw ValueError_exception("Argument 'flags' must be contiguous in last axis."); + + bool is_flags_2d = false; + if (flags_buf->ndim == 2) { + is_flags_2d = true; + } + + int* flags_data = (int*)flags_buf->buf; + + if (dtype == NPY_FLOAT) { + _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, + bin_counts, bins, (float)lower, (float)upper, + flags_data, is_flags_2d); + } + else if (dtype == NPY_DOUBLE) { + _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, + bin_counts, bins, (double)lower, (double)upper, + flags_data, is_flags_2d); + } + else { + throw TypeError_exception("Only float32 or float64 arrays are supported."); + } +} + PYBINDINGS("so3g") { @@ -1415,4 +1687,40 @@ PYBINDINGS("so3g") " linear_ncount: Number (int) of samples to use on each end, when measuring mean level for 'linear'" " detrend. Must be a positive integer or -1. If -1, nsamps / 2 will be used. Values " " larger than 1 suppress the influence of white noise.\n"); + bp::def("bin_signal", bin_signal, + "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper)" + "\n" + "Bin time-ordered data by the ``bin_by`` and return the binned signal and its standard deviation." + "Args:\n" + " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)" + " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" + " weight: array (float32/float64) of weights for the signal values. May have shapes " + " of (nsamps) or (ndets, nsamps)" + " binned_signal: binned signal array (float32/float64) with shape (ndet,nsamp). " + " Modified in place." + " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nsamp). " + " Modified in place." + " bin_counts: counts of binned samples (int32) with shape (ndet,nsamp). Modified in place." + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Modified in place." + " lower: lower bin range (float64)" + " upper: upper bin range (float64)\n"); + bp::def("bin_flagged_signal", bin_flagged_signal, + "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper, flags)" + "\n" + "Bin time-ordered data by the ``bin_by`` and return the binned signal and its standard deviation." + "Args:\n" + " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)" + " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" + " weight: array (float32/float64) of weights for the signal values. May have shapes " + " of (nsamps) or (ndets, nsamps)" + " binned_signal: binned signal array (float32/float64) with shape (ndet,nsamp). " + " Modified in place." + " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nsamp). " + " Modified in place." + " bin_counts: counts of binned samples (int32) with shape (ndet,nsamp). Modified in place." + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Modified in place." + " lower: lower bin range (float64)" + " upper: upper bin range (float64)" + " flags: array (int32) indicating whether to exclude flagged samples when binning the signal." + " Can be of shape (nsamp) or (ndet,nsamp).\n"); } \ No newline at end of file diff --git a/test/test_array_ops.py b/test/test_array_ops.py index 20c6036a..da06613b 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -395,7 +395,7 @@ def test_01_median_detrending(self): def test_02_linear_detrending(self): nsamps = 1000 - ndets = 10 + ndets = 3 dtype = "float32" order = "C" count = nsamps // 3 @@ -425,5 +425,147 @@ def test_02_linear_detrending(self): np.testing.assert_allclose(signal_copy, signal, rtol=rtol, atol=atol) +class TestBinning(unittest.TestCase): + """ + Test binning. + """ + + def test_00_binning_no_flags_float64(self): + nsamps = 1000 + ndets = 3 + dtype = "float64" + order = "C" + + x_min = 0.0 + x_max = 1.0 + bins = 10 + + x = np.linspace(x_min, x_max, nsamps, dtype=dtype) + signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], + dtype=dtype, order=order) + weight = np.ones((nsamps), dtype=dtype, order=order) + + bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) + bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center + nbins = len(bin_centers) + + def numpy_binning(): + binned_signal = np.full([ndets, nbins], np.nan) + binned_signal_squared_mean = np.full([ndets, nbins], np.nan) + binned_signal_sigma = np.full([ndets, nbins], np.nan) + + # get bin indices + bin_indices = np.digitize(x, bin_edges) - 1 + bin_indices = np.clip(bin_indices, 0, nbins-1) + + bin_counts, _ = np.histogram(x, bins=bins, range=[x_min,x_max], weights = weight) + mcnts = bin_counts > 0 + + for i in range(ndets): + binned_signal[i][mcnts] = np.bincount(bin_indices, weights=signal[i]*weight, minlength=nbins + )[mcnts]/bin_counts[mcnts] + binned_signal_squared_mean[i][mcnts] = np.bincount(bin_indices, weights=(signal[i]*weight)**2, minlength=nbins + )[mcnts]/bin_counts[mcnts] + + binned_signal_sigma[:, mcnts] = np.sqrt(np.abs(binned_signal_squared_mean[:,mcnts] - binned_signal[:,mcnts]**2) + ) / np.sqrt(bin_counts[mcnts]) + bin_counts_dets = np.tile(bin_counts, (ndets, 1)) + + return binned_signal, binned_signal_sigma, bin_counts_dets + + + binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() + + # so3g + binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) + + so3g.bin_signal(x, signal, weight, binned_signal_so3g, binned_signal_sigma_so3g, + bin_counts_so3g, bin_edges, x_min, x_max)#, flags) + + tolerance = 1e-10 + np.testing.assert_allclose(binned_signal, binned_signal, rtol=tolerance) + np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, rtol=tolerance) + np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, rtol=tolerance) + + def test_00_binning_flags_float64(self): + nsamps = 1000 + ndets = 3 + dtype = "float64" + order = "C" + + x_min = 0.0 + x_max = 1.0 + bins = 10 + + x = np.linspace(x_min, x_max, nsamps, dtype=dtype) + signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], + dtype=dtype, order=order) + weight = np.ones((nsamps), dtype=dtype, order=order) + flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) + flags[0,::2] = 1 + flags[1,::3] = 1 + flags[2,::4] = 1 + + bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) + bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center + nbins = len(bin_centers) + + def numpy_binning(): + binned_signal = np.full([ndets, nbins], np.nan) + binned_signal_squared_mean = np.full([ndets, nbins], np.nan) + binned_signal_sigma = np.full([ndets, nbins], np.nan) + + # get bin indices + bin_indices = np.digitize(x, bin_edges) - 1 + bin_indices = np.clip(bin_indices, 0, nbins-1) + bin_counts_dets = np.full([ndets, nbins], np.nan) + + if flags.shape == (ndets, nsamps): + flag_is_2d = True + m_2d = ~flags.astype(bool) + elif flags.shape == (nsamps, ): + flag_is_2d = False + m = ~flags.astype(bool) + + for i in range(ndets): + if flag_is_2d: + m = m_2d[i] + + if weight.shape == (ndets, nsamps): + weight_det = weight[i] + elif weight.shape == (nsamps, ): + weight_det = weight + + bin_counts_dets[i] = np.bincount(bin_indices[m], weights=weight_det[m], minlength=nbins) + mcnts = bin_counts_dets[i] > 0 + binned_signal[i][mcnts] = np.bincount(bin_indices[m], weights=signal[i][m]*weight_det[m], minlength=nbins + )[mcnts]/bin_counts_dets[i][mcnts] + binned_signal_squared_mean[i][mcnts] = np.bincount(bin_indices[m], weights=(signal[i][m]*weight_det[m])**2, minlength=nbins + )[mcnts]/bin_counts_dets[i][mcnts] + binned_signal_sigma[i][mcnts] = np.sqrt(np.abs(binned_signal_squared_mean[i,mcnts] - binned_signal[i,mcnts]**2) + ) / np.sqrt(bin_counts_dets[i][mcnts]) + + return binned_signal, binned_signal_sigma, bin_counts_dets + + + binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() + + # so3g + binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) + + so3g.bin_flagged_signal(x, signal, weight, binned_signal_so3g, + binned_signal_sigma_so3g, bin_counts_so3g, + bin_edges, x_min, x_max, flags) + + tolerance = 1e-10 + np.testing.assert_allclose(binned_signal, binned_signal, rtol=tolerance) + np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, rtol=tolerance) + np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, rtol=tolerance) + + if __name__ == "__main__": unittest.main() From 3adbc81e269c3d194d6a4d263fa0a47aba1ef466 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 19:42:52 -0800 Subject: [PATCH 02/14] add strides, more dim checking, more tests --- src/array_ops.cxx | 129 +++++++++++++++++++++++++++-------------- test/test_array_ops.py | 119 ++++++++++++++++++++++++++++++------- 2 files changed, 182 insertions(+), 66 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index 06170cf1..bf57f9ab 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1297,10 +1297,10 @@ int _find_bin_index(const T* bin_edges, T value, int nbins) { if (value >= bin_edges[mid] && value < bin_edges[mid + 1]) { return mid; - } + } else if (value >= bin_edges[mid + 1]) { left = mid + 1; - } + } else { right = mid; } @@ -1314,7 +1314,9 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, const bp::object & weight, bp::object & binned_sig, bp::object & binned_sig_sigma, bp::object & bin_counts, const bp::object & bin_edges, const T lower, const T upper, - const int* flags_data=nullptr, const bool is_flags_2d=false) + const int* flags_data=nullptr, + const std::vector & flags_shape={0,0}, + const int flags_stride=0) { // signal BufferWrapper signal_buf ("signal", signal, false, std::vector{-1, -1}); @@ -1324,6 +1326,20 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, const int nsamps = signal_buf->shape[1]; T* signal_data = (T*)signal_buf->buf; + // Check flags shape + bool is_flags_2d = false; + if (flags_data) { + if (flags_shape.size() == 2) { + is_flags_2d = true; + if (flags_shape[0] != ndets || flags_shape[1] != nsamps) { + throw ValueError_exception("2D 'flags' array has incorrect shape"); + } + } + else if (flags_shape[0] != nsamps) { + throw ValueError_exception("1D 'flags' array has incorrect shape"); + } + } + // bin_by BufferWrapper bin_by_buf ("bin_by", bin_by, false, std::vector{nsamps}); if (bin_by_buf->strides[0] != bin_by_buf->itemsize) @@ -1331,18 +1347,25 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, T* bin_by_data = (T*)bin_by_buf->buf; // weight - BufferWrapper weight_buf ("weight", weight, false); - if (weight_buf->ndim == 1 && weight_buf->strides[0] != weight_buf->itemsize) + BufferWrapper weight_buf_temp ("weight", weight, true); + if (weight_buf_temp->ndim == 1 && weight_buf_temp->strides[0] != weight_buf_temp->itemsize) throw ValueError_exception("Argument 'weight' must be a C-contiguous 1d array"); - else if (weight_buf->ndim == 2 && weight_buf->strides[1] != weight_buf->itemsize) + else if (weight_buf_temp->ndim == 2 && weight_buf_temp->strides[1] != weight_buf_temp->itemsize) throw ValueError_exception("Argument 'weight' must be contiguous in last axis."); - T* weight_data = (T*)weight_buf->buf; - + // Check weight dimensions bool is_weight_2d = false; - if (weight_buf->ndim == 2) { + std::vector weight_dims; + if (weight_buf_temp->ndim == 2) { + weight_dims.push_back(ndets); is_weight_2d = true; } + if (weight_buf_temp->ndim >= 1) { + weight_dims.push_back(nsamps); + } + + BufferWrapper weight_buf ("weight", weight, false, weight_dims); + T* weight_data = (T*)weight_buf->buf; // bin_edges BufferWrapper bin_edges_buf ("bin_edges", bin_edges, false, std::vector{-1}); @@ -1369,56 +1392,67 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, throw ValueError_exception("Argument 'bin_counts' must be contiguous in last axis."); int* bin_counts_data = (int*)bin_counts_buf->buf; + // Strides + int signal_stride = signal_buf->strides[0] / sizeof(T); + int weight_stride = 0; + if (is_weight_2d) { + weight_stride = weight_buf->strides[0] / sizeof(T); + } + int binned_sig_stride = binned_sig_buf->strides[0] / sizeof(T); + int binned_sig_sigma_stride = binned_sig_sigma_buf->strides[0] / sizeof(T); + int bin_counts_stride = bin_counts_buf->strides[0] / sizeof(int); + // Map from data column to bin index T* bin_indices = (T*) malloc(nsamps * sizeof(T)); for (int i = 0; i < nsamps; ++i) { bin_indices[i] = _find_bin_index(bin_edges_data, bin_by_data[i], nbins); } - + // Make the histogram up front if no flag array given + // Assumes weights is 1D if (!flags_data) { for (int i = 0; i < nbins; ++i) { bin_counts_data[i] = 0; } _histogram(bin_by_data, weight_data, bin_counts_data, bin_edges_data, nsamps, nbins, lower, upper); - + // Set all other detectors to first det bins if no flags #pragma omp parallel for for (int i = 1; i < ndets; ++i) { - int binned_ioff = i * nbins; + int binned_ioff = i * bin_counts_stride; int* bin_counts_row = bin_counts_data + binned_ioff; for (int j = 0; j < nbins; ++j) { bin_counts_row[j] = bin_counts_data[j]; } } } - + T* binned_sig_sq_mean = (T*) malloc(nbins * ndets * sizeof(T)); - + #pragma omp parallel for for (int i = 0; i < ndets; ++i) { - int ioff = i * nsamps; - int binned_ioff = i * nbins; + int ioff = i * signal_stride; + int binned_ioff = i * bin_counts_stride; int weight_ioff = 0; // Weights may be 1D or 2D if (is_weight_2d) { - weight_ioff = ioff; + weight_ioff = i * weight_stride; } - + T* signal_row = signal_data + ioff; - T* binned_sig_row = binned_sig_data + binned_ioff; - T* binned_sig_sq_mean_row = binned_sig_sq_mean + binned_ioff; - T* binned_sig_sigma_row = binned_sig_sigma_data + binned_ioff; - T* weight_row = weight_data + weight_ioff; + T* binned_sig_row = binned_sig_data + (i * binned_sig_stride); + T* binned_sig_sq_mean_row = binned_sig_sq_mean + (i * nbins); + T* binned_sig_sigma_row = binned_sig_sigma_data + (i * binned_sig_sigma_stride); + T* weight_row = weight_data + weight_stride; int* bin_counts_row = bin_counts_data + binned_ioff; - + // Zero out binned data for (int j = 0; j < nbins; ++j) { binned_sig_row[j] = 0; binned_sig_sq_mean_row[j] = 0; - + if (flags_data) { bin_counts_row[j] = 0; } @@ -1430,9 +1464,9 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, if (flags_data) { // Flags may be 1D or 2D int flags_ioff = 0; - + if (is_flags_2d) { - flags_ioff = ioff; + flags_ioff = i * flags_stride; } samp_flagged = flags_data[flags_ioff + j]; @@ -1442,7 +1476,7 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, if (flags_data) { bin_counts_row[bin] += weight_row[j]; } - + if (bin_counts_row[bin] > 0) { T ws = weight_row[j] * signal_row[j]; binned_sig_row[bin] += ws; @@ -1455,7 +1489,7 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, if (bin_counts_row[j] > 0) { binned_sig_row[j] /= bin_counts_row[j]; binned_sig_sq_mean_row[j] /= bin_counts_row[j]; - binned_sig_sigma_row[j] = + binned_sig_sigma_row[j] = std::sqrt(std::abs(binned_sig_sq_mean_row[j] - binned_sig_row[j] * binned_sig_row[j])) / std::sqrt(bin_counts_row[j]); @@ -1476,7 +1510,7 @@ void bin_signal(const bp::object & bin_by, const bp::object & signal, { // Get data type int dtype = get_dtype(signal); - + if (dtype == NPY_FLOAT) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bins, (float)lower, (float)upper); @@ -1506,22 +1540,25 @@ void bin_flagged_signal(const bp::object & bin_by, const bp::object & signal, else if (flags_buf->ndim == 2 && flags_buf->strides[1] != flags_buf->itemsize) throw ValueError_exception("Argument 'flags' must be contiguous in last axis."); - bool is_flags_2d = false; + std::vector flags_shape; + flags_shape.push_back(flags_buf->shape[0]); + if (flags_buf->ndim == 2) { - is_flags_2d = true; + flags_shape.push_back(flags_buf->shape[1]); } - + int* flags_data = (int*)flags_buf->buf; - + int flags_stride = flags_buf->strides[0] / sizeof(int); + if (dtype == NPY_FLOAT) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bins, (float)lower, (float)upper, - flags_data, is_flags_2d); + flags_data, flags_shape, flags_stride); } else if (dtype == NPY_DOUBLE) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bins, (double)lower, (double)upper, - flags_data, is_flags_2d); + flags_data, flags_shape, flags_stride); } else { throw TypeError_exception("Only float32 or float64 arrays are supported."); @@ -1690,35 +1727,37 @@ PYBINDINGS("so3g") bp::def("bin_signal", bin_signal, "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper)" "\n" - "Bin time-ordered data by the ``bin_by`` and return the binned signal and its standard deviation." + "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation. " + "This function uses OMP to parallelize over the dets (rows) axis.\n" "Args:\n" " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)" " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" " weight: array (float32/float64) of weights for the signal values. May have shapes " " of (nsamps) or (ndets, nsamps)" - " binned_signal: binned signal array (float32/float64) with shape (ndet,nsamp). " + " binned_signal: binned signal array (float32/float64) with shape (ndet,nbin). " " Modified in place." - " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nsamp). " + " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin). " " Modified in place." - " bin_counts: counts of binned samples (int32) with shape (ndet,nsamp). Modified in place." - " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Modified in place." + " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place." + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing." " lower: lower bin range (float64)" " upper: upper bin range (float64)\n"); bp::def("bin_flagged_signal", bin_flagged_signal, "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper, flags)" "\n" - "Bin time-ordered data by the ``bin_by`` and return the binned signal and its standard deviation." + "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation. " + "This function uses OMP to parallelize over the dets (rows) axis.\n" "Args:\n" " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)" " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" " weight: array (float32/float64) of weights for the signal values. May have shapes " " of (nsamps) or (ndets, nsamps)" - " binned_signal: binned signal array (float32/float64) with shape (ndet,nsamp). " + " binned_signal: binned signal array (float32/float64) with shape (ndet,nbin). " " Modified in place." - " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nsamp). " + " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin). " " Modified in place." - " bin_counts: counts of binned samples (int32) with shape (ndet,nsamp). Modified in place." - " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Modified in place." + " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place." + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing." " lower: lower bin range (float64)" " upper: upper bin range (float64)" " flags: array (int32) indicating whether to exclude flagged samples when binning the signal." diff --git a/test/test_array_ops.py b/test/test_array_ops.py index da06613b..d1551368 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -453,20 +453,20 @@ def numpy_binning(): binned_signal = np.full([ndets, nbins], np.nan) binned_signal_squared_mean = np.full([ndets, nbins], np.nan) binned_signal_sigma = np.full([ndets, nbins], np.nan) - + # get bin indices bin_indices = np.digitize(x, bin_edges) - 1 bin_indices = np.clip(bin_indices, 0, nbins-1) - + bin_counts, _ = np.histogram(x, bins=bins, range=[x_min,x_max], weights = weight) mcnts = bin_counts > 0 - - for i in range(ndets): + + for i in range(ndets): binned_signal[i][mcnts] = np.bincount(bin_indices, weights=signal[i]*weight, minlength=nbins )[mcnts]/bin_counts[mcnts] binned_signal_squared_mean[i][mcnts] = np.bincount(bin_indices, weights=(signal[i]*weight)**2, minlength=nbins )[mcnts]/bin_counts[mcnts] - + binned_signal_sigma[:, mcnts] = np.sqrt(np.abs(binned_signal_squared_mean[:,mcnts] - binned_signal[:,mcnts]**2) ) / np.sqrt(bin_counts[mcnts]) bin_counts_dets = np.tile(bin_counts, (ndets, 1)) @@ -475,7 +475,7 @@ def numpy_binning(): binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() - + # so3g binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) @@ -485,14 +485,14 @@ def numpy_binning(): bin_counts_so3g, bin_edges, x_min, x_max)#, flags) tolerance = 1e-10 - np.testing.assert_allclose(binned_signal, binned_signal, rtol=tolerance) - np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, rtol=tolerance) - np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, rtol=tolerance) + np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) + np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) - def test_00_binning_flags_float64(self): + def test_01_binning_flags_float32(self): nsamps = 1000 ndets = 3 - dtype = "float64" + dtype = "float32" order = "C" x_min = 0.0 @@ -502,7 +502,7 @@ def test_00_binning_flags_float64(self): x = np.linspace(x_min, x_max, nsamps, dtype=dtype) signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], dtype=dtype, order=order) - weight = np.ones((nsamps), dtype=dtype, order=order) + weight = np.ones((ndets, nsamps), dtype=dtype, order=order) flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) flags[0,::2] = 1 flags[1,::3] = 1 @@ -516,28 +516,105 @@ def numpy_binning(): binned_signal = np.full([ndets, nbins], np.nan) binned_signal_squared_mean = np.full([ndets, nbins], np.nan) binned_signal_sigma = np.full([ndets, nbins], np.nan) - + + # get bin indices + bin_indices = np.digitize(x, bin_edges) - 1 + bin_indices = np.clip(bin_indices, 0, nbins-1) + bin_counts_dets = np.full([ndets, nbins], np.nan) + + if flags.shape == (ndets, nsamps): + flag_is_2d = True + m_2d = ~flags.astype(bool) + elif flags.shape == (nsamps, ): + flag_is_2d = False + m = ~flags.astype(bool) + + for i in range(ndets): + if flag_is_2d: + m = m_2d[i] + + if weight.shape == (ndets, nsamps): + weight_det = weight[i] + elif weight.shape == (nsamps, ): + weight_det = weight + + bin_counts_dets[i] = np.bincount(bin_indices[m], weights=weight_det[m], minlength=nbins) + mcnts = bin_counts_dets[i] > 0 + binned_signal[i][mcnts] = np.bincount(bin_indices[m], weights=signal[i][m]*weight_det[m], minlength=nbins + )[mcnts]/bin_counts_dets[i][mcnts] + binned_signal_squared_mean[i][mcnts] = np.bincount(bin_indices[m], weights=(signal[i][m]*weight_det[m])**2, minlength=nbins + )[mcnts]/bin_counts_dets[i][mcnts] + binned_signal_sigma[i][mcnts] = np.sqrt(np.abs(binned_signal_squared_mean[i,mcnts] - binned_signal[i,mcnts]**2) + ) / np.sqrt(bin_counts_dets[i][mcnts]) + + return binned_signal, binned_signal_sigma, bin_counts_dets + + + binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() + + # so3g + binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) + + so3g.bin_flagged_signal(x, signal, weight, binned_signal_so3g, + binned_signal_sigma_so3g, bin_counts_so3g, + bin_edges, x_min, x_max, flags) + + tolerance = 1e-4 + np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) + np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) + + def test_02_binning_flags_float64(self): + nsamps = 1000 + ndets = 3 + dtype = "float64" + order = "C" + + x_min = 0.0 + x_max = 1.0 + bins = 10 + + x = np.linspace(x_min, x_max, nsamps, dtype=dtype) + signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], + dtype=dtype, order=order) + weight = np.ones((ndets, nsamps), dtype=dtype, order=order) + flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) + # flags[0,::2] = 1 + # flags[1,::3] = 1 + # flags[2,::4] = 1 + + bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) + bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center + nbins = len(bin_centers) + + def numpy_binning(): + binned_signal = np.full([ndets, nbins], np.nan) + binned_signal_squared_mean = np.full([ndets, nbins], np.nan) + binned_signal_sigma = np.full([ndets, nbins], np.nan) + # get bin indices bin_indices = np.digitize(x, bin_edges) - 1 bin_indices = np.clip(bin_indices, 0, nbins-1) bin_counts_dets = np.full([ndets, nbins], np.nan) - + if flags.shape == (ndets, nsamps): flag_is_2d = True m_2d = ~flags.astype(bool) elif flags.shape == (nsamps, ): flag_is_2d = False m = ~flags.astype(bool) - + for i in range(ndets): if flag_is_2d: m = m_2d[i] - + if weight.shape == (ndets, nsamps): weight_det = weight[i] elif weight.shape == (nsamps, ): weight_det = weight - + bin_counts_dets[i] = np.bincount(bin_indices[m], weights=weight_det[m], minlength=nbins) mcnts = bin_counts_dets[i] > 0 binned_signal[i][mcnts] = np.bincount(bin_indices[m], weights=signal[i][m]*weight_det[m], minlength=nbins @@ -551,7 +628,7 @@ def numpy_binning(): binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() - + # so3g binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) @@ -562,9 +639,9 @@ def numpy_binning(): bin_edges, x_min, x_max, flags) tolerance = 1e-10 - np.testing.assert_allclose(binned_signal, binned_signal, rtol=tolerance) - np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, rtol=tolerance) - np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, rtol=tolerance) + np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) + np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) if __name__ == "__main__": From 8908f58e8b7b4dd123cb79a0e045475530afdf99 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 20:22:56 -0800 Subject: [PATCH 03/14] fix bug with range --- src/array_ops.cxx | 4 ++-- test/test_array_ops.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index bf57f9ab..74bce1c5 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1272,7 +1272,7 @@ void _histogram(const T* data, const T* weights, int* histogram, } } // Edge case to match np.histogram - if (data[i] == bin_edges[nbins]) { + if (data[i] == bin_edges[nbins] && data[i] <= upper) { histogram[nbins - 1] += weights[i]; } } @@ -1762,4 +1762,4 @@ PYBINDINGS("so3g") " upper: upper bin range (float64)" " flags: array (int32) indicating whether to exclude flagged samples when binning the signal." " Can be of shape (nsamp) or (ndet,nsamp).\n"); -} \ No newline at end of file +} diff --git a/test/test_array_ops.py b/test/test_array_ops.py index d1551368..953b514d 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -431,14 +431,14 @@ class TestBinning(unittest.TestCase): """ def test_00_binning_no_flags_float64(self): - nsamps = 1000 - ndets = 3 + nsamps = 10 + ndets = 1 dtype = "float64" order = "C" x_min = 0.0 x_max = 1.0 - bins = 10 + bins = 2 x = np.linspace(x_min, x_max, nsamps, dtype=dtype) signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], @@ -482,7 +482,7 @@ def numpy_binning(): bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) so3g.bin_signal(x, signal, weight, binned_signal_so3g, binned_signal_sigma_so3g, - bin_counts_so3g, bin_edges, x_min, x_max)#, flags) + bin_counts_so3g, bin_edges, x_min, x_max) tolerance = 1e-10 np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) From 3cd9fffad4638a65cc516dcbf4e187b28d44e02e Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 20:24:20 -0800 Subject: [PATCH 04/14] undo debug changes --- test/test_array_ops.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_array_ops.py b/test/test_array_ops.py index 953b514d..feb4c287 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -431,14 +431,14 @@ class TestBinning(unittest.TestCase): """ def test_00_binning_no_flags_float64(self): - nsamps = 10 - ndets = 1 + nsamps = 1000 + ndets = 3 dtype = "float64" order = "C" x_min = 0.0 x_max = 1.0 - bins = 2 + bins = 10 x = np.linspace(x_min, x_max, nsamps, dtype=dtype) signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], From 6d6f92707dd2cf25adc2f3f346316e7f2f2767e1 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 20:29:05 -0800 Subject: [PATCH 05/14] undo more debug changes --- test/test_array_ops.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_array_ops.py b/test/test_array_ops.py index feb4c287..8fd7f87e 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -240,7 +240,7 @@ def test_02_linear_extrapolation(self): t_interp_end = 1009.0 t_interp_size = 2000 - ndet = 3 + ndet = 10 dtype = "float32" order = "C" @@ -581,9 +581,9 @@ def test_02_binning_flags_float64(self): dtype=dtype, order=order) weight = np.ones((ndets, nsamps), dtype=dtype, order=order) flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) - # flags[0,::2] = 1 - # flags[1,::3] = 1 - # flags[2,::4] = 1 + flags[0,::2] = 1 + flags[1,::3] = 1 + flags[2,::4] = 1 bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center From cea537f7e30f373feb1a99efeec7f59fa769df87 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 20:38:48 -0800 Subject: [PATCH 06/14] some doc changes, rename bins --- src/array_ops.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index 74bce1c5..c19fd9e7 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1505,7 +1505,7 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, void bin_signal(const bp::object & bin_by, const bp::object & signal, const bp::object & weight, bp::object & binned_sig, bp::object & binned_sig_sigma, bp::object & bin_counts, - const bp::object & bins, const double lower, + const bp::object & bin_edges, const double lower, const double upper) { // Get data type @@ -1513,11 +1513,11 @@ void bin_signal(const bp::object & bin_by, const bp::object & signal, if (dtype == NPY_FLOAT) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, - bin_counts, bins, (float)lower, (float)upper); + bin_counts, bin_edges, (float)lower, (float)upper); } else if (dtype == NPY_DOUBLE) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, - bin_counts, bins, (double)lower, (double)upper); + bin_counts, bin_edges, (double)lower, (double)upper); } else { throw TypeError_exception("Only float32 or float64 arrays are supported."); @@ -1527,7 +1527,7 @@ void bin_signal(const bp::object & bin_by, const bp::object & signal, void bin_flagged_signal(const bp::object & bin_by, const bp::object & signal, const bp::object & weight, bp::object & binned_sig, bp::object & binned_sig_sigma, bp::object & bin_counts, - const bp::object & bins, const double lower, + const bp::object & bin_edges, const double lower, const double upper, const bp::object & flags) { // Get data type @@ -1552,12 +1552,12 @@ void bin_flagged_signal(const bp::object & bin_by, const bp::object & signal, if (dtype == NPY_FLOAT) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, - bin_counts, bins, (float)lower, (float)upper, + bin_counts, bin_edges, (float)lower, (float)upper, flags_data, flags_shape, flags_stride); } else if (dtype == NPY_DOUBLE) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, - bin_counts, bins, (double)lower, (double)upper, + bin_counts, bin_edges, (double)lower, (double)upper, flags_data, flags_shape, flags_stride); } else { @@ -1734,7 +1734,7 @@ PYBINDINGS("so3g") " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" " weight: array (float32/float64) of weights for the signal values. May have shapes " " of (nsamps) or (ndets, nsamps)" - " binned_signal: binned signal array (float32/float64) with shape (ndet,nbin). " + " binned_sig: binned signal array (float32/float64) with shape (ndet,nbin). " " Modified in place." " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin). " " Modified in place." @@ -1752,7 +1752,7 @@ PYBINDINGS("so3g") " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" " weight: array (float32/float64) of weights for the signal values. May have shapes " " of (nsamps) or (ndets, nsamps)" - " binned_signal: binned signal array (float32/float64) with shape (ndet,nbin). " + " binned_sig: binned signal array (float32/float64) with shape (ndet,nbin). " " Modified in place." " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin). " " Modified in place." From 7fe9585b5b9e2d9f0f522bcbefb2b9279fd35f82 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 20:48:57 -0800 Subject: [PATCH 07/14] undo accidental changes --- test/test_array_ops.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_array_ops.py b/test/test_array_ops.py index 8fd7f87e..a80b69cd 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -240,7 +240,7 @@ def test_02_linear_extrapolation(self): t_interp_end = 1009.0 t_interp_size = 2000 - ndet = 10 + ndet = 3 dtype = "float32" order = "C" @@ -395,7 +395,7 @@ def test_01_median_detrending(self): def test_02_linear_detrending(self): nsamps = 1000 - ndets = 3 + ndets = 10 dtype = "float32" order = "C" count = nsamps // 3 From a57a92ec51e771a197ea196ca9b37eaad6bc1e57 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 20:50:59 -0800 Subject: [PATCH 08/14] fix typo --- test/test_array_ops.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/test_array_ops.py b/test/test_array_ops.py index a80b69cd..64446863 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -485,7 +485,7 @@ def numpy_binning(): bin_counts_so3g, bin_edges, x_min, x_max) tolerance = 1e-10 - np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) @@ -562,7 +562,7 @@ def numpy_binning(): bin_edges, x_min, x_max, flags) tolerance = 1e-4 - np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) @@ -639,7 +639,7 @@ def numpy_binning(): bin_edges, x_min, x_max, flags) tolerance = 1e-10 - np.testing.assert_allclose(binned_signal, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) From 955ddcba69ad2561e8b490348a192fd9482d3d7c Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Thu, 27 Feb 2025 21:19:05 -0800 Subject: [PATCH 09/14] fix docstring format --- src/array_ops.cxx | 52 +++++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 26 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index c19fd9e7..4a4f51a0 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1727,39 +1727,39 @@ PYBINDINGS("so3g") bp::def("bin_signal", bin_signal, "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper)" "\n" - "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation. " + "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation.\n" "This function uses OMP to parallelize over the dets (rows) axis.\n" "Args:\n" - " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)" - " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" - " weight: array (float32/float64) of weights for the signal values. May have shapes " - " of (nsamps) or (ndets, nsamps)" - " binned_sig: binned signal array (float32/float64) with shape (ndet,nbin). " - " Modified in place." - " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin). " + " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)\n" + " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)\n" + " weight: array (float32/float64) of weights for the signal values. May have shapes\n" + " of (nsamps) or (ndets, nsamps)\n" + " binned_sig: binned signal array (float32/float64) with shape (ndet,nbin).\n" + " Modified in place.\n" + " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin).\n" " Modified in place." - " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place." - " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing." - " lower: lower bin range (float64)" + " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing.\n" + " lower: lower bin range (float64)\n" " upper: upper bin range (float64)\n"); bp::def("bin_flagged_signal", bin_flagged_signal, - "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper, flags)" + "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper, flags)\n" "\n" - "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation. " + "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation.\n" "This function uses OMP to parallelize over the dets (rows) axis.\n" "Args:\n" - " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)" - " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)" - " weight: array (float32/float64) of weights for the signal values. May have shapes " - " of (nsamps) or (ndets, nsamps)" - " binned_sig: binned signal array (float32/float64) with shape (ndet,nbin). " - " Modified in place." - " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin). " - " Modified in place." - " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place." - " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing." - " lower: lower bin range (float64)" - " upper: upper bin range (float64)" - " flags: array (int32) indicating whether to exclude flagged samples when binning the signal." + " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)\n" + " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)\n" + " weight: array (float32/float64) of weights for the signal values. May have shapes\n" + " of (nsamps) or (ndets, nsamps)\n" + " binned_sig: binned signal array (float32/float64) with shape (ndet,nbin).\n" + " Modified in place.\n" + " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin).\n" + " Modified in place.\n" + " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing.\n" + " lower: lower bin range (float64)\n" + " upper: upper bin range (float64)\n" + " flags: array (int32) indicating whether to exclude flagged samples when binning the signal.\n" " Can be of shape (nsamp) or (ndet,nsamp).\n"); } From 874ae70dc0743a53c60de5c772667062521bc742 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Fri, 28 Feb 2025 06:27:43 -0800 Subject: [PATCH 10/14] fix bin index type, stride cleanup --- src/array_ops.cxx | 20 +++++--------------- 1 file changed, 5 insertions(+), 15 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index 4a4f51a0..6649a8a9 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1403,7 +1403,7 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, int bin_counts_stride = bin_counts_buf->strides[0] / sizeof(int); // Map from data column to bin index - T* bin_indices = (T*) malloc(nsamps * sizeof(T)); + int* bin_indices = (int*) malloc(nsamps * sizeof(int)); for (int i = 0; i < nsamps; ++i) { bin_indices[i] = _find_bin_index(bin_edges_data, bin_by_data[i], nbins); } @@ -1434,12 +1434,7 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, for (int i = 0; i < ndets; ++i) { int ioff = i * signal_stride; int binned_ioff = i * bin_counts_stride; - int weight_ioff = 0; - - // Weights may be 1D or 2D - if (is_weight_2d) { - weight_ioff = i * weight_stride; - } + int weight_ioff = i * weight_stride; T* signal_row = signal_data + ioff; T* binned_sig_row = binned_sig_data + (i * binned_sig_stride); @@ -1462,13 +1457,7 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, for (int j = 0; j < nsamps; ++j) { bool samp_flagged = false; if (flags_data) { - // Flags may be 1D or 2D - int flags_ioff = 0; - - if (is_flags_2d) { - flags_ioff = i * flags_stride; - } - + int flags_ioff = i * flags_stride; samp_flagged = flags_data[flags_ioff + j]; } if (!samp_flagged) { @@ -1543,12 +1532,13 @@ void bin_flagged_signal(const bp::object & bin_by, const bp::object & signal, std::vector flags_shape; flags_shape.push_back(flags_buf->shape[0]); + int flags_stride = 0; if (flags_buf->ndim == 2) { flags_shape.push_back(flags_buf->shape[1]); + flags_stride = flags_buf->strides[0] / sizeof(int); } int* flags_data = (int*)flags_buf->buf; - int flags_stride = flags_buf->strides[0] / sizeof(int); if (dtype == NPY_FLOAT) { _bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, From e022ced59f41c38b125add9d9b61ace492180083 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Mon, 3 Mar 2025 07:12:14 -0800 Subject: [PATCH 11/14] clarify allowing different bin widths in docstring --- src/array_ops.cxx | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index 6649a8a9..fc3d4302 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1718,7 +1718,7 @@ PYBINDINGS("so3g") "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper)" "\n" "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation.\n" - "This function uses OMP to parallelize over the dets (rows) axis.\n" + "This function uses OMP to parallelize over the dets (rows) axis. Supports unequal bin widths.\n" "Args:\n" " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)\n" " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)\n" @@ -1729,14 +1729,15 @@ PYBINDINGS("so3g") " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin).\n" " Modified in place." " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" - " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing.\n" + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing but\n" + " but may have different widths.\n" " lower: lower bin range (float64)\n" " upper: upper bin range (float64)\n"); bp::def("bin_flagged_signal", bin_flagged_signal, "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper, flags)\n" "\n" "Bin time-ordered data by ``bin_by`` and return the binned signal and its standard deviation.\n" - "This function uses OMP to parallelize over the dets (rows) axis.\n" + "This function uses OMP to parallelize over the dets (rows) axis. Supports unequal bin widths.\n" "Args:\n" " bin_by: the array (float32/float64) by which signal is binned with shape (nsamp)\n" " signal: the signal array (float32/float64) to be binned with shape (ndet,nsamp)\n" @@ -1747,7 +1748,8 @@ PYBINDINGS("so3g") " binned_sig_sigma: estimated sigma of binned signal (float32/float64) with shape (ndet,nbin).\n" " Modified in place.\n" " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" - " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing.\n" + " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing but\n" + " but may have different widths.\n" " lower: lower bin range (float64)\n" " upper: upper bin range (float64)\n" " flags: array (int32) indicating whether to exclude flagged samples when binning the signal.\n" From 112c5e120bc6cfc69e7cdf9abbe52089252bfdb5 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Wed, 9 Jul 2025 07:54:36 -0700 Subject: [PATCH 12/14] updates and fixes --- src/array_ops.cxx | 39 +++----- test/test_array_ops.py | 204 +++++++++++++++++------------------------ 2 files changed, 95 insertions(+), 148 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index fc3d4302..4033c00d 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1256,28 +1256,6 @@ void detrend(bp::object & tod, const std::string & method, const int linear_ncou } } -template -void _histogram(const T* data, const T* weights, int* histogram, - const T* bin_edges, const int nsamps, const int nbins, - const T lower, const T upper) -{ - for (int i = 0; i < nsamps; ++i) { - for (int j = 0; j < nbins; ++j) { - if (data[i] < lower || data[i] > upper) - continue; - - if ((data[i] >= bin_edges[j]) && (data[i] < bin_edges[j + 1])) { - histogram[j] += weights[i]; - break; - } - } - // Edge case to match np.histogram - if (data[i] == bin_edges[nbins] && data[i] <= upper) { - histogram[nbins - 1] += weights[i]; - } - } -} - template int _find_bin_index(const T* bin_edges, T value, int nbins) { int left = 0; @@ -1408,14 +1386,16 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, bin_indices[i] = _find_bin_index(bin_edges_data, bin_by_data[i], nbins); } - // Make the histogram up front if no flag array given - // Assumes weights is 1D if (!flags_data) { for (int i = 0; i < nbins; ++i) { bin_counts_data[i] = 0; } - _histogram(bin_by_data, weight_data, bin_counts_data, bin_edges_data, - nsamps, nbins, lower, upper); + for (int i = 0; i < nsamps; ++i) { + if (bin_by[i] < lower || bin_by[i] > upper) + continue; + int bin = bin_indices[i]; + bin_counts_data[bin] += weight_data[i]; + } // Set all other detectors to first det bins if no flags #pragma omp parallel for @@ -1489,6 +1469,9 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, } } } + + free(bin_indices); + free(binned_sig_sq_mean); } void bin_signal(const bp::object & bin_by, const bp::object & signal, @@ -1750,8 +1733,8 @@ PYBINDINGS("so3g") " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing but\n" " but may have different widths.\n" - " lower: lower bin range (float64)\n" - " upper: upper bin range (float64)\n" + " lower: lower bin range. Data points falling outside this range will be ignored. (float64)\n" + " upper: upper bin range. Data points falling outside this range will be ignored. (float64)\n" " flags: array (int32) indicating whether to exclude flagged samples when binning the signal.\n" " Can be of shape (nsamp) or (ndet,nsamp).\n"); } diff --git a/test/test_array_ops.py b/test/test_array_ops.py index 64446863..e38987eb 100644 --- a/test/test_array_ops.py +++ b/test/test_array_ops.py @@ -429,36 +429,43 @@ class TestBinning(unittest.TestCase): """ Test binning. """ - - def test_00_binning_no_flags_float64(self): - nsamps = 1000 - ndets = 3 - dtype = "float64" - order = "C" - - x_min = 0.0 - x_max = 1.0 - bins = 10 - + # make fake input data + def make_input(self, nsamps, ndets, x_min, x_max, dtype, order="C"): x = np.linspace(x_min, x_max, nsamps, dtype=dtype) signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], dtype=dtype, order=order) weight = np.ones((nsamps), dtype=dtype, order=order) - bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) + return x, signal, weight + + # allocate so3g containers + def make_so3g(self, nsamps, ndets, nbins, dtype, order): + binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) + bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) + + return binned_signal_so3g, binned_signal_sigma_so3g, bin_counts_so3g + + def get_bins(self, x, bins, ranges): + bin_edges = np.histogram_bin_edges(x, bins=bins, range=ranges,) bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center nbins = len(bin_centers) - def numpy_binning(): - binned_signal = np.full([ndets, nbins], np.nan) - binned_signal_squared_mean = np.full([ndets, nbins], np.nan) - binned_signal_sigma = np.full([ndets, nbins], np.nan) + return bin_edges, nbins + + # equivalent to sotodlib binning function + def numpy_binning(self, x, signal, weight, bin_edges, nbins, ranges, flags=None): + ndets, nsamps = signal.shape + binned_signal = np.full([ndets, nbins], np.nan) + binned_signal_squared_mean = np.full([ndets, nbins], np.nan) + binned_signal_sigma = np.full([ndets, nbins], np.nan) - # get bin indices - bin_indices = np.digitize(x, bin_edges) - 1 - bin_indices = np.clip(bin_indices, 0, nbins-1) + # get bin indices + bin_indices = np.digitize(x, bin_edges) - 1 + bin_indices = np.clip(bin_indices, 0, nbins-1) - bin_counts, _ = np.histogram(x, bins=bins, range=[x_min,x_max], weights = weight) + if flags is None: + bin_counts, _ = np.histogram(x, bins=nbins, range=ranges, weights = weight) mcnts = bin_counts > 0 for i in range(ndets): @@ -473,53 +480,7 @@ def numpy_binning(): return binned_signal, binned_signal_sigma, bin_counts_dets - - binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() - - # so3g - binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) - binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) - bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) - - so3g.bin_signal(x, signal, weight, binned_signal_so3g, binned_signal_sigma_so3g, - bin_counts_so3g, bin_edges, x_min, x_max) - - tolerance = 1e-10 - np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) - np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) - np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) - - def test_01_binning_flags_float32(self): - nsamps = 1000 - ndets = 3 - dtype = "float32" - order = "C" - - x_min = 0.0 - x_max = 1.0 - bins = 10 - - x = np.linspace(x_min, x_max, nsamps, dtype=dtype) - signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], - dtype=dtype, order=order) - weight = np.ones((ndets, nsamps), dtype=dtype, order=order) - flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) - flags[0,::2] = 1 - flags[1,::3] = 1 - flags[2,::4] = 1 - - bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) - bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center - nbins = len(bin_centers) - - def numpy_binning(): - binned_signal = np.full([ndets, nbins], np.nan) - binned_signal_squared_mean = np.full([ndets, nbins], np.nan) - binned_signal_sigma = np.full([ndets, nbins], np.nan) - - # get bin indices - bin_indices = np.digitize(x, bin_edges) - 1 - bin_indices = np.clip(bin_indices, 0, nbins-1) + else: bin_counts_dets = np.full([ndets, nbins], np.nan) if flags.shape == (ndets, nsamps): @@ -549,94 +510,97 @@ def numpy_binning(): return binned_signal, binned_signal_sigma, bin_counts_dets + def test_00_binning_no_flags_float64(self): + nsamps = 1000 + ndets = 3 + dtype = "float64" + order = "C" + + x_min = 0.0 + x_max = 1.0 + bins = 10 + ranges = (x_min, x_max) - binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() + x, signal, weight = self.make_input(nsamps, ndets, x_min, x_max, dtype, order="C") + bin_edges, nbins = self.get_bins(x, bins, ranges) + # numpy binning + binned_signal, binned_signal_sigma, bin_counts_dets = self.numpy_binning(x, signal, weight, bin_edges, nbins, ranges) # so3g - binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) - binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) - bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) + binned_signal_so3g, binned_signal_sigma_so3g, bin_counts_so3g = self.make_so3g(nsamps, ndets, nbins, dtype, order) - so3g.bin_flagged_signal(x, signal, weight, binned_signal_so3g, - binned_signal_sigma_so3g, bin_counts_so3g, - bin_edges, x_min, x_max, flags) + so3g.bin_signal(x, signal, weight, binned_signal_so3g, binned_signal_sigma_so3g, + bin_counts_so3g, bin_edges, ranges[0], ranges[1]) - tolerance = 1e-4 + tolerance = 1e-10 np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) - def test_02_binning_flags_float64(self): + def test_01_binning_flags_float32(self): nsamps = 1000 ndets = 3 - dtype = "float64" + dtype = "float32" order = "C" x_min = 0.0 x_max = 1.0 bins = 10 + ranges = (x_min, x_max) + + x, signal, weight = self.make_input(nsamps, ndets, x_min, x_max, + dtype, order="C") - x = np.linspace(x_min, x_max, nsamps, dtype=dtype) - signal = np.array([(i + 1) * np.sin(2*np.pi*x + i) for i in range(ndets)], - dtype=dtype, order=order) - weight = np.ones((ndets, nsamps), dtype=dtype, order=order) flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) flags[0,::2] = 1 flags[1,::3] = 1 flags[2,::4] = 1 - bin_edges = np.histogram_bin_edges(x, bins=bins, range=[x_min,x_max],) - bin_centers = (bin_edges[1] - bin_edges[0])/2. + bin_edges[:-1] # edge to center - nbins = len(bin_centers) - - def numpy_binning(): - binned_signal = np.full([ndets, nbins], np.nan) - binned_signal_squared_mean = np.full([ndets, nbins], np.nan) - binned_signal_sigma = np.full([ndets, nbins], np.nan) + bin_edges, nbins = self.get_bins(x, bins, ranges) + # numpy binning + binned_signal, binned_signal_sigma, bin_counts_dets = self.numpy_binning(x, signal, weight, bin_edges, nbins, ranges, flags) - # get bin indices - bin_indices = np.digitize(x, bin_edges) - 1 - bin_indices = np.clip(bin_indices, 0, nbins-1) - bin_counts_dets = np.full([ndets, nbins], np.nan) + # so3g + binned_signal_so3g, binned_signal_sigma_so3g, bin_counts_so3g = self.make_so3g(nsamps, ndets, nbins, dtype, order) - if flags.shape == (ndets, nsamps): - flag_is_2d = True - m_2d = ~flags.astype(bool) - elif flags.shape == (nsamps, ): - flag_is_2d = False - m = ~flags.astype(bool) + so3g.bin_flagged_signal(x, signal, weight, binned_signal_so3g, + binned_signal_sigma_so3g, bin_counts_so3g, + bin_edges, ranges[0], ranges[1], flags) - for i in range(ndets): - if flag_is_2d: - m = m_2d[i] + tolerance = 1e-4 + np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) + np.testing.assert_allclose(binned_signal_sigma_so3g, binned_signal_sigma, atol=tolerance) + np.testing.assert_allclose(bin_counts_so3g, bin_counts_dets, atol=tolerance) - if weight.shape == (ndets, nsamps): - weight_det = weight[i] - elif weight.shape == (nsamps, ): - weight_det = weight + def test_02_binning_flags_float64(self): + nsamps = 1000 + ndets = 3 + dtype = "float64" + order = "C" - bin_counts_dets[i] = np.bincount(bin_indices[m], weights=weight_det[m], minlength=nbins) - mcnts = bin_counts_dets[i] > 0 - binned_signal[i][mcnts] = np.bincount(bin_indices[m], weights=signal[i][m]*weight_det[m], minlength=nbins - )[mcnts]/bin_counts_dets[i][mcnts] - binned_signal_squared_mean[i][mcnts] = np.bincount(bin_indices[m], weights=(signal[i][m]*weight_det[m])**2, minlength=nbins - )[mcnts]/bin_counts_dets[i][mcnts] - binned_signal_sigma[i][mcnts] = np.sqrt(np.abs(binned_signal_squared_mean[i,mcnts] - binned_signal[i,mcnts]**2) - ) / np.sqrt(bin_counts_dets[i][mcnts]) + x_min = 0.0 + x_max = 1.0 + bins = 10 + ranges = (x_min, x_max) - return binned_signal, binned_signal_sigma, bin_counts_dets + x, signal, weight = self.make_input(nsamps, ndets, x_min, x_max, + dtype, order="C") + flags = np.zeros((ndets, nsamps), dtype=np.int32, order=order) + flags[0,::2] = 1 + flags[1,::3] = 1 + flags[2,::4] = 1 - binned_signal, binned_signal_sigma, bin_counts_dets = numpy_binning() + bin_edges, nbins = self.get_bins(x, bins, ranges) + # numpy binning + binned_signal, binned_signal_sigma, bin_counts_dets = self.numpy_binning(x, signal, weight, bin_edges, nbins, ranges, flags) # so3g - binned_signal_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) - binned_signal_sigma_so3g = np.zeros((ndets, nbins), dtype=dtype, order=order) - bin_counts_so3g = np.zeros((ndets, nbins), dtype=np.int32, order=order) + binned_signal_so3g, binned_signal_sigma_so3g, bin_counts_so3g = self.make_so3g(nsamps, ndets, nbins, dtype, order) so3g.bin_flagged_signal(x, signal, weight, binned_signal_so3g, binned_signal_sigma_so3g, bin_counts_so3g, - bin_edges, x_min, x_max, flags) + bin_edges, ranges[0], ranges[1], flags) tolerance = 1e-10 np.testing.assert_allclose(binned_signal_so3g, binned_signal, atol=tolerance) From a74675d962d76083795d7adb34854ebd04027939 Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Wed, 9 Jul 2025 08:11:20 -0700 Subject: [PATCH 13/14] fix docstring --- src/array_ops.cxx | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index f82e0096..fa555337 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1718,8 +1718,8 @@ PYBINDINGS("so3g") " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing but\n" " but may have different widths.\n" - " lower: lower bin range (float64)\n" - " upper: upper bin range (float64)\n"); + " lower: lower bin range (float64). Data points falling outside this range will be ignored.\n" + " upper: upper bin range (float64). Data points falling outside this range will be ignored.\n"); bp::def("bin_flagged_signal", bin_flagged_signal, "bin_signal(bin_by, signal, weight, binned_sig, binned_sig_sigma, bin_counts, bin_edges, lower, upper, flags)\n" "\n" @@ -1737,8 +1737,8 @@ PYBINDINGS("so3g") " bin_counts: counts of binned samples (int32) with shape (ndet,nbin). Modified in place.\n" " bin_edges: array (float32/float64) of bin edges with length=nbins+1. Must be monotonically increasing but\n" " but may have different widths.\n" - " lower: lower bin range. Data points falling outside this range will be ignored. (float64)\n" - " upper: upper bin range. Data points falling outside this range will be ignored. (float64)\n" + " lower: lower bin range (float64). Data points falling outside this range will be ignored.\n" + " upper: upper bin range (float64). Data points falling outside this range will be ignored.\n" " flags: array (int32) indicating whether to exclude flagged samples when binning the signal.\n" " Can be of shape (nsamp) or (ndet,nsamp).\n"); } From b8dd35d8842610e539c7a9d103b6d13a7b5253df Mon Sep 17 00:00:00 2001 From: Michael McCrackan Date: Wed, 9 Jul 2025 08:31:13 -0700 Subject: [PATCH 14/14] add comment on up-front bin count calc --- src/array_ops.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/array_ops.cxx b/src/array_ops.cxx index fa555337..d4dd17a4 100644 --- a/src/array_ops.cxx +++ b/src/array_ops.cxx @@ -1390,6 +1390,8 @@ void _bin_signal(const bp::object & bin_by, const bp::object & signal, bin_indices[i] = _find_bin_index(bin_edges_data, bin_by_data[i], nbins); } + // Populate bin_counts_data up front if no flag array given since its + // faster if (!flags_data) { for (int i = 0; i < nbins; ++i) { bin_counts_data[i] = 0;