Changes look good but how about using .item() instead of .flat[0]

Sure, looks cleaner indeed. Let me make the changes.

Done.

@nschloe Thanks! The new error messages need tests. Can you please add them too?

Instead of fiddling around with the individual optimization methods, let me try and edit the function wrappers; perhaps the the changes will be more concise.

I think the PR is in better shape now. I essentially moved the scalar conversion to into ScalarFunction, and by that could save a couple of checks here and there.

I've not had much time to spend on scipy over the last few months, and don't think I'll be able to look before the branch point.

I don't think my concern has been addressed: If NumPy is redefining what is acceptable to treat as a scalar, shouldn't it be up to the user-provided functions to conform to that definition?

A good argument against would be if scipy.optimize documents support for treating rank-1 arrays as scalars. In that case I would prefer a deprecation matching NumPy.

@peterbell10 I agree, but from what I recall scipy is rather conservative when it comes to making changes to the user interface. That's why I "fixed" some scipy functions to only hand shape-0 arrays to numpy. To illustrate, an example would be goal functionals like

def f(x):
     return np.array([x ** 2 - 1.0])

This works in scipy without warning now. It would have to be changed to

def f(x):
     return x ** 2 - 1.0

to avoid the warning. This PR makes sure both versions pass without warning.

Yes, I know. I would like to see evidence that array-scalar conversion is a documented part of the current API, not just a coincidence due to NumPy's behavior. If it's not documented, then I see no reason to continue supporting it with an ugly workaround.

Just to be sure I understand correctly: You'd like to know if SciPy actually documents that you can do

def f(x):
     return np.array([x ** 2 - 1.0])

and if it's not documented, it'd be okay to just emit a warning like numpy does?

Perhaps an argument in favor of suppressing the warning:
Currently (and unfortunately), the input points x to an optimization method are always converted to shape (1,) if they are scalar. (See this issue report of mine which has been met with disagreement from scipy devs.) This means that in simple codes like

from scipy import optimize
import numpy as np


def f(x):
    out = -np.exp(-(x - 0.7)**2)
    return out

result = optimize.minimize(f, x0=0.3)

x and out will have shape (1,). (result.x is also of shape (1,) which I think is a mistake, but anyway.) With the change in numpy, the above code will print a warning. To avoid, the user would have to add [0]

from scipy import optimize
import numpy as np


def f(x):
    out = -np.exp(-(x - 0.7)**2)
    return out[0]

result = optimize.minimize(f, x0=0.3)

Not sure if that's desirable.

and consolidating handling into fewer places would make any future deprecation easier if we go down that route.

Indeed, that's what this PR is does, too (-> _wrap_scalar_function and ScalarFunction).

I'll leave this open for a bit in case an optimize expert wants to comment.

@peterbell10 Perhaps now that 1.7.0 is released is a good time to merge?

Are you sure np.isscalar(x) is the right choice in this PR (vs np.ndim(x) == 0(? Its documentation indicates it is often the wrong choice.

@eric-wieser Good point. I've checked the use of isscalar in scipy and found that in almost all cases, one should use ndim == 0. I've opened a discussion here.

np.ndim would be terrible here because it converts scalars into an array internally. So it's much slower in the case we want to be fast. I also don't see why we would care about object scalars here; everything should be numeric. So what are the downsides to np.isscalar then?

So what are the downsides to np.isscalar then?
everything should be numeric

Many of the isscalar conditionals are taking user input (e.g., from goal functionals in optimization), so we have no control over whether something is a Python scalar or not. The downside of isscalar is that "object scalars" (np.array(3.14)) are not isscalars. That's although they behave (almost) exactly like them.

So it's much slower in the case we want to be fast

It's faster indeed:

In [4]: %timeit np.ndim(3.14)
819 ns ± 0.784 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [5]: %timeit np.isscalar(3.14)
174 ns ± 0.705 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

We're talking about saving a few hundred nanoseconds here though which is less than the product of an integer with an object scalar:

In [6]: a = np.array(1.0)
In [7]: %timeit 2 * a
1.1 µs ± 8.02 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

A stronger argument for np.isscalar imho is readability. np.isscalar(x) is clear, np.dim(x) == 0 can be quite confusing. This leads me to believe that the better fix might be in numpy, not client code.

By "object scalar" I mean re.compile or Exception which are the examples in numpy/numpy#11882. I'm fine with zero-d arrays not passing isscalar because they literally are not scalars...

I'm fine with zero-d arrays not passing isscalar because they literally are not scalars...

Hehe. – I would argue they are scalars because they small, taste, and look like scalars. You can't tell the difference. In numerical computation, they behave exactly the same. How the data is stored in memory is just an implementational detail.

Edit: The only difference I can think of right now are Python-list computations which fail for Python floats, 3.14 * [1.0] but succeed for np-scalars np.asarray(3.14) * [1.0] because the list is converted to an array first.

Scalars are immutable and 0d arrays are mutable. That's a pretty big difference.

At least this is the case for python int, float and NumPy scalars like np.float32.