1

I am trying to use Python/Numpy vectorized functions to reduce for loops.

My function call looks like this

out_vectors = v_calculation(
                            in_vectors,
                            p
)

The vectorized function definition is

v_calculation = np.vectorize(
                     my_calculation,
                     signature='(j,i),(i)->()'
)

in_vectors is an array of shape

(3,6,200,3)

But the first 2 dimensions (3,6) can be anything. These are the loop dimensions.

p is an array of shape

(3,)

My calculation is this

def my_calculation(in_vector, p):
    """
    total magnetic field from Biot-Savart's law
    """
    out_vector = np.zeros((3,))

    l_vector = in_vector[1:, :] - in_vector[:-1, :]
    r_vector = (in_vector[:-1, :] + l_vector / 2) - p

    out_vector = np.sum(np.cross(l_vector, r_vector) / \
                        np.linalg.norm(r_vector) ** 3,
                        axis=0
    )

    return out_vector

In this function, in_vector is an array of shape (200, 3) and p is the same shape (3,). out_vector shape is (3,). This is correct.

out_vectors, the result of the vectorized function should be (6,3). This should be the results of my_calculation summed over the first dimension of the input_vectors (3 in this case) for each of the second dimension of the input_vectors (6 in this case). The second dimension of the result is 3 (x, y, z components for the vector), same as the dimension of p and the fourth dimension of input_vectors. I hope this is all clear.

My code is failing in the vectorized function call

Stacktrace

~/path/to/my/code.py in calculate_vectors(mgr)
    588         out_vectors = v_calculation(
    589                                 in_vectors,
--> 590                                 p
    591         )

~/miniconda/lib/python3.7/site-packages/numpy/lib/function_base.py in __call__(self, *args, **kwargs)
   1970             vargs.extend([kwargs[_n] for _n in names])
   1971 
-> 1972         return self._vectorize_call(func=func, args=vargs)
   1973 
   1974     def _get_ufunc_and_otypes(self, func, args):

~/miniconda/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call(self, func, args)
   2036         """Vectorized call to `func` over positional `args`."""
   2037         if self.signature is not None:
-> 2038             res = self._vectorize_call_with_signature(func, args)
   2039         elif not args:
   2040             res = func()

~/miniconda/lib/python3.7/site-packages/numpy/lib/function_base.py in _vectorize_call_with_signature(self, func, args)
   2100 
   2101             for output, result in zip(outputs, results):
-> 2102                 output[index] = result
   2103 
   2104         if outputs is None:

ValueError: setting an array element with a sequence.
gordon macmillan
  • 123
  • 3
  • 13
  • A quick comment - `np.vectorize` is slower than a more direct iteration. And with a signature it is even slower. So use this for convenience, not speed. Also I have not seen, or answered, many questions involving the `signature` option. It's a relatively new feature. – hpaulj Mar 05 '19 at 19:34
  • The error means that `output[index]` is a single element slot of the `output` array, but `result` is some sort of array. Which suggests that the signature does not match your function. (But I haven't examined either). – hpaulj Mar 05 '19 at 19:36
  • `np.array((3,))` does not create a shape (3,) array. However that statement is useless. We don't 'define' variables ahead of time in Python. `out_vector` gets the value assigned to it from the `np.sum(...)` expression. – hpaulj Mar 05 '19 at 21:57
  • Oh yeah, that is true. I meant to pre-allocate it using np.zeros. https://stackoverflow.com/questions/33333683/should-i-preallocate-a-numpy-array shows the use of pre-allocation. i was under the impression it is used quite commonly. Is it not? – gordon macmillan Mar 05 '19 at 22:07
  • You preallocate only if using the `out` parameter of a `ufunc` (as illustrated in the link), or if doing indexed assignment, e.g. `output[i] = fun(i)`. `a = np.zeros((3,))` creates a 3 element array. `output=...` assigns a new object to the variable, wiping out any previous assignments. – hpaulj Mar 05 '19 at 22:11

1 Answers1

0

This works for me. Note I changed the return signature, to match the shared final dimension of both inputs.

In [54]: A = np.arange(12).reshape(4,3); b = np.arange(3)                       
In [55]: my_calculation(A,b)                                                    
Out[55]: array([0., 0., 0.])
In [56]: f = np.vectorize(my_calculation, signature='(j,i),(i)->(i)')           
In [57]: f(A,b)                                                                 
Out[57]: array([0., 0., 0.])
In [58]: f([A,A,A],b)                                                           
Out[58]: 
array([[0., 0., 0.],
       [0., 0., 0.],
       [0., 0., 0.]])
hpaulj
  • 221,503
  • 14
  • 230
  • 353
  • Thank you hpaulj. This worked exactly like I wanted it to. I apparently was thinking about the output wrong in that I should have been expecting a size (3,6,3) array. – gordon macmillan Mar 06 '19 at 19:50