BUG: Fix Akima1DInterpolator returning nans, issue scipy#5683

scipy/interpolate/_monotone.py

@@ -297,11 +297,12 @@ def __init__(self, x, y, axis=0):
         f1 = dm[2:]
         f2 = dm[:-2]
         f12 = f1 + f2
-        # These are the indices where the the slope at breakpoint is defined:
-        id_ = np.nonzero(f12 > 1e-9 * np.max(f12))[0]
-        # set the slope at breakpoint
-        t[id_] = (f1[id_] * m[id_ + 1] + f2[id_] * m[id_ + 2]) / f12[id_]
-
+        # These are the mask of where the the slope at breakpoint is defined:
+        ind = np.nonzero(f12 > 1e-9 * np.max(f12))
+        x_ind, y_ind = ind[0], ind[1:]
+        # Set the slope at breakpoint
+        t[ind] = (f1[ind] * m[(x_ind + 1,) + y_ind] +
+                  f2[ind] * m[(x_ind + 2,) + y_ind]) / f12[ind]
         # calculate the higher order coefficients
         c = (3. * m[2:-2] - 2. * t[:-1] - t[1:]) / dx
         d = (t[:-1] + t[1:] - 2. * m[2:-2]) / dx ** 2

scipy/interpolate/tests/test_interpolate.py

@@ -667,6 +667,15 @@ def test_eval_3d(self):
         yi[:, 1, 1] = 4. * yi_
         assert_allclose(ak(xi), yi)
 
+    def test_degenerate_case_broadcasting(self):
+        # This test is for issue #5683.
+        x = np.array([0, 1, 2])
+        y = np.vstack((x, x**2)).T
+        ak = Akima1DInterpolator(x, y)
+        x_eval = np.array([0.5, 1.5])
+        y_eval = ak(x_eval)
+        assert_allclose(y_eval, np.vstack((x_eval, x_eval**2)).T)
+
     def test_extend(self):
         x = np.arange(0., 11.)
         y = np.array([0., 2., 1., 3., 2., 6., 5.5, 5.5, 2.7, 5.1, 3.])