Evaluating a function using a numpy array returns inf and nan

Question

I have the following array and numpy function:

import numpy import scipy.special i = numpy.linspace(0, 500, 501, dtype = int) def func(x, n): return scipy.special.eval_hermite(n, x) 

I am evaluating a function for each element in my numpy i array using two different evaluations:

Approach 1:

 hermites = func(0, i) 

Approach 2:

 hermites = [func(0, idx) for idx in i] 

Two approaches lead to two different answers. They are different, because after element 50 , approaching 1 starts to return inf and nan . Approach 2 also does not provide the correct value for each element i . However, he can calculate more of them. Approach 2 fails for i >= 64 .

Both approaches give me an answer at about the same time (0.7 s for len(i) = 15000 , determined using timeit ). What I do not understand are different results. This is because I learned to avoid for loops in python as much as possible. This time it is not.

The idea that this is connected with memory also crossed my mind. However, the assessment of one single element, i.e. print func(0, 64) also returns 0. (equal to the output of approach 2 ).

What's happening?

Answer 1

This is a bug in scipy, created occasionally by the amazing numpy "ufuncs" casting rules. The problem is that in scipy version 0.16 and later, when the first argument, eval_hermite is an integer array and the second is an integer scalar, the return value data type is a single-precision floating point ( numpy.float32 ). When the second argument is a 64-bit floating point value, the return type is numpy.float64 . The largest value that can be represented using float32 is much less than float64 , so when the second argument is an integer, eval_hermite will go to infinity much faster.

For example, this is with scipy 0.16.0 and numpy 1.10.1:

 In [26]: from scipy.special import eval_hermite 

Note that the data type of the return value is float32 :

 In [27]: eval_hermite([20, 50, 100, 200], 0) Out[27]: array([ 6.70442586e+11, -inf, inf, inf], dtype=float32) 

If the second argument is a floating point, the return type is float64 , and large values ​​can be represented:

 In [28]: eval_hermite([20, 50, 100, 200], 0.0) Out[28]: array([ 6.70442573e+011, -1.96078147e+039, 3.06851876e+093, 8.45055019e+216]) 

The job for your code is to always ensure that the second argument to eval_hermite is a floating point value. For example,

 hermites = func(0.0, i) 

This problem has been fixed in scipy development version (see https://github.com/scipy/scipy/pull/4896 ), so scipy 0.17 should not have this problem when it is released.