4

I am using the PLegendre function from the SHTOOLS package. It returns an array of Legendre polynomials for a particular argument. PLegendre(lmax,x) returns an array of Legendre polynomials P_0(x) to P_lmax(x). It works like this:

In [1]: from pyshtools import PLegendre
loading shtools documentation

In [2]: import numpy as np

In [3]: PLegendre(3,0.5)
Out[3]: array([ 1.    ,  0.5   , -0.125 , -0.4375])

I would like to pass an array as a parameter, so I use frompyfunc.

In [4]: legendre=np.frompyfunc(PLegendre,2,1)

In [5]: legendre(3,np.linspace(0,1,4))
Out[5]: 
array([array([ 1. ,  0. , -0.5, -0. ]),
   array([ 1.        ,  0.33333333, -0.33333333, -0.40740741]),
   array([ 1.        ,  0.66666667,  0.16666667, -0.25925926]),
   array([ 1.,  1.,  1.,  1.])], dtype=object)

The output is an array of arrays. I understand that I can create an array of elements from this by slicing the array.

In [6]: a=legendre(3,np.linspace(0,1,4))

In [7]: array([a[i][:] for i in xrange(4)])
Out[7]: 
array([[ 1.        ,  0.        , -0.5       , -0.        ],
   [ 1.        ,  0.33333333, -0.33333333, -0.40740741],
   [ 1.        ,  0.66666667,  0.16666667, -0.25925926],
   [ 1.        ,  1.        ,  1.        ,  1.        ]])

But.. is there a way to get to this directly, instead of having to slice the array of arrays?

Jishnu
  • 769
  • 8
  • 17

2 Answers2

3

I think it can not be done directly as already pointed out here in the case of np.vectorize which does almost the same thing. Note that your code is not faster then an ordinary for loop by using np.frompyfunc ... the code only looks nicer.

However, what you can do is using np.vstack instead of the list comprehension

a = legendre(3,np.linspace(0,1,4))
np.vstack(a)
Community
  • 1
  • 1
plonser
  • 3,323
  • 2
  • 18
  • 22
2

np.frompyfunc is compiled, so I'd have to dig into the source to see exactly what it is doing. But it appears to assume that the func output is an (inscrutable) Python object.

foo1 = np.frompyfunc(np.arange,1,1)
foo2 = np.vectorize(np.arange,otypes='O')

These 2 functions produce the same outputs, though foo1 is faster.

foo1(np.arange(4))

produces arrays with different sizes

array([array([], dtype=int32), array([0]), array([0, 1]), array([0, 1, 2])], dtype=object)

where as foo1(np.ones((4,)) are all the same, and in theory could be stacked.

There's no attempt, during the loop or after, to test whether the objects are arrays (or lists) and whether they can be combined into a single higher dimensional array.

plonser's use of vstack is a good idea. In fact frompyfunc plus vstack is faster than the more common list comprehension plus vstack.

In [54]: timeit np.vstack([np.arange(i) for i in 10*np.ones((10,))])
10000 loops, best of 3: 169 µs per loop

In [55]: timeit np.vstack(foo1(10*np.ones((10,))))
10000 loops, best of 3: 127 µs per loop
hpaulj
  • 221,503
  • 14
  • 230
  • 353