odlgroup / odl

Operator Discretization Library https://odlgroup.github.io/odl/
Mozilla Public License 2.0
375 stars 105 forks source link

NumPy dtype=object problem in linear resampling operator #1648

Closed leftaroundabout closed 2 months ago

leftaroundabout commented 3 months ago

The internally used array types, in particular norm_distances in ResamplingOperator are currently inconsistent. These should mostly be np.ndarray(dtype=float) or similar, but some of them are in fact lists or NumPy arrays with dtype=object.

In older versions of NumPy, this kind of discrepancy merely caused silent performance degradation. In newer ones, adding to an object-array is actually an error, which is caught by the test suite:

___________________________________________ [doctest] odl.discr.discr_ops.Resampling.__init__ _______________________________________________________
064 
065         Apply the corresponding resampling operator to an element:
066 
067         >>> print(resampling([0, 1, 0]))
068         [ 0.,  0.,  1.,  1.,  0.,  0.]
069 
070         With linear interpolation:
071 
072         >>> resampling = odl.Resampling(coarse_discr, fine_discr, 'linear')
073         >>> print(resampling([0, 1, 0]))
UNEXPECTED EXCEPTION: numpy.core._exceptions._UFuncOutputCastingError: Cannot cast ufunc 'add' output from dtype('O') to dtype('float64') with casting rule 'same_kind'
Traceback (most recent call last):
  File "/home/justussa/miniconda3/envs/odl_olddeps/lib/python3.7/doctest.py", line 1330, in __run
    compileflags, 1), test.globs)
  File "<doctest odl.discr.discr_ops.Resampling.__init__[8]>", line 1, in <module>
  File "/home/justussa/progwrit/python/odl/odl/operator/operator.py", line 694, in __call__
    out = self._call_out_of_place(x, **kwargs)
  File "/home/justussa/progwrit/python/odl/odl/discr/discr_ops.py", line 116, in _call
    interpolator, self.range.meshgrid, out=out_arr
  File "/home/justussa/progwrit/python/odl/odl/discr/discr_utils.py", line 172, in point_collocation
    out = func(points, **kwargs)
  File "/home/justussa/progwrit/python/odl/odl/discr/discr_utils.py", line 495, in per_axis_interp
    res = interpolator(x, out=out)
  File "/home/justussa/progwrit/python/odl/odl/discr/discr_utils.py", line 606, in __call__
    values = self._evaluate(indices, norm_distances, out)
  File "/home/justussa/progwrit/python/odl/odl/discr/discr_utils.py", line 815, in _evaluate
    out += np.asarray(self.values[edge]) * weight[vslice]
numpy.core._exceptions._UFuncOutputCastingError: Cannot cast ufunc 'add' output from dtype('O') to dtype('float64') with casting rule 'same_kind'
leftaroundabout commented 2 months ago

Fixed by #1655.