Skip to content

Commit

Permalink
Merge branch 'master' into detrend
Browse files Browse the repository at this point in the history
  • Loading branch information
@mmccrackan
mmccrackan authored Jan 16, 2025
2 parents 98a4e9a + 695c1bd commit b7b4e09
Show file tree
Hide file tree
Showing 7 changed files with 364 additions and 60 deletions.
4 changes: 2 additions & 2 deletions python/proj/wcs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -766,15 +766,15 @@ def compute_nside_tile(self, assembly, nActivePerThread=5, nThreads=None):
"""
if self.nside_tile == 'auto':
## Estimate fsky
nside_tile0 = 4 # ntile = 192, for estimating fsky
nside_tile0 = min(4, self.nside) # ntile = 192, for estimating fsky
self.nside_tile = nside_tile0
nActive = len(self.get_active_tiles(assembly)['active_tiles'])
fsky = nActive / (12 * nside_tile0**2)
if nThreads is None:
nThreads = so3g.useful_info()['omp_num_threads']
# nside_tile is smallest power of 2 satisfying nTile >= nActivePerThread * nthread / fsky
self.nside_tile = int(2**np.ceil(0.5 * np.log2(nActivePerThread * nThreads / (12 * fsky))))
# print('Setting nside_tile =', self.nside_tile)
self.nside_tile = min(self.nside_tile, self.nside)
return self.nside_tile

def get_active_tiles(self, assembly, assign=False):
Expand Down
5 changes: 5 additions & 0 deletions src/Projection.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -596,6 +596,11 @@ class Pixelizor_Healpix<Tiled> {
npix_per_tile = npix / ntiles;
check_nside(nside); // check validity
check_nside(nside_tile);
if (nside_tile > nside) {
std::ostringstream err;
err << "Invalid nside_tile " << nside_tile << " > nside " << nside;
throw ValueError_exception(err.str());
}
}
~Pixelizor_Healpix() {};

Expand Down
165 changes: 156 additions & 9 deletions src/array_ops.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include <vector>
#include <string>
#include <cmath>
#include <algorithm>
extern "C" {
#include <cblas.h>
// Additional prototypes for Fortran LAPACK routines.
Expand Down Expand Up @@ -637,6 +638,64 @@ void block_moment(const bp::object & tod, const bp::object & out, int bsize, int
_block_moment(tod_data, output, bsize, moment, central, ndet, nsamp, shift);
}

template <typename T>
void _minmax(T* data, T* output, int mode, int start, int stop)
{
T val;
if(mode == 0){ // get the min
val = *(std::min_element(data+start, data+stop));
}
else if(mode == 1){ // get the max
val = *(std::max_element(data+start, data+stop));
}
else{ // get the peak to peak
auto min = std::min_element(data+start, data+stop);
auto max = std::max_element(data+start, data+stop);
val = *max - *min;
}
for(int si = start; si < stop; si++) {
output[si] = val;
}

}

template <typename T>
void _block_minmax(T* tod_data, T* output, int bsize, int mode, int ndet, int nsamp, int shift)
{
int nblock = ((nsamp - shift)/bsize) + 1;
#pragma omp parallel for
for(int di = 0; di < ndet; di++)
{
int ioff = di*nsamp;
// do the the pre-shift portion
if(shift > 0){
_minmax(tod_data, output, mode, ioff, ioff+shift);
}

for(int bi = 0; bi < nblock; bi++) {
int start = (bi * bsize) + shift;
int stop = min(start + bsize, nsamp);
_minmax(tod_data, output, mode, ioff+start, ioff+stop);
}
}
}

template <typename T>
void block_minmax(const bp::object & tod, const bp::object & out, int bsize, int mode, int shift)
{
BufferWrapper<T> tod_buf ("tod", tod, false, std::vector<int>{-1, -1});
int ndet = tod_buf->shape[0];
int nsamp = tod_buf->shape[1];
T* tod_data = (T*)tod_buf->buf;
if (tod_buf->strides[1] != tod_buf->itemsize || tod_buf->strides[0] != tod_buf->itemsize*nsamp)
throw buffer_exception("tod must be C-contiguous along last axis");
BufferWrapper<T> out_buf ("out", out, false, std::vector<int>{ndet, nsamp});
if (out_buf->strides[1] != out_buf->itemsize || out_buf->strides[0] != out_buf->itemsize*nsamp)
throw buffer_exception("out must be C-contiguous along last axis");
T* output = (T*)out_buf->buf;
_block_minmax(tod_data, output, bsize, mode, ndet, nsamp, shift);
}


void _clean_flag(int* flag_data, int width, int ndet, int nsamp)
{
Expand Down Expand Up @@ -826,6 +885,53 @@ void find_quantized_jumps(const bp::object & tod, const bp::object & out, const
}
}

template <typename T>
void subtract_jump_heights(const bp::object & tod, const bp::object & out, const bp::object & heights, const bp::object & jumps) {
BufferWrapper<T> tod_buf ("tod", tod, false, std::vector<int>{-1, -1});
int ndet = tod_buf->shape[0];
int nsamp = tod_buf->shape[1];
if (tod_buf->strides[1] != tod_buf->itemsize || tod_buf->strides[0] != tod_buf->itemsize*nsamp)
throw buffer_exception("tod must be C-contiguous along last axis");
T* tod_data = (T*)tod_buf->buf;
BufferWrapper<T> out_buf ("out", out, false, std::vector<int>{ndet, nsamp});
if (out_buf->strides[1] != out_buf->itemsize || out_buf->strides[0] != out_buf->itemsize*nsamp)
throw buffer_exception("out must be C-contiguous along last axis");
T* output = (T*)out_buf->buf;
BufferWrapper<T> h_buf ("heights", heights, false, std::vector<int>{ndet, nsamp});
if (h_buf->strides[1] != h_buf->itemsize || h_buf->strides[0] != h_buf->itemsize*nsamp)
throw buffer_exception("heights must be C-contiguous along last axis");
T* jump_heights = (T*)h_buf->buf;
auto ranges = extract_ranges<int32_t>(jumps);

#pragma omp parallel for
for(int di = 0; di < ranges.size(); di++) {
int start = 0;
int stop = 0;
T min_h;
T max_h;
T height;
T to_sub = 0;
for (auto const &r: ranges[di].segments) {
start = di*nsamp + r.first;
for(int j = stop; j < start && to_sub != 0; j++) {
output[j] = tod_data[j] - to_sub;
}
stop = di*nsamp + r.second;
min_h = *(std::min_element(jump_heights+start, jump_heights+stop));
max_h = *(std::max_element(jump_heights+start, jump_heights+stop));
// Decide whether this is a negative or positive jump.
height = (abs(min_h) > abs(max_h)) ? min_h : max_h;
to_sub = to_sub + height;
for(int j = start; j < stop && to_sub != 0; j++) {
output[j] = tod_data[j] - to_sub;
}
}
for(int j = stop; j < di*nsamp + nsamp && to_sub != 0; j++) {
output[j] = tod_data[j] - to_sub;
}
}
}

template <typename T>
using _interp_func_pointer = void (*)(const double* x, const double* y,
const double* x_interp, T* y_interp,
Expand Down Expand Up @@ -903,7 +1009,7 @@ void _interp1d(const bp::object & x, const bp::object & y, const bp::object & x_
gsl_interp_accel* acc = gsl_interp_accel_alloc();
gsl_spline* spline = gsl_spline_alloc(interp_type, n_x);

#pragma omp parallel for
#pragma omp for
for (int row = 0; row < n_rows; ++row) {

int y_row_start = row * y_data_stride;
Expand All @@ -913,10 +1019,10 @@ void _interp1d(const bp::object & x, const bp::object & y, const bp::object & x_
T* y_row = y_data + y_row_start;
T* y_interp_row = y_interp_data + y_interp_row_start;

interp_func(x_data, y_row, x_interp_data, y_interp_row,
interp_func(x_data, y_row, x_interp_data, y_interp_row,
n_x, n_x_interp, spline, acc);
}

// Free gsl objects
gsl_spline_free(spline);
gsl_interp_accel_free(acc);
Expand All @@ -933,7 +1039,7 @@ void _interp1d(const bp::object & x, const bp::object & y, const bp::object & x_
std::transform(x_data, x_data + n_x, x_dbl,
[](float value) { return static_cast<double>(value); });

std::transform(x_interp_data, x_interp_data + n_x_interp, x_interp_dbl,
std::transform(x_interp_data, x_interp_data + n_x_interp, x_interp_dbl,
[](float value) { return static_cast<double>(value); });

#pragma omp parallel
Expand All @@ -942,7 +1048,7 @@ void _interp1d(const bp::object & x, const bp::object & y, const bp::object & x_
gsl_interp_accel* acc = gsl_interp_accel_alloc();
gsl_spline* spline = gsl_spline_alloc(interp_type, n_x);

#pragma omp parallel for
#pragma omp for
for (int row = 0; row < n_rows; ++row) {

int y_row_start = row * y_data_stride;
Expand All @@ -952,13 +1058,13 @@ void _interp1d(const bp::object & x, const bp::object & y, const bp::object & x_
// Transform y row to double array for gsl
double y_dbl[n_x];

std::transform(y_data + y_row_start, y_data + y_row_end, y_dbl,
std::transform(y_data + y_row_start, y_data + y_row_end, y_dbl,
[](float value) { return static_cast<double>(value); });

T* y_interp_row = y_interp_data + y_interp_row_start;

// Don't copy y_interp to doubles as it is cast during assignment
interp_func(x_dbl, y_dbl, x_interp_dbl, y_interp_row,
interp_func(x_dbl, y_dbl, x_interp_dbl, y_interp_row,
n_x, n_x_interp, spline, acc);
}

Expand All @@ -980,15 +1086,15 @@ void interp1d_linear(const bp::object & x, const bp::object & y,
const gsl_interp_type* interp_type = gsl_interp_linear;
// Pointer to interpolation function
_interp_func_pointer<float> interp_func = &_linear_interp<float>;

_interp1d<float>(x, y, x_interp, y_interp, interp_type, interp_func);
}
else if (dtype == NPY_DOUBLE) {
// GSL interpolation type
const gsl_interp_type* interp_type = gsl_interp_linear;
// Pointer to interpolation function
_interp_func_pointer<double> interp_func = &_linear_interp<double>;

_interp1d<double>(x, y, x_interp, y_interp, interp_type, interp_func);
}
else {
Expand Down Expand Up @@ -1183,6 +1289,7 @@ PYBINDINGS("so3g")
"Args:\n"
" tod: data array (float32) with shape (ndet, nsamp)\n"
" out: output array (float32) with shape (ndet, nsamp)\n"
" can be the same as tod\n"
" bsize: number of samples in each block\n"
" moment: moment to compute, should be >= 1\n"
" central: whether to compute the central moment in each block\n"
Expand All @@ -1193,6 +1300,24 @@ PYBINDINGS("so3g")
"Compute the nth moment in blocks on a float32 array.\n"
"\n"
"See details in docstring for block_moment.\n");
bp::def("block_minmax", block_minmax<float>,
"block_minmax(tod, out, bsize, mode, shift)\n"
"\n"
"Compute the minimum, maximum, or peak to peak in blocks on a float32 array.\n"
"\n"
"Args:\n"
" tod: data array (float32) with shape (ndet, nsamp)\n"
" out: output array (float32) with shape (ndet, nsamp)\n"
" can be the same as tod\n"
" bsize: number of samples in each block\n"
" mode: if 0 compute the block minimum, if 1 the maximum, anything else will compute the peak to peak\n"
" shift: sample to start block at, used in each row\n");
bp::def("block_minmax64", block_minmax<double>,
"block_minmax64(tod, out, bsize, mode, shift)\n"
"\n"
"Compute the minimum, maximum, or peak to peak in blocks on a float64 array.\n"
"\n"
"See details in docstring for block_minmax.\n");
bp::def("matched_jumps", matched_jumps<float>,
"matched_jumps(tod, out, min_size, bsize)\n"
"\n"
Expand Down Expand Up @@ -1229,6 +1354,28 @@ PYBINDINGS("so3g")
"Output will be 0 where jumps are not found and the assumed jump height where jumps are found.\n"
"\n"
"See details in docstring for find_quantized_jumps.\n");
bp::def("subtract_jump_heights", subtract_jump_heights<float>,
"subtract_jump_heights(tod, out, heights, jumps)"
"\n"
"For each detector, compute the cumulative effect of the jumps identified by the array 'heights' and the RangesMatrix 'jumps'."
"For each range in 'jumps', the values from 'heights' are checked and the size of the jump is either the largest positive"
"or the largest negative number (whichever has the largest absolute value)."
"The 'output' value is the difference of 'tod' and the accumulated jump vector."
"\n"
"Args:\n"
" tod: data array (float32) with shape (ndet, nsamp)\n"
" out: output array (float32) with shape (ndet, nsamp)\n"
" can be the same as tod\n"
" heights: the height of the jump at each samples\n"
" should be an array (float32) with shape (ndet, nsamp)\n"
" jumps: RangesMatrix with the jump locations and shape (ndet, nsamp).\n");
bp::def("subtract_jump_heights64", subtract_jump_heights<double>,
"subtract_jump_heights64(tod, out, heights, jumps)"
"\n"
"Subtract fit jump heights from known jump locatations in a float64 array."
"If multiple samples in a jump have different heights, the largest height is used.\n"
"\n"
"See details in docstring for subtract_jump_heights.\n");
bp::def("clean_flag", clean_flag,
"clean_flag(flag, width)"
"\n"
Expand Down
Loading

0 comments on commit b7b4e09

Please sign in to comment.