0

I would like to fill a numpy array with values using a function. I want the array to start with one initial value and be filled to a given length, using each previous value in the array as the input to the function.

Each array value i should be (i-1)*x**(y/z).

After a bit of work, I have got to:

import numpy as np
f = np.zeros([31,1])
f[0] = 20
fun = lambda i, j: i*2**(1/3)
f[1:] = np.fromfunction(np.vectorize(fun), (len(f)-1,1), dtype = int)

This fills an array with

[firstvalue=20, 0, i-1 + 1*2**(1/3),...]

I have arrived here having read

https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.fromfunction.html

Use of numpy fromfunction

Most efficient way to map function over numpy array

Fastest way to populate a matrix with a function on pairs of elements in two numpy vectors?

How do I create a numpy array using a function?

But I'm just not getting how to translate it to my function.

Mike
  • 21
  • 8
  • Show us how you would do this calculation with iteration. As my answer shows (and your links) `fromfunction` does not help you do iterative code. It just applies the function to the whole set of indices (with 1 call). `vectorize` also isn't helpful. It's just a convenient way of broadcasting multiple arrays to a function that only takes scalars. Neither of these promise speed. – hpaulj Mar 22 '18 at 21:04
  • Thanks, I realised when I woke up today that it's impossible to vectorise this type of operation, as it would violate causality; the function cannot know what the i-1th value is until it has already calculated it, so it must be a series (iterative) operation to calculate the next value of i. I feel very stupid. I will have a go at iterating it and repost. – Mike Mar 23 '18 at 15:28
  • `for i in np.linspace(1,len(f)-1,len(f)-1, dtype=int): f[i] = f[i-1]*2**(1/3)` – Mike Mar 23 '18 at 16:12
  • Doesn't seem very pythonic, more matlabby, but that's where I'm coming from... – Mike Mar 23 '18 at 16:14
  • Iteration is very pythonic; that's how lists are processed all the time. MATLAB used to require 'vectorization' in the same sense as `numpy`; but they've since added `jit` compilation, that reduces the need to think in terms of whole arrays. ``numba` and `cython` can be used to the same effect. – hpaulj Mar 23 '18 at 19:44
  • See my edits for a `cumprod` solution. – hpaulj Mar 23 '18 at 20:04
  • Very helpful and thorough answer – Mike Mar 25 '18 at 15:51

1 Answers1

1

Except for the initial 20, this produces the same values

np.arange(31)*2**(1/3)

Your iterative version (slightly modified)

def foo0(n):
    f = np.zeros(n)
    f[0] = 20
    for i in range(1,n): 
        f[i] = f[i-1]*2**(1/3)
    return f

An alternative:

def foo1(n):
    g = [20]
    for i in range(n-1):
        g.append(g[-1]*2**(1/3))
    return np.array(g)

They produce the same thing:

In [25]: np.allclose(foo0(31), foo1(31))
Out[25]: True

Mine is a bit faster:

In [26]: timeit foo0(100)
35 µs ± 75 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [27]: timeit foo1(100)
23.6 µs ± 83.6 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

But we don't need to evaluate 2**(1/3) every time

def foo2(n):
    g = [20]
    const = 2**(1/3)
    for i in range(n-1):
        g.append(g[-1]*const)
    return np.array(g)

minor time savings. But that's just multiplying each entry by the same const. So we can use cumprod for a bigger time savings:

def foo3(n):
    g = np.ones(n)*(2**(1/3))
    g[0]=20
    return np.cumprod(g)

In [37]: timeit foo3(31)
14.9 µs ± 14.8 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [40]: np.allclose(foo0(31), foo3(31))
Out[40]: True
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • Thankyou, but even if the first value was 20, this doesn't use the i-1th value to calculate i. – Mike Mar 22 '18 at 17:22