Say that I have array x and y:
x = numpy.array([1,2,3,4,5,6,7,8,9,10]) # actual content is the a result of another calculation step
There's a formula for y, and each element is based on the previous element, let i denote the index of y, each element is:
y[i] = y[i-1] * 2 + x[i]
When calculating the first element, let y[i-1] = 50. In other words, y should be:
[101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236]
How do I calculate y with numpy?
Lets build a few of the items in your sequence:
y[0] = 2*y[-1] + x[0]
y[1] = 2*y[0] + x[1] = 4*y[-1] + 2*x[0] + x[1]
y[2] = 2*y[1] + x[2] = 8*y[-1] + 4*x[0] + 2*x[1] + x[2]
...
y[n] = 2**(n+1)*y[-1] + 2**n*x[0] + 2**(n-1)*x[1] + ... + x[n]
It may not be immediately obvious, but you can build the above sequence with numpy doing something like:
n = len(x)
y_1 = 50
pot = 2**np.arange(n-1, -1, -1)
y = np.cumsum(pot * x) / pot + y_1 * 2**np.arange(1, n+1)
>>> y
array([ 101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236])
The down side to this type of solutions is that they are not very general: a small change in your problem may render the whole approach useless. But whenever you can solve a problem with a little algebra, it is almost certainly going to beat any algorithmic approach by a far margin.
If you need a recursive computation, if your y[i] should depend on the computed y[i-1] from the same run, then there seems to be no built-in solution in numpy, and you will need to compute it using a simple for loop:
y = np.empty(x.size)
last = 50
for i in range(x.size):
y[i] = last = last * 2 + x[i]
See this question: Is a "for" loop necessary if elements of the a numpy vector are dependant upon the previous element?
Otherwise, you can implement your formula in one line using numpy:
y = np.concatenate(([50], y[:-1])) * 2 + x
Explanation:
y[:-1]
Creates a N-1-sized array: y_0, y_1, ... y_N-1.
np.concatenate(([50], y[:-1]))
Creates a N-sized array with the first element your starting value 50. So this expression basically is your y[i-1].
Then you can do the math element-wise using numpy array arithmetics.
Perhaps the fastest and most concise way is to use scipy.signal.lfilter, which implements exactly the kind of recursive relationship you described:
from scipy.signal import lfilter
import numpy as np
x = np.array([1,2,3,4,5,6,7,8,9,10])
b = [1., 0.]
a = [1., -2.]
zi = np.array([2*50]) # initial condition
y, _ = lfilter(b, a, x, zi=zi)
Result will be np.float64, but you can cast to e.g. np.int32 if that's what you need:
>>> y.astype(np.int32)
array([ 101, 204, 411, 826, 1657, 3320, 6647, 13302, 26613, 53236])
Here is how you do it with Python:
x = [1, 2, 3, 4, 5, 6, 7, 8 ,9, 10]
y = [50]
for i in range(len(x)):
y.append(y[-1] * 2 + x[i])
y = y[1:]
You might want to calculate it from the last element to prevent using the new values for the next i.
This is how to do it with numpy:
import numpy as np
x = np.array([ 1, 2, 3, 4, 5, 6, 7, 8 ,9, 10 ])
y = np.array([ 50 ])
for i in np.arange(len(x)):
y = np.append(
y,
( y[-1] * 2 + x[i] )
)
y = y[1:]
print(y)
