0

I would like to get a matrix of values given two ndarray's from a ufunc, for example:

degs = numpy.array(range(5))
pnts = numpy.array([0.0, 0.1, 0.2])
values = scipy.special.eval_chebyt(degs, pnts)

The above code doesn't work (it gives a ValueError because it tries to broadcast two arrays and fails since they have different shapes: (5,) and (3,)); I would like to get a matrix of values with rows corresponding to degrees and columns to points at which polynomials are evaluated (or vice versa, it doesn't matter).

Currently my workaround is simply to use for-loop:

values = numpy.zeros((5,3))
for j in range(5):
    values[j] = scipy.special.eval_chebyt(j, pnts)

Is there a way to do that? In general, how would you let a ufunc know you want an n-dimensional array if you have n array_like arguments?

I know about numpy.vectorize, but that seems neither faster nor more elegant than just a simple for-loop (and I'm not even sure you can apply it to an existent ufunc).

UPDATE What about ufunc's that receive 3 or more parameters? trying outer method gives a ValueError: outer product only supported for binary functions. For example, scipy.special.eval_jacobi.

gsarret
  • 160
  • 2
  • 7

2 Answers2

3

What you need is exactly the outer method of ufuncs:

ufunc.outer(A, B, **kwargs)

  Apply the ufunc op to all pairs (a, b) with a in A and b in B.

values = scipy.special.eval_chebyt.outer(degs, pnts)
#array([[ 1.    ,  1.    ,  1.    ],
#      [ 0.    ,  0.1   ,  0.2   ],
#      [-1.    , -0.98  , -0.92  ],
#      [-0.    , -0.296 , -0.568 ],
#      [ 1.    ,  0.9208,  0.6928]])

UPDATE

For more parameters, you must broadcast by hand. meshgrid often help for that,spanning each parameter in a dimension. For exemple :

n=3
alpha = numpy.array(range(5))
beta =  numpy.array(range(3))
x = numpy.array(range(2))

data = numpy.meshgrid(n,alpha,beta,x)
values = scipy.special.eval_jacobi(*data)
Community
  • 1
  • 1
B. M.
  • 18,243
  • 2
  • 35
  • 54
  • I searched every possible search query except for "outer" keyword, although it definitely reminded me of the outer product operation. I should read more docs. Thank you! – gsarret Oct 07 '18 at 16:39
  • Whoops, sorry, the problem seems harder than I thought. Please see the updated version of the question. – gsarret Oct 07 '18 at 16:53
1

Reshape the input arguments for broadcasting. In this case, change the shape of degs to be (5, 1) instead of just (5,). The shape (5, 1) broadcast with the shape (3,) results in the shape (5, 3):

In [185]: import numpy as np

In [186]: import scipy.special

In [187]: degs = np.arange(5).reshape(-1, 1)  # degs has shape (5, 1)

In [188]: pnts = np.array([0.0, 0.1, 0.2])

In [189]: values = scipy.special.eval_chebyt(degs, pnts)

In [190]: values
Out[190]: 
array([[ 1.    ,  1.    ,  1.    ],
       [ 0.    ,  0.1   ,  0.2   ],
       [-1.    , -0.98  , -0.92  ],
       [-0.    , -0.296 , -0.568 ],
       [ 1.    ,  0.9208,  0.6928]])
Warren Weckesser
  • 110,654
  • 19
  • 194
  • 214