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