[Numpy-discussion] Consistency of numpy’s behavior with arrays of dimension 0

[Sylvain Corlay]

Hello everyone,

0d arrays are very convenient because they allow us to write functions that
are valid for both arrays and scalar-like argument, thanks to Boolean
indexing.

However, most special functions in numpy (and scipy) and most Boolean
operators, when called on 0d arrays, return scalars rather than 0d arrays,
which is quite annoying.

For example, numpy.exp called on a 0d array containing a float number
returns a float, rather than a 0d array, and if x is a 0d array, x > 0
returns a Boolean, rather than a 0d array containing a Boolean.

What I would expect is the following.

-If x is a 0d array containing float, I expect numpy.exp(x) to return a 0d
array, and x>0 to return a 0d array containing a Boolean.

-If x is a scalar, numpy.exp(x) returns the expected scalar type, and x>0
returns the expected Boolean.

Here is an example of a simple function that suffers from this issue (a
corrected version is proposed later)
import numpy
from scipy.stats import norm

def normal_time_value(sig, m, strikes):
    """
    The three arguments are array-like and have the same shape.
    Consider a random variable G ~ N (m , sig^2)
    The function returns          E[(G-K)+] - (E[G]-K)_+
    which is also equal to         E[(K-G)+] - (K-E[G])_+
    """
    sig = numpy.array(sig)
    strikes = numpy.array(strikes)
    m = numpy.array(m)
    tv = numpy.zeros(strikes.shape)
    tv[sig < 0] = numpy.nan         # sig must be nonnegative
    non0 = sig > 0.0
    dev = numpy.where(non0, (m - strikes) / sig, numpy.nan)
    tv[non0] = numpy.where(strikes > m, (m - strikes) * norm.cdf(dev) + sig
* norm.pdf(dev), (strikes - m) * norm.cdf(-dev) + sig * norm.pdf(dev))[non0]
    return tv

This function does not work with scalars or 0d arrays. To make it work, we
need to modify it in the following fashion: reconvert intermediate results
to 0d array to take advantage of the Boolean indexing.

import numpy
from scipy.stats import norm

def normal_time_value(sig, m, strikes):
    """
    The three arguments are array-like and have the same shape.
    Consider a random variable G ~ N (m , sig^2)
    The function returns          E[(G-K)+] - (E[G]-K)_+
    which is also equal to        E[(K-G)+] - (K-E[G])_+
    """
    sig = numpy.array(sig)
    strikes = numpy.array(strikes)
    m = numpy.array(m)
    tv = numpy.zeros(strikes.shape)
    tv[numpy.array(sig < 0)] = numpy.nan         # sig must be nonnegative
    non0 = numpy.array(sig > 0.0)
    dev = numpy.where(non0, (m - strikes) / sig, numpy.nan)
    tv[non0] = numpy.where(numpy.array(strikes > m), (m - strikes) *
norm.cdf(dev) + sig * norm.pdf(dev), (strikes - m) * norm.cdf(-dev) + sig *
norm.pdf(dev))[non0]
    return tv

This problem also affects functions like logical_and, logical_or and
logical_not, which all return numpy.bool_ type rather than 0d array of
dtype bool.
Travis Oliphant opened a ticket in track regarding the same issue some time
ago:  http://projects.scipy.org/numpy/ticket/823 It was transfered to
github https://github.com/numpy/numpy/issues/1421.

Best,

Sylvain
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/numpy-discussion/attachments/20130508/9c98d86b/attachment.html>