[WIP] fix some dimensionality inconsistencies

scipy/integrate/tests/test_bvp.py

@@ -319,6 +319,7 @@ def test_compute_global_jac():
     J = J.toarray()
 
     def J_block(h, p):
+        p = p[0]
         return np.array([
             [h**2*p**2/12 - 1, -0.5*h, -h**2*p**2/12 + 1, -0.5*h],
             [0.5*h*p**2, h**2*p**2/12 - 1, 0.5*h*p**2, 1 - h**2*p**2/12]

scipy/io/_fortran.py

@@ -133,7 +133,7 @@ def _read_size(self, eof_ok=False):
         elif len(b) < n:
             raise FortranFormattingError(
                 "End of file in the middle of the record size")
-        return int(np.frombuffer(b, dtype=self._header_dtype, count=1))
+        return int(np.frombuffer(b, dtype=self._header_dtype, count=1)[0])
 
     def write_record(self, *items):
         """

scipy/io/matlab/mio4.py

@@ -286,10 +286,10 @@ def shape_from_header(self, hdr):
 
             # Read only the row and column counts
             self.mat_stream.seek(dt.itemsize * (dims[0] - 1), 1)
-            rows = np.ndarray(shape=(1,), dtype=dt,
+            rows = np.ndarray(shape=(), dtype=dt,
                               buffer=self.mat_stream.read(dt.itemsize))
             self.mat_stream.seek(dt.itemsize * (dims[0] - 1), 1)
-            cols = np.ndarray(shape=(1,), dtype=dt,
+            cols = np.ndarray(shape=(), dtype=dt,
                               buffer=self.mat_stream.read(dt.itemsize))
 
             shape = (int(rows), int(cols))

scipy/linalg/tests/test_basic.py

@@ -1555,7 +1555,7 @@ def test_scaling(self):
         # Pre/post LAPACK 3.5.0 gives the same result up to an offset
         # since in each case col norm is x1000 greater and
         # 1000 / 32 ~= 1 * 32 hence balanced with 2 ** 5.
-        assert_allclose(int(np.diff(np.log2(np.diag(y)))), 5)
+        assert_allclose(int(np.diff(np.log2(np.diag(y)))[0]), 5)
 
     def test_scaling_order(self):
         A = np.array([[1, 0, 1e-4], [1, 1, 1e-2], [1e4, 1e2, 1]])
@@ -1565,7 +1565,7 @@ def test_scaling_order(self):
     def test_separate(self):
         _, (y, z) = matrix_balance(np.array([[1000, 1], [1000, 0]]),
                                    separate=1)
-        assert_equal(int(np.diff(np.log2(y))), 5)
+        assert_equal(int(np.diff(np.log2(y))[0]), 5)
         assert_allclose(z, np.arange(2))
 
     def test_permutation(self):

scipy/optimize/_basinhopping.py

@@ -310,7 +310,9 @@ def accept_reject(self, energy_new, energy_old):
         If new is higher than old, there is a chance it will be accepted,
         less likely for larger differences.
         """
-        w = math.exp(min(0, -float(energy_new - energy_old) * self.beta))
+        # The energy difference can be 0. If self.beta is inf, this emits a warning and
+        # returns nan.
+        w = math.exp(min(0, -(energy_new - energy_old) * self.beta))
         rand = self.random_state.rand()
         return w >= rand
 

scipy/optimize/_differentialevolution.py

@@ -874,8 +874,11 @@ def _calculate_population_energies(self, population):
 
         parameters_pop = self._scale_parameters(population)
         try:
-            calc_energies = list(self._mapwrapper(self.func,
-                                                  parameters_pop[0:nfevs]))
+            # Instead of casting this into an array and reshaping it,
+            # there's probably a nicer way to handle this.
+            calc_energies = np.array(list(self._mapwrapper(self.func,
+                                                           parameters_pop[0:nfevs])))
+            calc_energies = calc_energies.reshape(-1)
             energies[0:nfevs] = calc_energies
         except (TypeError, ValueError):
             # wrong number of arguments for _mapwrapper

scipy/optimize/_dual_annealing.py

@@ -94,7 +94,7 @@ def visiting(self, x, step, temperature):
             # Changing only one coordinate at a time based on strategy
             # chain step
             x_visit = np.copy(x)
-            visit = self.visit_fn(temperature, 1)
+            visit = self.visit_fn(temperature, 1)[0]
             if visit > self.TAIL_LIMIT:
                 visit = self.TAIL_LIMIT * self.rand_state.random_sample()
             elif visit < -self.TAIL_LIMIT:

scipy/optimize/_linprog_ip.py

@@ -218,29 +218,22 @@ def _get_delta(A, b, c, x, y, z, tau, kappa, gamma, eta, sparse=False,
         # Reference [4] Eq. 8.6
         rhatp = eta(gamma) * r_P
         rhatd = eta(gamma) * r_D
-        rhatg = np.array(eta(gamma) * r_G).reshape((1,))
+        rhatg = eta(gamma) * r_G
 
         # Reference [4] Eq. 8.7
         rhatxs = gamma * mu - x * z
-        rhattk = np.array(gamma * mu - tau * kappa).reshape((1,))
+        rhattk = gamma * mu - tau * kappa
 
         if i == 1:
             if ip:  # if the correction is to get "initial point"
                 # Reference [4] Eq. 8.23
                 rhatxs = ((1 - alpha) * gamma * mu -
                           x * z - alpha**2 * d_x * d_z)
-                rhattk = np.array(
-                    (1 -
-                     alpha) *
-                    gamma *
-                    mu -
-                    tau *
-                    kappa -
-                    alpha**2 *
-                    d_tau *
-                    d_kappa).reshape(
-                    (1,
-                     ))
+                rhattk = (
+                    (1 - alpha) * gamma * mu
+                    - tau * kappa
+                    - alpha**2 * d_tau * d_kappa
+                )
             else:  # if the correction is for "predictor-corrector"
                 # Reference [4] Eq. 8.13
                 rhatxs -= d_x * d_z

scipy/optimize/_trustregion_constr/report.py

@@ -17,11 +17,11 @@ def print_header(cls):
 
     @classmethod
     def print_iteration(cls, *args):
+        args = list(args)
         # args[3] is obj func. It should really be a float. However,
         # trust-constr typically provides a length 1 array. We have to coerce
         # it to a float otherwise the string format doesn't work.
-        args = list(args)
-        args[3] = float(args[3])
+        args[3] = args[3][0]
 
         iteration_format = ["{{:{}}}".format(x) for x in cls.ITERATION_FORMATS]
         fmt = "|" + "|".join(iteration_format) + "|"

scipy/optimize/lbfgsb.py

@@ -34,6 +34,7 @@
 ## Modifications by Travis Oliphant and Enthought, Inc. for inclusion in SciPy
 
 from __future__ import division, print_function, absolute_import
+import warnings
 
 import numpy as np
 from numpy import array, asarray, float64, int32, zeros
@@ -343,6 +344,12 @@ def func_and_grad(x):
             # until the completion of the current minimization iteration.
             # Overwrite f and g:
             f, g = func_and_grad(x)
+            # f should really always be a scalar
+            f = np.asarray(f)
+            if f.ndim > 0:
+                warnings.warn("f should return a scalar, not a vector of shape {}.".format(f.shape),
+                              category=UserWarning)
+            f = f.reshape(())
         elif task_str.startswith(b'NEW_X'):
             # new iteration
             n_iterations += 1

scipy/optimize/slsqp.py

@@ -397,8 +397,10 @@ def cjac(x, *args):
 
             # Compute objective function
             fx = func(x)
+            # TODO enforce returning an actual scalar, not a np.array of shape [1, 1, 1]
+            # or something like that?
             try:
-                fx = float(np.asarray(fx))
+                fx = np.asarray(fx).reshape([])
             except (TypeError, ValueError):
                 raise ValueError("Objective function must return a scalar")
             # Compute the constraints
@@ -445,7 +447,7 @@ def cjac(x, *args):
         # Call SLSQP
         slsqp(m, meq, x, xl, xu, fx, c, g, a, acc, majiter, mode, w, jw,
               alpha, f0, gs, h1, h2, h3, h4, t, t0, tol,
-              iexact, incons, ireset, itermx, line, 
+              iexact, incons, ireset, itermx, line,
               n1, n2, n3)
 
         # call callback if major iteration has incremented

scipy/optimize/tests/test__basinhopping.py

@@ -5,6 +5,7 @@
 import copy
 
 from numpy.testing import assert_almost_equal, assert_equal, assert_
+from scipy._lib._numpy_compat import suppress_warnings
 from pytest import raises as assert_raises
 import numpy as np
 from numpy import cos, sin
@@ -15,7 +16,8 @@
 
 
 def func1d(x):
-    f = cos(14.5 * x - 0.3) + (x + 0.2) * x
+    # f is a scalar, df a vector
+    f = cos(14.5 * x[0] - 0.3) + (x[0] + 0.2) * x[0]
     df = np.array(-14.5 * sin(14.5 * x - 0.3) + 2. * x + 0.2)
     return f, df
 
@@ -304,8 +306,10 @@ def callback2(x, f, accepted):
     def test_monotonic_basin_hopping(self):
         # test 1d minimizations with gradient and T=0
         i = 0
-        res = basinhopping(func1d, self.x0[i], minimizer_kwargs=self.kwargs,
-                           niter=self.niter, disp=self.disp, T=0)
+        with suppress_warnings() as sup:
+            sup.filter(RuntimeWarning, "invalid value encountered in double_scalars")
+            res = basinhopping(func1d, self.x0[i], minimizer_kwargs=self.kwargs,
+                               niter=self.niter, disp=self.disp, T=0)
         assert_almost_equal(res.x, self.sol[i], self.tol)
 
 

scipy/optimize/tests/test__differential_evolution.py

@@ -258,6 +258,8 @@ def test_args_tuple_is_passed(self):
         def quadratic(x, *args):
             if type(args) != tuple:
                 raise ValueError('args should be a tuple')
+            # the energy is a scalar
+            x = x[0]
             return args[0] + args[1] * x + args[2] * x**2.
 
         result = differential_evolution(quadratic,

scipy/optimize/tests/test__numdiff.py

@@ -1,8 +1,8 @@
 from __future__ import division
 
-import math
 from itertools import product
 
+import math
 import numpy as np
 from numpy.testing import assert_allclose, assert_equal, assert_
 from pytest import raises as assert_raises
@@ -185,6 +185,8 @@ def fun_non_numpy(self, x):
         return math.exp(x)
 
     def jac_non_numpy(self, x):
+        # Input can be float x or numpy.array([x]). Cast it all into numpy.array(x).
+        x = np.asarray(x).reshape(())
         return math.exp(x)
 
     def test_scalar_scalar(self):

scipy/optimize/tests/test__shgo.py

@@ -56,7 +56,8 @@ class StructTest2(StructTestFunction):
     """
 
     def f(self, x):
-        return (x - 30) * numpy.sin(x)
+        # f must return a scalar
+        return (x[0] - 30) * numpy.sin(x[0])
 
     def g(x):
         return 58 - numpy.sum(x, axis=0)

scipy/optimize/tests/test_optimize.py

@@ -132,6 +132,7 @@ def test_cg(self):
 
     def test_cg_cornercase(self):
         def f(r):
+            r = r[0]
             return 2.5 * (1 - np.exp(-1.5*(r - 0.5)))**2
 
         # Check several initial guesses. (Too far away from the
@@ -177,8 +178,9 @@ def test_bfgs(self):
 
     def test_bfgs_infinite(self):
         # Test corner case where -Inf is the minimum.  See gh-2019.
-        func = lambda x: -np.e**-x
-        fprime = lambda x: -func(x)
+        func = lambda x: -np.e**-x[0]
+        fprime = lambda x: np.e**-x
+
         x0 = [0]
         olderr = np.seterr(over='ignore')
         try:
@@ -544,6 +546,7 @@ def test_bfgs_numerical_jacobian(self):
 
     def test_bfgs_gh_2169(self):
         def f(x):
+            x = x[0]
             if x < 0:
                 return 1.79769313e+308
             else:
@@ -554,7 +557,7 @@ def f(x):
     def test_bfgs_double_evaluations(self):
         # check bfgs does not evaluate twice in a row at same point
         def f(x):
-            xp = float(x)
+            xp = x[0]
             assert xp not in seen
             seen.add(xp)
             return 10*x**2, 20*x
@@ -851,6 +854,10 @@ def test_no_increase(self, method):
         # initial point.
 
         def func(x):
+            # x can be scalar or np.array([x]), output must return a scalar
+            # TODO there should be some consistency in x. np.array([x]) makes sense for
+            #      an optimization problem
+            x = np.asarray(x).reshape(())
             return (x - 1)**2
 
         def bad_grad(x):
@@ -1433,7 +1440,7 @@ def test_1D(self):
         def f(x):
             assert_(len(x.shape) == 1)
             assert_(x.shape[0] == 1)
-            return x ** 2
+            return x[0] ** 2
 
         optimize.brute(f, [(-1, 1)], Ns=3, finish=None)
 

scipy/signal/filter_design.py

@@ -4158,13 +4158,15 @@ def cheb2ap(N, rs):
 
 
 def _vratio(u, ineps, mp):
-    [s, c, d, phi] = special.ellipj(u, mp)
+    # u is the corner of a 1D simplex; _vratio must return a scalar.
+    x = u[0]
+    [s, c, d, phi] = special.ellipj(x, mp)
     ret = abs(ineps - s / c)
     return ret
 
 
 def _kratio(m, k_ratio):
-    m = float(m)
+    m = m[0]
     if m < 0:
         m = 0.0
     if m > 1:

scipy/sparse/tests/test_base.py

@@ -780,7 +780,7 @@ def check_setdiag(a, b, k):
             # correct length, and too short or too long vectors
             for r in [-1, len(np.diag(a, k)), 2, 30]:
                 if r < 0:
-                    v = int(np.random.randint(1, 20, size=1))
+                    v = np.random.randint(1, 20, size=1)[0]
                 else:
                     v = np.random.randint(1, 20, size=r)
 

scipy/stats/_distn_infrastructure.py

@@ -2414,7 +2414,9 @@ def fit_loc_scale(self, data, *args):
     def _entropy(self, *args):
         def integ(x):
             val = self._pdf(x, *args)
-            return entr(val)
+            ret = entr(val)
+            # cast return value into a scalar
+            return np.asarray(ret).reshape(())
 
         # upper limit is often inf, so suppress warnings when integrating
         _a, _b = self._get_support(*args)
@@ -2494,9 +2496,9 @@ def expect(self, func=None, args=(), loc=0, scale=1, lb=None, ub=None,
 
         Examples
         --------
-
-        To understand the effect of the bounds of integration consider
 
+        To understand the effect of the bounds of integration consider
+
         >>> from scipy.stats import expon
         >>> expon(1).expect(lambda x: 1, lb=0.0, ub=2.0)
         0.6321205588285578

scipy/stats/mstats_extras.py

@@ -171,8 +171,7 @@ def _hdsd_1D(data, prob):
             mx_ = np.fromiter([np.dot(w,xsorted[np.r_[list(range(0,k)),
                                                       list(range(k+1,n))].astype(int_)])
                                   for k in range(n)], dtype=float_)
-            mx_var = np.array(mx_.var(), copy=False, ndmin=1) * n / float(n-1)
-            hdsd[i] = float(n-1) * np.sqrt(np.diag(mx_var).diagonal() / float(n))
+            hdsd[i] = float(n-1) * np.sqrt(mx_.var() / float(n-1))
         return hdsd
 
     # Initialization & checks