numpy.fromfunction with specified non-arange-like array

Question

I have defined a function;

def f(x,y):
    return (1+x)**3.2 * np.cos(y**2.31+x)

and have the following two arrays:

x = np.linspace(10,20,10)
y = np.linspace(5,8,5)

Now I wish to create a 10 x 5 matrix where each (i,j) element is given by f(x[i],x[j]). This is easily achieved with a for loop:

result = np.zeros(len(x)*len(y)).reshape((len(x),len(y)))
for i in range(len(x)):
    for j in range(len(y)):
        result[i][j] = f(x[i],y[j])

However, I want to significantly increase computation time so I am looking for a non-for loop approach. It would seem that np.fromfunction brings me partly there:

result = np.fromfunction(lambda x, y: f(x,y), (10, 5), dtype=int)

However, this takes x and y to be arange-like, i.e. this takes elements x = 0,1,...,9 and y = 0,1,...,4. Is there a way to tell this function that I want it to take my previously defined arrays for x and y. If this is not possible with np.fromfunction, is there another way to achieve the same result I get with my for loop, but with less computation time?

Certainly. See the edited question. It's a complicated function with cosines and powers. This is just an example as the real function I use is even more complicated. — The One
– The One, Commented Apr 3, 2020 at 17:11
fromfunction is not a substitute for loops! It only works if your function already works with whole arrays. I know the docs can be confusing, especially if you come with preconceived notions, but when in doubt, look at the code as well. — hpaulj
– hpaulj, Commented Apr 3, 2020 at 19:21

hpaulj · Accepted Answer · 2020-04-03 22:10:42Z

If your function is written so it accepts whole arrays, you can just pass the x and y so they broadcast together.

In [472]: x = np.linspace(10,20,10) 
     ...: y = np.linspace(5,8,5)                                                               
In [473]: def f(x,y): 
     ...:     return (1-x)**3.2 * np.cos(y**2.31+x) 
     ...:

Make x (10,1) array which then broadcasts with the (5,) y to give a (10,5) result:

In [474]: f(x[:,None], y)                                                                      
/usr/local/bin/ipython3:2: RuntimeWarning: invalid value encountered in power

Out[474]: 
array([[nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan],
       [nan, nan, nan, nan, nan]])

All fromfunction does is generate indices and pass them to your function:

In [478]: np.indices((3,4))                                                                    
Out[478]: 
array([[[0, 0, 0, 0],
        [1, 1, 1, 1],
        [2, 2, 2, 2]],

       [[0, 1, 2, 3],
        [0, 1, 2, 3],
        [0, 1, 2, 3]]])
#function(*indices...)

If your function expects values, not the i,j indices fromfunction does nothing for you!. it just forces you to use the extra layer of indexing f(x[i],y[j]). If your function handle arrays, fromfunction does not speed up the calculations. And if your function can only accept scalar, fromfunction does not work at all.

Why the nan?

In [489]: f(x[0],y[0])                                                                         
/usr/local/bin/ipython3:2: RuntimeWarning: invalid value encountered in double_scalars

Out[489]: nan

===

With the corrected f:

In [494]: def f(x,y): 
     ...:     return (1+x)**3.2 * np.cos(y**2.31+x) 
     ...:                                                                                      
In [495]: x,y                                                                                  
Out[495]: 
(array([10.        , 11.11111111, 12.22222222, 13.33333333, 14.44444444,
        15.55555556, 16.66666667, 17.77777778, 18.88888889, 20.        ]),
 array([5.  , 5.75, 6.5 , 7.25, 8.  ]))

ix_ can be used to add the broadcasting dimension to x:

In [496]: np.ix_(x,y)                                                                          
Out[496]: 
(array([[10.        ],
        [11.11111111],
        [12.22222222],
        [13.33333333],
        [14.44444444],
        [15.55555556],
        [16.66666667],
        [17.77777778],
        [18.88888889],
        [20.        ]]), 
 array([[5.  , 5.75, 6.5 , 7.25, 8.  ]]))

And those arrays can then be passed to f (as I did in [474]):

In [497]: f(*np.ix_(x,y))                                                                      
Out[497]: 
array([[  1322.50123614,  -1353.39115671,  -1702.96842965,
          2039.59858593,   2149.99822839],
       [ -1268.79441897,   1220.20345147,    572.47038806,
           401.57578676,   1322.04810644],
       [ -3873.89288421,   3872.45293387,   3741.19176848,
         -3203.15416244,  -2320.43820171],
       [ -2274.83634272,   2356.4885057 ,   3316.20257628,
         -4368.07920439,  -4932.15975126],
       [  3805.32723978,  -3710.95441159,  -2413.75855025,
           343.98082438,  -1742.79381001],
       [  7825.16462628,  -7850.07869311,  -7934.59775682,
          7309.08041012,   5891.07144208],
       [  2696.9980712 ,  -2869.31389573,  -4956.11491324,
          7455.17663401,   9114.69519122],
       [ -8800.47513177,   8651.89940176,   6527.58767232,
         -2896.13626079,   1015.66697767],
       [-13326.47899275,  13419.77502789,  14202.98227539,
        -13981.08503081, -12233.58338808],
       [ -1484.09855587,   1795.12675239,   5660.63236478,
        -10620.54747276, -14370.53381538]])

Applied to 3 arrays:

In [503]: np.ix_([100,200],[10,20,30],[1,2,3,4])                                               
Out[503]: 
(array([[[100]],

        [[200]]]), array([[[10],
         [20],
         [30]]]), array([[[1, 2, 3, 4]]]))

produces a (2,1,1), (1,3,1) and (1,1,4) arrays, which would broadcast together to produce a (2,3,4) result.

def ff(x,y,z): return x+y+z  
In [504]: ff(*_)                                                                               
Out[504]: 
array([[[111, 112, 113, 114],
        [121, 122, 123, 124],
        [131, 132, 133, 134]],

       [[211, 212, 213, 214],
        [221, 222, 223, 224],
        [231, 232, 233, 234]]])

The nan is due to 1-x being negative in my example, making the (1-x)**3.2 complex. Changing this to (1+x)**3.2 fixes this. I shall update my post account for to this.
Also, this is a simplified version of my problem. In truth, my function looks more complicated than this (e.g. it takes actually 4 input values instead of 2), leading to much more complications that make np.fromfunction perhaps less appealing. And it isn't as evident as to how I could use the broadcasting you describe for more than 2 input values. I think I'd best just explain these complications in a new post.
@TheOne. You should still figure out how to properly vectorize your code. Post the real function if you need to. Calling fromfunction will run a glorified for loop for you under the hood.

kederrac · Accepted Answer · 2020-04-03 17:19:54Z

2

you can use:

np.fromfunction(lambda i, j:  f(x[i],y[j]), (len(x),len(y)), dtype=int)

you have to specify dtype=int otherwise the data-type of the coordinate will be float which can not be used as indices in f(x[i],y[j])

edited Apr 3, 2020 at 17:19

answered Apr 3, 2020 at 17:12

kederrac

17.4k6 gold badges36 silver badges58 bronze badges

Collectives™ on Stack Overflow

numpy.fromfunction with specified non-arange-like array

2 Answers 2

4 Comments

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

4 Comments

Comments

Your Answer

Sign up or log in

Post as a guest

Related