Fix out-of-bounds access in bland::compute_broadcast_shape()

[daniestevez]

The bland::compute_broadcast_shape() function is not implemented correctly and can perform out of bound accesses for some inputs, which can lead to aborts depending on the C++ library. In particular, I ran into this problem under Arch Linux when running bliss_find_hits with the single_coarse_guppi_59046_80036_DIAG_VOYAGER-1_0011.rawspec.0000.h5 example file. According to the backtrace, the crash happens during bliss:normalize(), which calls to bland::divide().

The problem is that compute_broadcast_shape() greedily consumes elements from shape by incrementing input_shape_index whenever the current components of shape and out_shape are equal or the current component of shape is 1. The following inputs all lead to out of bound accesses on shape:

• shape = [1], out_shape = [10, 1024]
• shape = [10], out_shape = [10, 1024]

I had Claude debugging this problem, and it said that Numpy-style broadcasting (which is most likely what is intended here) works by right aligning both shapes and implicitly padding with ones on the left when needed, it suggested the following simpler implementation, which simply obtains broadcast_shape by padding shape on the left with ones to make its length equal to out_shape.size().

However, some more thought is probably needed about what is the API contract for this broadcast_shape() function. That is, what inputs are acceptable. For instance, what Claude suggested does not treat correctly the case out_shape.size() < shape.size(), because then offset is negative and we also run into out of bounds accesses. Someone with more knowledge of bliss should take a closer look at this.

All I can say for now is that this patch allowed me to run bliss correctly with the Voyager test data, and with other data files, in my Arch Linux system.

[daniestevez]

For more context, this is the crash I was running into

$ HDF5_PLUGIN_PATH=/usr/lib/python3.14/site-packages/hdf5plugin/plugins/ ./bliss/bliss_find_hits /tmp/voyager/single_coarse_guppi_59046_80036_DIAG_VOYAGER-1_0011.rawspec.0000.h5 -d cpu
Found a usable device
INFO: HDF5 looking for filter plugins in:
    '/usr/lib/python3.14/site-packages/hdf5plugin/plugins'
INFO: Got source name VOYAGER-1
WARN: writing to file single_coarse_guppi_59046_80036_DIAG_VOYAGER-1_0011.rawspec.0000.capnp that already exists. Existing contents will be overwritten.
/usr/include/c++/16/bits/stl_vector.h:1253: std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](size_type) [with _Tp = long int; _Alloc = std::allocator<long int>; reference = long int&; size_type = long unsigned int]: Assertion '__n < this->size()' failed.
zsh: abort (core dumped)  HDF5_PLUGIN_PATH=/usr/lib/python3.14/site-packages/hdf5plugin/plugins/   -d

I was using -d cpu because when using the GPU I was also running into a CUDA out-of-memory error with my RTX 3070, but at some point I discovered that I could avoid this by narrowing down the Doppler drift range. I think that the same issue also happens when using the CPU, but I don't remember if I tested this.

This is the gdb backtrace

