2

I'm encountering a scipy function that seems to return a numpy array no matter what's passed to it. In my application I need to be able to pass scalars and lists only, so the only "problem" is that when I pass a scalar to the function an array with one element is returned (when I would expect a scalar). Should I ignore this behaviour, or hack up the function to ensure that when a scalar is passed a scalar is returned?

Example code:

#! /usr/bin/env python

import scipy
import scipy.optimize
from numpy import cos

# This a some function we want to compute the inverse of
def f(x):
    y = x + 2*cos(x)
    return y

# Given y, this returns x such that f(x)=y
def f_inverse(y):

    # This will be zero if f(x)=y
    def minimize_this(x):
        return y-f(x)

    # A guess for the solution is required
    x_guess = y
    x_optimized = scipy.optimize.fsolve(minimize_this, x_guess) # THE PROBLEM COMES FROM HERE
    return x_optimized

# If I call f_inverse with a list, a numpy array is returned
print f_inverse([1.0, 2.0, 3.0])
print type( f_inverse([1.0, 2.0, 3.0]) )

# If I call f_inverse with a tuple, a numpy array is returned
print f_inverse((1.0, 2.0, 3.0))
print type( f_inverse((1.0, 2.0, 3.0)) )

# If I call f_inverse with a scalar, a numpy array is returned
print f_inverse(1.0)
print type( f_inverse(1.0) )

# This is the behaviour I expected (scalar passed, scalar returned).
# Adding [0] on the return value is a hackey solution (then thing would break if a list were actually passed).
print f_inverse(1.0)[0] # <- bad solution
print type( f_inverse(1.0)[0] )

On my system the output of this is:

[ 2.23872989  1.10914418  4.1187546 ]
<type 'numpy.ndarray'>
[ 2.23872989  1.10914418  4.1187546 ]
<type 'numpy.ndarray'>
[ 2.23872989]
<type 'numpy.ndarray'>
2.23872989209
<type 'numpy.float64'>

I'm using SciPy 0.10.1 and Python 2.7.3 provided by MacPorts.

SOLUTION

After reading the answers below I settled on the following solution. Replace the return line in f_inverse with:

if(type(y).__module__ == np.__name__):
    return x_optimized
else:
    return type(y)(x_optimized)

Here return type(y)(x_optimized) causes the return type to be the same as the type the function was called with. Unfortunately this does not work if y is a numpy type, so if(type(y).__module__ == np.__name__) is used to detect numpy types using the idea presented here and exclude them from the type conversion.

Community
  • 1
  • 1
Douglas B. Staple
  • 10,510
  • 8
  • 31
  • 58

4 Answers4

3

The first line of the implementation in scipy.optimize.fsolve is:

x0 = array(x0, ndmin=1)

This means that your scalar will be turned into a 1-element sequence, and your 1-element sequence will be essentially unchanged.

The fact that it seems to work is an implementation detail, and I would refactor your code to not allow sending a scalar into fsolve. I know this might seem to go against duck-typing, but the function asks for an ndarray for that argument, so you should respect the interface to be robust to changes in implementation. I don't, however, see any problem with conditionally using x_guess = array(y, ndmin=1) for converting scalars into an ndarray in your wrapper function and converting the result back to scalar when necessary.

Here is the relevant part of docstring of fsolve function:

def fsolve(func, x0, args=(), fprime=None, full_output=0,
           col_deriv=0, xtol=1.49012e-8, maxfev=0, band=None,
           epsfcn=0.0, factor=100, diag=None):
    """
    Find the roots of a function.

    Return the roots of the (non-linear) equations defined by
    ``func(x) = 0`` given a starting estimate.

    Parameters
    ----------
    func : callable f(x, *args)
        A function that takes at least one (possibly vector) argument.
    x0 : ndarray
        The starting estimate for the roots of ``func(x) = 0``.

    ----SNIP----

    Returns
    -------
    x : ndarray
        The solution (or the result of the last iteration for
        an unsuccessful call).

    ----SNIP----
wim
  • 338,267
  • 99
  • 616
  • 750
  • Quite a few SciPy functions accept "array-like" input and will convert that to arrays implicitly. +1 for the docstring, though. – Fred Foo Sep 24 '12 at 13:27
  • Thanks. I'm still not sure what to do with this though, because I don't know if it's safe to ignore the difference between a scalar and a numpy array with one element. – Douglas B. Staple Sep 24 '12 at 13:34
  • In this case I think it's not safe, so my answer is you should indeed hack up your `f_inverse` function if you _really_ want the interface of scalar-in scalar-out. – wim Sep 24 '12 at 13:39
  • Yes rereading your answer I see that's implied. I'll accept this answer if nothing "nicer" comes up. – Douglas B. Staple Sep 24 '12 at 13:43
2

Here's how you can convert Numpy arrays to lists and Numpy scalars to Python scalars:

>>> x = np.float32(42)
>>> type(x)
<type 'numpy.float32'>
>>> x.tolist()
42.0

In other words, the tolist method on np.ndarray handles scalars specially.

That still leaves you with single-element lists, but those are easy enough to handle in the usual way.

Fred Foo
  • 355,277
  • 75
  • 744
  • 836
1

I guess wims answer really already says it mostly, but maybe this makes the differences clearer.

The scalar returned by numpy should with array[0] should be (almost?) fully compatible to the standard python float:

a = np.ones(2, dtype=float)
isinstance(a[0], float) == True # even this is true.

For the most part already the 1 sized array is compatible to both a scalar and list, though for example it is a mutable object while the float is not:

a = np.ones(1, dtype=float)
import math
math.exp(a) # works
# it is not isinstance though
isinstance(a, float) == False
# The 1-sized array works sometimes more like number:
bool(np.zeros(1)) == bool(np.asscalar(np.zeros(1)))
# While lists would be always True if they have more then one element.
bool([0]) != bool(np.zeros(1))

# And being in place might create confusion:
a = np.ones(1); c = a; c += 3
b = 1.; c = b; c += 3
a != b

So if the user should not know about it, I think the first is fine the second it is dangerous.

You can also use np.asscalar(result) to convert a size 1 array (of any dimension) to the correct python scalar:

In [29]: type(np.asscalar(a[0])) Out[29]: float

If you want to make sure there are no surprises for a user who is not supposed to know about numpy, you will have to at least get the 0's element if a scalar was passed in. If the user should be numpy aware, just documentation is probably as good.

seberg
  • 8,785
  • 2
  • 31
  • 30
  • 1
    The `isinstance` check works because the OP is passing in Python floats: `isinstance(np.float64(1), float)` is true, but `isinstance(np.float32(1), float)` is false. – Fred Foo Sep 24 '12 at 14:33
  • True, but a python float is always double, if someone actually passes in `np.float32(1)` they should know to expect numpy objects back. – seberg Sep 24 '12 at 14:39
1

As @wim pointed out, fsolve transforms your scalar into a ndarray of shape (1,) and returns an array of shape (1,).

If you really want to get a scalar as output, you could try to put the following at the end of your function:

if solution.size == 1:
    return solution.item()
return solution

(The item method copies an element of an array and return a standard Python scalar)

Pierre GM
  • 19,809
  • 3
  • 56
  • 67