From 2b3edef59617607a2f61489e3610ae29fe962e41 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 8 Aug 2024 14:44:48 +0200
Subject: [PATCH 01/18] Ensure the distances to computed linear-interpolation
 weights from are stored as arrays.

Without this, NumPy implicitly generates arrays but makes them ragged (dtype=object),
which is bad for performance and disabled in newer versions.
---
 odl/discr/discr_utils.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 4d8d855d60f..50498fa1585 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -706,6 +706,8 @@ def _compute_nearest_weights_edge(idcs, ndist):
 
 def _compute_linear_weights_edge(idcs, ndist):
     """Helper for linear interpolation."""
+    ndist = np.asarray(ndist)
+
     # Get out-of-bounds indices from the norm_distances. Negative
     # means "too low", larger than or equal to 1 means "too high"
     lo = np.where(ndist < 0)

From 80b1d8e8b8873bafc88eb38b3c940eafde01762f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 8 Aug 2024 15:32:11 +0200
Subject: [PATCH 02/18] Failure to convert to array implies it is not a
 suitable input array.

In older NumPy, this would silently create an array-of-object, but that
has a wrong shape too.
---
 odl/util/vectorization.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/odl/util/vectorization.py b/odl/util/vectorization.py
index 6d668a8c0fc..460fa62e81c 100644
--- a/odl/util/vectorization.py
+++ b/odl/util/vectorization.py
@@ -21,7 +21,10 @@
 
 def is_valid_input_array(x, ndim=None):
     """Test if ``x`` is a correctly shaped point array in R^d."""
-    x = np.asarray(x)
+    try:
+        x = np.asarray(x)
+    except ValueError:
+        return False
 
     if ndim is None or ndim == 1:
         return x.ndim == 1 and x.size > 1 or x.ndim == 2 and x.shape[0] == 1

From 30626fe24e55a9a437255e10756d9813dc210f82 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 8 Aug 2024 15:36:51 +0200
Subject: [PATCH 03/18] Disable a test that currently fails and is probably not
 worth supporting.

Giving an inhomogeneous-arrays list as the result to a function may sometimes
be convenient, but the automatic conversion is both prone to hiding bugs and
always bad for performance, so better is to require that the result actually
has the correct shape or at least one that NumPy can directly broadcast to the
correct one.
---
 odl/test/discr/discr_utils_test.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/odl/test/discr/discr_utils_test.py b/odl/test/discr/discr_utils_test.py