(gdb) bt full
#0  0x00007fffe6c9a09c in ?? () from /usr/lib/libc.so.6
No symbol table info available.
#1  0x00007fffe6c3e5d0 in raise () from /usr/lib/libc.so.6
No symbol table info available.
#2  0x00007fffe6c25685 in abort () from /usr/lib/libc.so.6
No symbol table info available.
#3  0x00007fffe709d5bd in std::__glibcxx_assert_fail(char const*, int, char const*, char const*) () from /usr/lib/libstdc++.so.6
No symbol table info available.
#4  0x00005555555e3f50 in std::vector<long, std::allocator<long> >::operator[](unsigned long) ()
No symbol table info available.
#5  0x0000555555c5eae9 in bland::compute_broadcast_shape(std::vector<long, std::allocator<long> >, std::vector<long, std::allocator<long> >) ()
No symbol table info available.
#6  0x00005555557814f8 in bland::ndarray bland::elementwise_binary_op<float, float, float, bland::elementwise_divide_op_ts>(bland::ndarray&, bland::ndarray const&, bland::ndarray const&) ()
No symbol table info available.
#7  0x00005555556c6c23 in bland::ndarray bland::elementwise_binary_op_impl_wrapper<bland::elementwise_divide_op_ts>::call<float, float, float>(bland::ndarray, bland::ndarray const&, bland::ndarray const&) ()
No symbol table info available.
#8  0x00005555556a9145 in bland::ndarray bland::dispatch<bland::elementwise_binary_op_impl_wrapper<bland::elementwise_divide_op_ts>, float, float>(bland::ndarray&, bland::ndarray const&, bland::ndarray const&) ()
No symbol table info available.
#9  0x0000555555699ee1 in bland::ndarray bland::dispatch<bland::elementwise_binary_op_impl_wrapper<bland::elementwise_divide_op_ts>, float>(bland::ndarray&, bland::ndarray const&, bland::ndarray const&) ()
No symbol table info available.
#10 0x0000555555697704 in bland::ndarray bland::dispatch<bland::elementwise_binary_op_impl_wrapper<bland::elementwise_divide_op_ts>>(bland::ndarray&, bland::ndarray const&, bland::ndarray const&) ()
No symbol table info available.
#11 0x0000555555696b1e in bland::ndarray bland::cpu::divide_cpu<bland::ndarray>(bland::ndarray, bland::ndarray) ()
No symbol table info available.
#12 0x0000555555695330 in bland::ndarray device_dispatch<bland::ndarray (*)(bland::ndarray, bland::ndarray), bland::ndarray (*)(bland::ndarray, bland::ndarray), bland::ndarray, bland::ndarray>(bland::ndarray (*)(bland::ndarray, bland::ndarray), bland::ndarray (*)(bland::ndarray, bland::ndarray), bland::ndarray, bland::ndarray) ()
No symbol table info available.
#13 0x0000555555694c76 in bland::ndarray bland::divide<bland::ndarray>(bland::ndarray, bland::ndarray) ()
No symbol table info available.
#14 0x00005555555f21eb in bliss::normalize(bliss::coarse_channel) ()
No symbol table info available.
#15 0x00005555555f23c2 in bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}::operator()(bliss::coarse_channel) const ()
No symbol table info available.
#16 0x00005555555f2c57 in bliss::coarse_channel std::__invoke_impl<bliss::coarse_channel, bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}&, bliss::coarse_channel>(std::__invoke_other, bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}&, bliss::coarse_channel&&) ()
No symbol table info available.
#17 0x00005555555f2ade in std::enable_if<is_invocable_r_v<bliss::coarse_channel, bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}&, bliss::coarse_channel>, bliss::coarse_channel>::type std::__invoke_r<bliss::coarse_channel, bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}&, bliss::coarse_channel>(bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}&, bliss::coarse_channel&&) ()
No symbol table info available.
#18 0x00005555555f297e in std::_Function_handler<bliss::coarse_channel (bliss::coarse_channel), bliss::normalize(bliss::scan)::{lambda(bliss::coarse_channel)#1}>::_M_invoke(std::_Any_data const&, bliss::coarse_channel&&) ()
No symbol table info available.
#19 0x00005555555e45b5 in std::function<bliss::coarse_channel (bliss::coarse_channel)>::operator()(bliss::coarse_channel) const ()
No symbol table info available.
#20 0x00005555555def81 in bliss::scan::read_coarse_channel(int) ()
No symbol table info available.
#21 0x0000555555605bc9 in bliss::detail::bliss_scan_to_capnp_scan(Scan::Builder&, bliss::scan) ()
No symbol table info available.
#22 0x0000555555602178 in bliss::write_scan_hits_to_capnp_file(bliss::scan, std::basic_string_view<char, std::char_traits<char> >) ()
No symbol table info available.
#23 0x0000555555601502 in bliss::write_scan_hits_to_file(bliss::scan, std::basic_string_view<char, std::char_traits<char> >, std::__cxx11::basic_string<char, std::char_traits<char>, std::allocator<char> >, double) ()
No symbol table info available.
#24 0x00005555555776e1 in main ()

I also tried building bliss in Debug mode, and by doing that it is possible to see which inputs cause the crash in compute_broadcast_shape(), but I didn't take note of that.