4

There has been some related questions for over 7 years, but I raise this issue again as I could see no 'numpy' way iteration method provided.

The task is as follows: If I have an numpy array 'arr' and have a custom function 'fn', how could I iteratively apply 'fn' over the 'arr'? 'fn' cannot be constructed by ufunc tools.

Below is the toy_code I come up with this:

import numpy as np

r_list = np.arange(1,6,dtype=np.float32)
# r_list = [1. 2. 3. 4. 5.]
r_list_extended = np.append([0.],r_list)
R_list_extended = np.zeros_like(r_list_extended)
print(r_list)
gamma = 0.99
pv_mc = lambda a, x: x+ a*gamma

# no cumsum, accumulate available
for i in range(len(r_list_extended)):
    if i ==0: continue
    else: R_list_extended[i] = pv_mc(R_list_extended[i-1],r_list_extended[i])

R_list = R_list_extended[1:]
print(R_list)
# R_list == [ 1.          2.99        5.9601      9.900499   14.80149401]

r_list is an array of r for each time. R_list is a cumulative sum of discounted r. Assume r_list and R_list are reverted beforehand. The loop in above does R[t] : r[t] + gamma * R[t-1]

I do not think this is the best way to utilize numpy.... If one can utilize tensorflow, then tf.scan() does the job as below:

import numpy as np
import tensorflow as tf

r_list = np.arange(1,6,dtype=np.float32)
# r_list = [1. 2. 3. 4. 5.]
gamma = 0.99
pv_mc = lambda a, x: x+ a*gamma
R_list_graph = tf.scan(pv_mc, r_list, initializer=np.array(0,dtype=np.float32))

with tf.Session() as sess:
    R_list = sess.run(R_list_graph, feed_dict={})
    print(R_list)
    # R_list = [ 1.        2.99      5.9601    9.900499 14.801495]

Thanks in advance for your help!

P-Gn
  • 23,115
  • 9
  • 87
  • 104
sdr2002
  • 542
  • 2
  • 6
  • 16

1 Answers1

4

You could use np.frompyfunc, whose documentation is somewhat obscure.

import numpy as np

r_list = np.arange(1,6,dtype=np.float32)
# r_list = [1. 2. 3. 4. 5.]
r_list_extended = np.append([0.],r_list)
R_list_extended = np.zeros_like(r_list_extended)
print(r_list)
gamma = 0.99
pv_mc = lambda a, x: x+ a*gamma
ufunc = np.frompyfunc(pv_mc, 2, 1)
R_list = ufunc.accumulate(r_list, dtype=np.object).astype(float)
print(R_list)
P-Gn
  • 23,115
  • 9
  • 87
  • 104
  • 1
    I've recommended `frompyfunc` a number of times, but didn't realize that it returns a `ufunc` with methods like `accumulate`. Elsewhere I found it gave, at best, a 2x speedup compared to explicit loops. Here it seems to give a 1.5 speedup, and it is still 20x slower than `cumsum` (`gamma=1` case). – hpaulj May 21 '18 at 21:07
  • The `dtype=object` is an important parameter that easy to miss. `frompyfunc` returns a object array. – hpaulj May 21 '18 at 21:19
  • Oh this is great. Yes I wandered around sth like 'vectorize' or sth but np.frompyfunc is the best bet I can do with this. As @hpaulj mentioned, I do not understand/like why np.object is only allowed to accumulate because it makes the calculation slow & need another property changer .astype(float) in the end. But again, dtype only accept np.object (not np.float) so I think this answer may be the way to go. But let me wait a bit for what others think – sdr2002 May 22 '18 at 09:32
  • 1
    "The returned ufunc always returns PyObject arrays." - it's what they said,,, cannot avoid of using 'object' dtype – sdr2002 May 22 '18 at 11:26
  • *update: for this specific calculation, scipy.signal.lfilter make your life easier – sdr2002 Jun 18 '18 at 12:54
  • scipy.signal.lfilter([1.0], [1.0, -gamma], r_list[::-1])[::-1] # refer to patcoady/trpo/train.py github – sdr2002 Jun 18 '18 at 12:54