Is this a Python/Numpy bug or a subtle gotcha?

Question

Consider the following two implemtations of the same piece of code. I would have thought they are identical but they are not.

Is this a Python/Numpy bug or a subtle gotcha? If the latter, what rule would make it obvious why it does not work as expected?

I was working with multiple arrays of data and having to process each array item by item, with each array manipulated by a table depending on it's metadata.

In the real world example 'n' is multiple factors and offsets but the following code still demonstrates the issue that I was getting the wrong result in all but one case.

import numpy as np

# Change the following line to True to show different behaviour
NEEDS_BUGS = False  # Changeme

# Create some data
data = np.linspace(0, 1, 10)
print(data)

# Create an array of vector functions each of which does a different operation on a set of data
vfuncd = dict()

# Two implementations
if NEEDS_BUGS:
    # Lets do this in a loop because we like loops - However WARNING this does not work!!
    for n in range(10):
        vfuncd[n] = np.vectorize(lambda x: x * n)
else:
    # Unwrap the loop - NOTE: Spoiler - this works
    vfuncd[0] = np.vectorize(lambda x: x * 0)
    vfuncd[1] = np.vectorize(lambda x: x * 1)
    vfuncd[2] = np.vectorize(lambda x: x * 2)
    vfuncd[3] = np.vectorize(lambda x: x * 3)
    vfuncd[4] = np.vectorize(lambda x: x * 4)
    vfuncd[5] = np.vectorize(lambda x: x * 5)
    vfuncd[6] = np.vectorize(lambda x: x * 6)
    vfuncd[7] = np.vectorize(lambda x: x * 7)
    vfuncd[8] = np.vectorize(lambda x: x * 8)
    vfuncd[9] = np.vectorize(lambda x: x * 9)

# Prove we have multiple different vectorised functions
for k, vfunc in vfuncd.items():
    print(k, vfunc)

# Do the work
res = {k: vfuncd[k](data) for k in vfuncd.keys()}

# Show the result
for k, r in res.items():
    print(k, r)

Answer 1

I don't know what exactly you're trying to achieve and if it's a bad idea or not (in terms of np.vectorize), but the issue you're facing is because of the way python makes closures. Quoting from an answer to the linked question:

Scoping in Python is lexical. A closure will always remember the name and scope of the variable, not the object it's pointing to. Since all the functions in your example are created in the same scope and use the same variable name, they always refer to the same variable.

in other words when you make that closure over n, you're not actually closing off the state of n, just the name. So when n changes, the value in your closure also changes. This is quite unexpected to me, but others find it natural.

Here is one fix using partial:

from functools import partial
.
.
.

def func(x, n):
    return x * n

for n in range(10):
    vfuncd[n] = np.vectorize(partial(func, n=n))

Or another using a factory method

def func_factory(n):
    return lambda x: x * n

for n in range(10):
    vfuncd[n] = np.vectorize(func_factory(n))

Answer 2

It seems that the python variable n is bound to the vectorized expression:

for n in range(10):
    vfuncd[n] = np.vectorize(lambda x: x * n)

This fixes it as it creates a new object with which to bind:

for n in range(10):
    vfuncd[n] = np.vectorize(lambda x: x * np.scalar(n))

In fact this has implications in terms of performance as I assume the value of the python variable would have to be fetched repeatedly.

Answer 3

In [13]: data = np.linspace(0,1,11)                                                     

Since the data array can be multiplied with a simple:

In [14]: data*3                                                                         
Out[14]: array([0. , 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3. ])

we don't need the complication of np.vectorize to see the closure issue. A simple lambda is enough.

In [15]: vfuncd = {} 
    ...: for n in range(3): 
    ...:     vfuncd[n] = lambda x:x*n 
    ...:                                                                                
In [16]: vfuncd                                                                         
Out[16]: 
{0: <function __main__.<lambda>(x)>,
 1: <function __main__.<lambda>(x)>,
 2: <function __main__.<lambda>(x)>}
In [17]: {k:v(data) for k,v in vfuncd.items()}                                          
Out[17]: 
{0: array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ]),
 1: array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ]),
 2: array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ])}

We won't get the closure problem if we use a proper numpy "vectorization":

In [18]: data * np.arange(3)[:,None]                                                    
Out[18]: 
array([[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ],
       [0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ]])

Or a simple iteration is we need a dictionary:

In [20]: {k:data*k for k in range(3)}                                                   
Out[20]: 
{0: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
 1: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 2: array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ])}

np.vectorize has a speed disclaimer. But it is justified where the function only takes scalar inputs, and we want the flexibility of numpy broadcasting - i.e. for 2 or more arguments.

Creating multiple vectorize is clearly an 'anti-pattern'. I'd rather see one vectorize with the appropriate arguments:

In [25]: f = np.vectorize(lambda x,n: x*n)                                              
In [26]: {n: f(data,n) for n in range(3)}                                               
Out[26]: 
{0: array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]),
 1: array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ]),
 2: array([0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ])}

That f can also produce the array Out[18] (but is slower):

In [27]: f(data, np.arange(3)[:,None])                                                  
Out[27]: 
array([[0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. ],
       [0. , 0.2, 0.4, 0.6, 0.8, 1. , 1.2, 1.4, 1.6, 1.8, 2. ]])