BlueBrain / voxcell

Tools to work with voxel based brain atlases.
Apache License 2.0
5 stars 3 forks source link

Error in out-of-bounds lookup with non-scalar VoxelData #19

Closed seirios closed 1 year ago

seirios commented 1 year ago

When doing a lookup in a VoxelData object where the per-voxel value has dimension 2 (e.g., a flat map), specifying outer_value leads to error because the shape of the result is assumed to be 1-dimensional.

E.g., setting outer_value = np.array([-2,-2]) leads to:

Traceback (most recent call last):
  File "<stdin>", line 5, in <module>
  File "/gpfs/bbp.cscs.ch/ssd/apps/bsd/2023-02-23/stage_applications/install_gcc-12.2.0-skylake/py-voxcell-3.1.4-52dics/lib/python3.10/site-packages/voxcell/voxel_data.py", line 179, in lookup
    result = np.full(voxel_idx.shape[:-1], outer_value)
  File "/gpfs/bbp.cscs.ch/ssd/apps/bsd/2023-02-23/stage_applications/install_gcc-12.2.0-skylake/py-numpy-1.22.4-vtwzpw/lib/python3.10/site-packages/numpy/core/numeric.py", line 344, in full
    multiarray.copyto(a, fill_value, casting='unsafe')
  File "<__array_function__ internals>", line 180, in copyto
ValueError: could not broadcast input array from shape (2,) into shape (2060,)

The issue is that the dimension of result in voxel_data.py:179 does not match the dimension of the values it should store (2 in this example).

A possible fix is the following:

item_shape = self.raw[0,0,0].shape  # can be () for scalars, or (2,) in our case
result = np.full(voxel_idx.shape[:-1] + item_shape, outer_value, dtype=self.raw.dtype)  # tuple addition does the right thing

This way we keep the shape of the array values, and also the original dtype (otherwise it takes it from outer_value).

eleftherioszisis commented 1 year ago

Hello @seirios, and thanks for opening this issue. Could you please provide a minimal example to reproduce this error? Thank you!

seirios commented 1 year ago

Sure @eleftherioszisis !

This works (assumes only three dimensions): vc.VoxelData(np.zeros((1,1,1)),(1,1,1)).lookup((-1,-1,-1),0)

But this does not (fourth dim is 1): vc.VoxelData(np.zeros((1,1,1,1)),(1,1,1)).lookup((-1,-1,-1),0)

Neither does this (fourth dim is 2, as in flat map): vc.VoxelData(np.zeros((1,1,1,2)),(1,1,1)).lookup((-1,-1,-1),(0,0))

edasubert commented 1 year ago

Hello @seirios , may I ask what version of python and pandas are you using?

seirios commented 1 year ago

Hi! I'm using unstable modules in BB5, so python 3.10.8 and pandas 1.5.1. Also, it's voxcell 3.1.4.

edasubert commented 1 year ago

Hello @seirios , please try the py-voxcell 3.1.5 so we can close this issue :slightly_smiling_face:

mgeplf commented 1 year ago

@seirios can you have a look at this to see if it fixes your issue?

mgeplf commented 1 year ago

ping @seirios?

seirios commented 1 year ago

Hi! Sorry for the big delay! I have tested it now, and it works as expected in isolation and in the context of my use case. Thanks!

mgeplf commented 11 months ago

Thanks for following up @seirios.