index 464af4228e9..4d019156a6f 100644
--- a/odl/test/discr/discr_utils_test.py
+++ b/odl/test/discr/discr_utils_test.py
@@ -355,7 +355,8 @@ def func_tens_dual(x, out=None):
 
 
 func_tens_params = [(func_tens_ref, f)
-                    for f in [func_tens_oop, func_tens_oop_seq,
+                    for f in [ # func_tens_oop,
+                              func_tens_oop_seq,
                               func_tens_ip, func_tens_ip_seq]]
 func_tens = simple_fixture('func_tens', func_tens_params,
                            fmt=' {name} = {value[1].__name__} ')

From 380ee09d38bea4e4fb490a08cbbd1d58a87a7d09 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 8 Aug 2024 15:46:01 +0200
Subject: [PATCH 04/18] Ensure interpolation weight is computed with the
 correct array type.

---
 odl/discr/discr_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 50498fa1585..1c0dade8e19 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -801,7 +801,7 @@ def _evaluate(self, indices, norm_distances, out=None):
         # axis, resulting in a loop of length 2**ndim
         for lo_hi, edge in zip(product(*([['l', 'h']] * len(indices))),
                                product(*edge_indices)):
-            weight = 1.0
+            weight = np.array([1.0], dtype = self.values.dtype)
             # TODO(kohr-h): determine best summation order from array strides
             for lh, w_lo, w_hi in zip(lo_hi, low_weights, high_weights):
 

From c90044a9be5a3d76555491573ef65bb87d3acff2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 8 Aug 2024 16:44:36 +0200
Subject: [PATCH 05/18] Ensure mesh coordinate calculations are not carried out
 with dtype=object.

This would previously happen because meshgrids are passed in as ragged arrays,
and NumPy does not convert the rows to float dtype as should happen.
---
 odl/discr/discr_utils.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 1c0dade8e19..4d5e731d747 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -621,6 +621,8 @@ def _find_indices(self, x):
 
         # iterate through dimensions
         for xi, cvec in zip(x, self.coord_vecs):
+            xi = np.asarray(xi, dtype = self.values.dtype)
+
             idcs = np.searchsorted(cvec, xi) - 1
 
             idcs[idcs < 0] = 0

From bd73b42b646cdadbb0989ffd633564e276a0cf49 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 8 Aug 2024 17:21:30 +0200
Subject: [PATCH 06/18] Linter whitespace rules.

---
 odl/discr/discr_utils.py           | 4 ++--
 odl/test/discr/discr_utils_test.py | 6 +++---
 2 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 4d5e731d747..e75a3c5abff 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -621,7 +621,7 @@ def _find_indices(self, x):
 
         # iterate through dimensions
         for xi, cvec in zip(x, self.coord_vecs):
-            xi = np.asarray(xi, dtype = self.values.dtype)
+            xi = np.asarray(xi, dtype=self.values.dtype)
 
             idcs = np.searchsorted(cvec, xi) - 1
 
@@ -803,7 +803,7 @@ def _evaluate(self, indices, norm_distances, out=None):
         # axis, resulting in a loop of length 2**ndim
         for lo_hi, edge in zip(product(*([['l', 'h']] * len(indices))),
                                product(*edge_indices)):
-            weight = np.array([1.0], dtype = self.values.dtype)
+            weight = np.array([1.0], dtype=self.values.dtype)
             # TODO(kohr-h): determine best summation order from array strides
             for lh, w_lo, w_hi in zip(lo_hi, low_weights, high_weights):
 
diff --git a/odl/test/discr/discr_utils_test.py b/odl/test/discr/discr_utils_test.py
index 4d019156a6f..02019a248fc 100644
--- a/odl/test/discr/discr_utils_test.py
+++ b/odl/test/discr/discr_utils_test.py
@@ -355,9 +355,9 @@ def func_tens_dual(x, out=None):
 
 
 func_tens_params = [(func_tens_ref, f)
-                    for f in [ # func_tens_oop,
-                              func_tens_oop_seq,
-                              func_tens_ip, func_tens_ip_seq]]
+                    for f in [  # func_tens_oop,
+                             func_tens_oop_seq,
+                             func_tens_ip, func_tens_ip_seq]]
 func_tens = simple_fixture('func_tens', func_tens_params,
                            fmt=' {name} = {value[1].__name__} ')
 

From 600fecd89eae318dc414184a4573751db6abbb07 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Fri, 9 Aug 2024 16:15:31 +0200
Subject: [PATCH 07/18] For integral data, interpolating on integral points may
 not be appropriate.

One of the tests samples from an integral grid at non-integral points.
This failed after the explicit conversion introduced in c90044a.
Falling back to `float` fixes the test case, though perhaps it would
be better to select a dedicated meshing dtype.
---
 odl/discr/discr_utils.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index e75a3c5abff..2a3b77d3444 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -621,7 +621,10 @@ def _find_indices(self, x):
 
         # iterate through dimensions
         for xi, cvec in zip(x, self.coord_vecs):
-            xi = np.asarray(xi, dtype=self.values.dtype)
+            try:
+                xi = np.asarray(xi).astype(self.values.dtype, casting='safe')
+            except TypeError:
+                xi = np.asarray(xi, dtype=float)
 
             idcs = np.searchsorted(cvec, xi) - 1
 

From 3e3ea87003cdc457bc8b535843a02343baf9a5a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Mon, 12 Aug 2024 16:44:19 +0200
Subject: [PATCH 08/18] NumPy changed what exception is raised for ufunc-call
 with wrong number of arguments.

The ODL tests checked against the previously raised `ValueError`, but in newer
versions `TypeError` is raised instead, which these tests could not handle.
---
 odl/test/discr/discr_space_test.py | 6 ++++--
 odl/test/space/tensors_test.py     | 6 ++++--
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/odl/test/discr/discr_space_test.py b/odl/test/discr/discr_space_test.py
index 4b254e5d494..366b497d56e 100644
--- a/odl/test/discr/discr_space_test.py
+++ b/odl/test/discr/discr_space_test.py
@@ -897,7 +897,9 @@ def test_ufunc_corner_cases(odl_tspace_impl):
 
     # --- UFuncs with nin = 1, nout = 1 --- #
 
-    with pytest.raises(ValueError):
+    wrong_argcount_error = ValueError if np.__version__<"1.21" else TypeError
+
+    with pytest.raises(wrong_argcount_error):
         # Too many arguments
         x.__array_ufunc__(np.sin, '__call__', x, np.ones((2, 3)))
 
@@ -928,7 +930,7 @@ def test_ufunc_corner_cases(odl_tspace_impl):
 
     # --- UFuncs with nin = 2, nout = 1 --- #
 
-    with pytest.raises(ValueError):
+    with pytest.raises(wrong_argcount_error):
         # Too few arguments
         x.__array_ufunc__(np.add, '__call__', x)
 
diff --git a/odl/test/space/tensors_test.py b/odl/test/space/tensors_test.py
index 1af171b20e3..ef8d99e8825 100644
--- a/odl/test/space/tensors_test.py
+++ b/odl/test/space/tensors_test.py
@@ -1530,7 +1530,9 @@ def test_ufunc_corner_cases(odl_tspace_impl):
 
     # --- Ufuncs with nin = 1, nout = 1 --- #
 
-    with pytest.raises(ValueError):
+    wrong_argcount_error = ValueError if np.__version__<"1.21" else TypeError
+
+    with pytest.raises(wrong_argcount_error):
         # Too many arguments
         x.__array_ufunc__(np.sin, '__call__', x, np.ones((2, 3)))
 
@@ -1561,7 +1563,7 @@ def test_ufunc_corner_cases(odl_tspace_impl):
 
     # --- Ufuncs with nin = 2, nout = 1 --- #
 
-    with pytest.raises(ValueError):
+    with pytest.raises(wrong_argcount_error):
         # Too few arguments
         x.__array_ufunc__(np.add, '__call__', x)
 

From 16474de5e3fec8cd98705094aa6d8cb075fe613a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Tue, 13 Aug 2024 16:35:49 +0200
Subject: [PATCH 09/18] Don't rely on obsolete NumPy inhomogenous arrays for
 boundary specification.

The partition builders are quite liberal (arguably, unnecessarily) in the format of the
boolean lists specifying along which dimensions the boundary should be included in
the grid and along which not. Previously, these free-form lists were passed through
NumPy as a first step of analysing their form, but NumPy itself is now more strict
in what can be in an array (unless explicitly asked with `dtype=object`, in which
case however a plain list could be used just as well).
These checks can also be done in a more direct fashion, presuming that the caller
actually follows the specification.
---
 odl/discr/grid.py     |  4 ++--
 odl/util/normalize.py | 10 +++++-----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/odl/discr/grid.py b/odl/discr/grid.py
index ef0d9bd2c35..0317629fe93 100644
--- a/odl/discr/grid.py
+++ b/odl/discr/grid.py
@@ -1111,8 +1111,8 @@ def uniform_grid_fromintv(intv_prod, shape, nodes_on_bdry=True):
 
     shape = normalized_scalar_param_list(shape, intv_prod.ndim, safe_int_conv)
 
-    if np.shape(nodes_on_bdry) == ():
-        nodes_on_bdry = ([(bool(nodes_on_bdry), bool(nodes_on_bdry))] *
+    if isinstance(nodes_on_bdry, bool):
+        nodes_on_bdry = ([(nodes_on_bdry, nodes_on_bdry)] *
                          intv_prod.ndim)
     elif intv_prod.ndim == 1 and len(nodes_on_bdry) == 2:
         nodes_on_bdry = [nodes_on_bdry]
diff --git a/odl/util/normalize.py b/odl/util/normalize.py
index e28911d32bf..434e6456b2d 100644
--- a/odl/util/normalize.py
+++ b/odl/util/normalize.py
@@ -278,11 +278,11 @@ def normalized_nodes_on_bdry(nodes_on_bdry, length):
     >>> normalized_nodes_on_bdry([[True, False], False, True], length=3)
     [(True, False), (False, False), (True, True)]
     """
-    shape = np.shape(nodes_on_bdry)
-    if shape == ():
-        out_list = [(bool(nodes_on_bdry), bool(nodes_on_bdry))] * length
-    elif length == 1 and shape == (2,):
-        out_list = [(bool(nodes_on_bdry[0]), bool(nodes_on_bdry[1]))]
+    if isinstance(nodes_on_bdry, bool):
+        return [(bool(nodes_on_bdry), bool(nodes_on_bdry))] * length
+    elif (length == 1 and len(nodes_on_bdry) == 2
+          and all(isinstance(d, bool) for d in nodes_on_bdry)):
+        return [nodes_on_bdry[0], nodes_on_bdry[1]]
     elif len(nodes_on_bdry) == length:
         out_list = []
 

From b6a912a87d082f20ad76b8316fb41a55b31e9e63 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Tue, 13 Aug 2024 18:10:28 +0200
Subject: [PATCH 10/18] Explicitly go through elements that may need to be
 broadcasted in fn. eval for points_collocation.

Not doing this lead to some of the familiar implicit dtype=object fallbacks that
are illegal in modern NumPy.
---
 odl/discr/discr_utils.py | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 2a3b77d3444..900e422241e 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -1339,12 +1339,12 @@ def dual_use_func(x, out=None, **kwargs):
             elif tensor_valued:
                 # The out object can be any array-like of objects with shapes
                 # that should all be broadcastable to scalar_out_shape.
-                results = np.array(out)
-                if results.dtype == object or scalar_in:
+
+                if any(res.shape != scalar_out_shape for res in out) or scalar_in:
                     # Some results don't have correct shape, need to
                     # broadcast
                     bcast_res = []
-                    for res in results.ravel():
+                    for res in out:
                         if ndim == 1:
                             # As usual, 1d is tedious to deal with. This
                             # code deals with extra dimensions in result
@@ -1355,17 +1355,17 @@ def dual_use_func(x, out=None, **kwargs):
                             if shp and shp[0] == 1:
                                 res = res.reshape(res.shape[1:])
                         bcast_res.append(
-                            np.broadcast_to(res, scalar_out_shape))
+                            np.broadcast_to(res, scalar_out_shape).astype(scalar_out_dtype))
 
                     out_arr = np.array(bcast_res, dtype=scalar_out_dtype)
-                elif results.dtype != scalar_out_dtype:
+                else:
+                    out_arr = np.asarray(out)
+                if out_arr.dtype != scalar_out_dtype:
                     raise ValueError(
                         'result is of dtype {}, expected {}'
-                        ''.format(dtype_repr(results.dtype),
+                        ''.format(dtype_repr(out_arr.dtype),
                                   dtype_repr(scalar_out_dtype))
                     )
-                else:
-                    out_arr = results
 
                 out = out_arr.reshape(out_shape)
 

From 2828057f36eb9962a39bd8f9eea3695215598a38 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Wed, 14 Aug 2024 12:49:53 +0200
Subject: [PATCH 11/18] Coding style details criticised by pep8speaks.

---
 odl/discr/discr_utils.py           | 3 ++-
 odl/test/discr/discr_space_test.py | 2 +-
 odl/test/space/tensors_test.py     | 2 +-
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 900e422241e..05154101f71 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -1355,7 +1355,8 @@ def dual_use_func(x, out=None, **kwargs):
                             if shp and shp[0] == 1:
                                 res = res.reshape(res.shape[1:])
                         bcast_res.append(
-                            np.broadcast_to(res, scalar_out_shape).astype(scalar_out_dtype))
+                            np.broadcast_to(res, scalar_out_shape)
+                              .astype(scalar_out_dtype))
 
                     out_arr = np.array(bcast_res, dtype=scalar_out_dtype)
                 else:
diff --git a/odl/test/discr/discr_space_test.py b/odl/test/discr/discr_space_test.py
index 366b497d56e..ebfda76ae1e 100644
--- a/odl/test/discr/discr_space_test.py
+++ b/odl/test/discr/discr_space_test.py
@@ -897,7 +897,7 @@ def test_ufunc_corner_cases(odl_tspace_impl):
 
     # --- UFuncs with nin = 1, nout = 1 --- #
 
-    wrong_argcount_error = ValueError if np.__version__<"1.21" else TypeError
+    wrong_argcount_error = ValueError if np.__version__ < "1.21" else TypeError
 
     with pytest.raises(wrong_argcount_error):
         # Too many arguments
diff --git a/odl/test/space/tensors_test.py b/odl/test/space/tensors_test.py
index ef8d99e8825..422a5450fa2 100644
--- a/odl/test/space/tensors_test.py
+++ b/odl/test/space/tensors_test.py
@@ -1530,7 +1530,7 @@ def test_ufunc_corner_cases(odl_tspace_impl):
 
     # --- Ufuncs with nin = 1, nout = 1 --- #
 
-    wrong_argcount_error = ValueError if np.__version__<"1.21" else TypeError
+    wrong_argcount_error = ValueError if np.__version__ < "1.21" else TypeError
 
     with pytest.raises(wrong_argcount_error):
         # Too many arguments

From 53c77602eeececf8b17532c45f166bb9dae0b43f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Tue, 27 Aug 2024 17:55:33 +0200
Subject: [PATCH 12/18] Add warning when falling back to `float` for
 interpolation coefficients.

When this happens it is likely that the used did something like starting
from an integer mesh, but in that case linear interpolation does not seem
very appropriate.
---
 odl/discr/discr_utils.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 05154101f71..9276e3725e9 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -20,6 +20,7 @@
 from builtins import object
 from functools import partial
 from itertools import product
+from warnings import warn
 
 import numpy as np
 
@@ -624,6 +625,8 @@ def _find_indices(self, x):
             try:
                 xi = np.asarray(xi).astype(self.values.dtype, casting='safe')
             except TypeError:
+                warn("Unable to infer accurate dtype for"
+                  +" interpolation coefficients, defaulting to `float`.")
                 xi = np.asarray(xi, dtype=float)
 
             idcs = np.searchsorted(cvec, xi) - 1

From 77e17f3ac311a42312ed8e5777fed08bc3a7c2cf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 29 Aug 2024 16:27:11 +0200
Subject: [PATCH 13/18] Manual broadcasting for the general
 tensor-with-inconsistent-dimensions case.

Previously, this was more or less automatically done in NumPy, but not anymore.
Normally, this is more likely to be user mistake. However, ODL actually
relies somewhat on such broadcasts when defining functions on "sparse"
mesh grids, so I added this functionality back be recursively
transversing lists of different-shape arrays (like in the old version
NumPy did, manually generating ragged dtype=object arrays).
---
 odl/discr/discr_utils.py | 51 ++++++++++++++++++++--------------------
 1 file changed, 25 insertions(+), 26 deletions(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 9276e3725e9..d0f4f26e832 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -933,6 +933,24 @@ def _func_out_type(func):
 
     return has_out, out_optional
 
+def _broadcast_nested_list(arr_lists, element_shape):
+    """ A generalisation of `np.broadcast_to`, applied to an arbitrarily
+    deep list (or tuple) eventually containing arrays or scalars. """
+    ndim = len(element_shape)
+    if isinstance(arr_lists, np.ndarray) or np.isscalar(arr_lists):
+        if ndim == 1:
+            # As usual, 1d is tedious to deal with. This
+            # code deals with extra dimensions in result
+            # components that stem from using x instead of
+            # x[0] in a function.
+            # Without this, broadcasting fails.
+            shp = getattr(arr_lists, 'shape', ())
+            if shp and shp[0] == 1:
+                arr_lists = arr_lists.reshape(arr_lists.shape[1:])
+        return np.broadcast_to(arr_lists, element_shape)
+    else:
+        return [_broadcast_nested_list(row, element_shape) for row in arr_lists]
+
 
 def sampling_function(func_or_arr, domain, out_dtype=None):
     """Return a function that can be used for sampling.
@@ -997,10 +1015,11 @@ def _default_oop(func_ip, x, **kwargs):
 
     def _default_ip(func_oop, x, out, **kwargs):
         """Default in-place variant of an out-of-place-only function."""
-        result = np.array(func_oop(x, **kwargs), copy=False)
-        if result.dtype == object:
+        result = func_oop(x, **kwargs)
+        try:
+            result = np.array(result, copy=False)
+        except ValueError:
             # Different shapes encountered, need to broadcast
-            flat_results = result.ravel()
             if is_valid_input_array(x, domain.ndim):
                 scalar_out_shape = out_shape_from_array(x)
             elif is_valid_input_meshgrid(x, domain.ndim):
@@ -1008,8 +1027,8 @@ def _default_ip(func_oop, x, out, **kwargs):
             else:
                 raise TypeError('invalid input `x`')
 
-            bcast_results = [np.broadcast_to(res, scalar_out_shape)
-                             for res in flat_results]
+            bcast_results = _broadcast_nested_list(result, scalar_out_shape)
+
             # New array that is flat in the `out_shape` axes, reshape it
             # to the final `out_shape + scalar_shape`, using the same
             # order ('C') as the initial `result.ravel()`.
@@ -1342,28 +1361,8 @@ def dual_use_func(x, out=None, **kwargs):
             elif tensor_valued:
                 # The out object can be any array-like of objects with shapes
                 # that should all be broadcastable to scalar_out_shape.
+                out_arr = np.asarray(_broadcast_nested_list(out, scalar_out_shape))
 
-                if any(res.shape != scalar_out_shape for res in out) or scalar_in:
-                    # Some results don't have correct shape, need to
-                    # broadcast
-                    bcast_res = []
-                    for res in out:
-                        if ndim == 1:
-                            # As usual, 1d is tedious to deal with. This
-                            # code deals with extra dimensions in result
-                            # components that stem from using x instead of
-                            # x[0] in a function.
-                            # Without this, broadcasting fails.
-                            shp = getattr(res, 'shape', ())
-                            if shp and shp[0] == 1:
-                                res = res.reshape(res.shape[1:])
-                        bcast_res.append(
-                            np.broadcast_to(res, scalar_out_shape)
-                              .astype(scalar_out_dtype))
-
-                    out_arr = np.array(bcast_res, dtype=scalar_out_dtype)
-                else:
-                    out_arr = np.asarray(out)
                 if out_arr.dtype != scalar_out_dtype:
                     raise ValueError(
                         'result is of dtype {}, expected {}'

From e80c43dffcfa09a4a1ad1443bc19b6c9347996c9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 29 Aug 2024 16:41:01 +0200
Subject: [PATCH 14/18] Revert "Disable a test that currently fails and is
 probably not worth supporting."

This reverts commit 30626fe24e55a9a437255e10756d9813dc210f82.

The test in question now succeeds after 77e17f3.
---
 odl/test/discr/discr_utils_test.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/odl/test/discr/discr_utils_test.py b/odl/test/discr/discr_utils_test.py
index 02019a248fc..464af4228e9 100644
--- a/odl/test/discr/discr_utils_test.py
+++ b/odl/test/discr/discr_utils_test.py
@@ -355,9 +355,8 @@ def func_tens_dual(x, out=None):
 
 
 func_tens_params = [(func_tens_ref, f)
-                    for f in [  # func_tens_oop,
-                             func_tens_oop_seq,
-                             func_tens_ip, func_tens_ip_seq]]
+                    for f in [func_tens_oop, func_tens_oop_seq,
+                              func_tens_ip, func_tens_ip_seq]]
 func_tens = simple_fixture('func_tens', func_tens_params,
                            fmt=' {name} = {value[1].__name__} ')
 

From 3b38e7fb83ee6fa6773009c5de3299fdd78f7c31 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 29 Aug 2024 16:38:07 +0200
Subject: [PATCH 15/18] Change doc test to have floating zero.

There is a check that the result must have the correct dtype, even
if it does not need to have the correct shape.
---
 odl/discr/discr_utils.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index d0f4f26e832..7d20d8087b9 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -120,7 +120,7 @@ def point_collocation(func, points, out=None, **kwargs):
     >>> ys = [3, 4]
     >>> mesh = sparse_meshgrid(xs, ys)
     >>> def vec_valued(x):
-    ...     return (x[0] - 1, 0, x[0] + x[1])  # broadcasting
+    ...     return (x[0] - 1., 0., x[0] + x[1])  # broadcasting
     >>> # For a function with several output components, we must specify the
     >>> # shape explicitly in the `out_dtype` parameter
     >>> func1 = sampling_function(

From 7c5e31f3fae922694361417ebb5c39bd01c0f204 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 29 Aug 2024 16:31:11 +0200
Subject: [PATCH 16/18] Omit unnecessary type coercion.

Pointed out by Jevgenija in https://github.com/odlgroup/odl/pull/1649#pullrequestreview-2268796987.
---
 odl/util/normalize.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/odl/util/normalize.py b/odl/util/normalize.py
index 434e6456b2d..27984f0640a 100644
--- a/odl/util/normalize.py
+++ b/odl/util/normalize.py
@@ -279,7 +279,7 @@ def normalized_nodes_on_bdry(nodes_on_bdry, length):
     [(True, False), (False, False), (True, True)]
     """
     if isinstance(nodes_on_bdry, bool):
-        return [(bool(nodes_on_bdry), bool(nodes_on_bdry))] * length
+        return [(nodes_on_bdry, nodes_on_bdry)] * length
     elif (length == 1 and len(nodes_on_bdry) == 2
           and all(isinstance(d, bool) for d in nodes_on_bdry)):
         return [nodes_on_bdry[0], nodes_on_bdry[1]]

From 9b3f4f42e3e8916c135bf0a9891537e44f316569 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 29 Aug 2024 18:04:43 +0200
Subject: [PATCH 17/18] Only perform manual broadcasting in case NumPy fails to
 broadcast homogeneously.

This might give slightly better performance in the expected homogeneous case.
---
 odl/discr/discr_utils.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 7d20d8087b9..1c0ea06cbca 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -1361,7 +1361,10 @@ def dual_use_func(x, out=None, **kwargs):
             elif tensor_valued:
                 # The out object can be any array-like of objects with shapes
                 # that should all be broadcastable to scalar_out_shape.
-                out_arr = np.asarray(_broadcast_nested_list(out, scalar_out_shape))
+                try:
+                    out_arr = np.asarray(out)
+                except ValueError:
+                    out_arr = np.asarray(_broadcast_nested_list(out, scalar_out_shape))
 
                 if out_arr.dtype != scalar_out_dtype:
                     raise ValueError(

From de6bc52f8321cb9afa7fb4928b680fdc5a38a685 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Justus=20Sagem=C3=BCller?= <justussa@kth.se>
Date: Thu, 29 Aug 2024 18:08:52 +0200
Subject: [PATCH 18/18] Explicitly pass domain dimension as an argument to
 inhomogeneous-broadcasting utility.

---
 odl/discr/discr_utils.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/odl/discr/discr_utils.py b/odl/discr/discr_utils.py
index 1c0ea06cbca..a2b49d317a3 100644
--- a/odl/discr/discr_utils.py
+++ b/odl/discr/discr_utils.py
@@ -933,10 +933,9 @@ def _func_out_type(func):
 
     return has_out, out_optional
 
-def _broadcast_nested_list(arr_lists, element_shape):
+def _broadcast_nested_list(arr_lists, element_shape, ndim):
     """ A generalisation of `np.broadcast_to`, applied to an arbitrarily
     deep list (or tuple) eventually containing arrays or scalars. """
-    ndim = len(element_shape)
     if isinstance(arr_lists, np.ndarray) or np.isscalar(arr_lists):
         if ndim == 1:
             # As usual, 1d is tedious to deal with. This
@@ -949,7 +948,7 @@ def _broadcast_nested_list(arr_lists, element_shape):
                 arr_lists = arr_lists.reshape(arr_lists.shape[1:])
         return np.broadcast_to(arr_lists, element_shape)
     else:
-        return [_broadcast_nested_list(row, element_shape) for row in arr_lists]
+        return [_broadcast_nested_list(row, element_shape, ndim) for row in arr_lists]
 
 
 def sampling_function(func_or_arr, domain, out_dtype=None):
@@ -1027,7 +1026,7 @@ def _default_ip(func_oop, x, out, **kwargs):
             else:
                 raise TypeError('invalid input `x`')
 
-            bcast_results = _broadcast_nested_list(result, scalar_out_shape)
+            bcast_results = _broadcast_nested_list(result, scalar_out_shape, domain.ndim)
 
             # New array that is flat in the `out_shape` axes, reshape it
             # to the final `out_shape + scalar_shape`, using the same
@@ -1364,7 +1363,7 @@ def dual_use_func(x, out=None, **kwargs):
                 try:
                     out_arr = np.asarray(out)
                 except ValueError:
-                    out_arr = np.asarray(_broadcast_nested_list(out, scalar_out_shape))
+                    out_arr = np.asarray(_broadcast_nested_list(out, scalar_out_shape, ndim=ndim))
 
                 if out_arr.dtype != scalar_out_dtype:
                     raise ValueError(