1

I'm fairly new to NumPy, so it's quite possible that I'm missing something fundamental. Don't hesitate to ask "stupid" questions about "basic" things!

I'm trying to write some functions that manipulate vectors. I'd like them to work on single vectors, as well as on arrays of vectors, like most of NumPy's ufuncs:

import math
import numpy

def func(scalar, x, vector):
    # arbitrary function
    # I'm NOT looking to replace this with numpy.magic_sum_multiply()
    # I'm trying to understand broadcasting & dtypes
    return scalar * x + vector

print(func(
    scalar=numpy.array(2),
    x=numpy.array([1, 0, 0]),
    vector=numpy.array([1, 0, 0]),
))
# => [3 0 0], as expected

print(func(
    scalar=numpy.array(2),
    x=numpy.array([1, 0, 0]),
    vector=numpy.array([[1, 0, 0], [0, 1, 0]]),
))
# => [[3 0 0], [2 1 0]], as expected. x & scalar are broadcasted out to match the multiple vectors

However, when trying to use multiple scalars, things go wrong:

print(func(
    scalar=numpy.array([1, 2]),
    x=numpy.array([1, 0, 0]),
    vector=numpy.array([1, 0, 0]),
))
# => ValueError: operands could not be broadcast together with shapes (2,) (3,)
# expected: [[2 0 0], [3 0 0]]

I'm not entirely surprised be this. After all, NumPy has no idea that I'm working with vectors that are an single entity, and not an arbitrary dimension.

I can solve this ad-hoc with some expand_dims() and/or squeeze() to add/remove axes, but that feels hacky...

So I figured that, since I'm working with vectors that are a single "entity", dtypes may be what I'm looking for:

vector_dtype = numpy.dtype([
    ('x', numpy.float64),
    ('y', numpy.float64),
    ('z', numpy.float64),
])
_ = numpy.array([(1, 0, 0), (0, 1, 0)], dtype=vector_dtype)
print(_.shape)  # => (2,), good, we indeed have 2 vectors!

_ = numpy.array((1, 0, 0, 7), dtype=vector_dtype)
# Good, basic checking that I'm staying in 3D
# => ValueError: could not assign tuple of length 4 to structure with 3 fields.

However, I seem to loose basic math capabilities:

print(2 * _)
# => TypeError: The DTypes <class 'numpy.dtype[void]'> and <class 'numpy.dtype[uint8]'> do not have a common DType. For example they cannot be stored in a single array unless the dtype is `object`.

So my main question is: How do I solve this?

  • Is there some numpy.magic_broadcast_that_understands_what_I_mean() function?
  • Can I define math-operators (such as addition, ...) on the vector-dtype?
Niobos
  • 880
  • 4
  • 15

2 Answers2

0

ufuncs obey the same broadcasting rules as the operators. And your own function, written with numpy operators and ufuncs have to work with those as well. Your function could tweak the dimensions to translate inputs to something works with the rest of numpy. (Writing your own ufuncs is an advanced topic.)

In [64]: scalar=numpy.array([1, 2])
    ...: x=numpy.array([1, 0, 0])
    ...: vector=numpy.array([1, 0, 0])
In [65]: scalar * x + vector
Traceback (most recent call last):
  File "<ipython-input-65-ad4a73833616>", line 1, in <module>
    scalar * x + vector
ValueError: operands could not be broadcast together with shapes (2,) (3,) 

The problem is the multiplication; regardless of what you call it, scalar is a (2,) shape array, which does not work with a (3,) array.

In [68]: scalar*x
Traceback (most recent call last):
  File "<ipython-input-68-0d21729ffa15>", line 1, in <module>
    scalar*x
ValueError: operands could not be broadcast together with shapes (2,) (3,) 

But what do you expect to happen? What shape should the result have?

If scalar is a (2,1) shaped array, then by broadcasting this result is (2,3) - taking the 2 from scalar and 3 from the other arrays:

In [76]: scalar[:,None] * x + vector
Out[76]: 
array([[2, 0, 0],
       [3, 0, 0]])

This is standard numpy broadcasting, and there's nothing "hacky" about it.

I don't know what you mean by calling scalar a 'single entity'.

Structured array is a convenient way of putting arrays with diverse dtypes into one structure. Or to access "columns" of convenient 'names'. But you can't perform math across the fields of such an array.

In [70]: z=np.array([(1, 0, 0), (0, 1, 0)], dtype=vector_dtype)
In [71]: z
Out[71]: 
array([(1., 0., 0.), (0., 1., 0.)],
      dtype=[('x', '<f8'), ('y', '<f8'), ('z', '<f8')])
In [72]: z.shape
Out[72]: (2,)
In [73]: z.dtype
Out[73]: dtype([('x', '<f8'), ('y', '<f8'), ('z', '<f8')])
In [74]: z['x']
Out[74]: array([1., 0.])
In [75]: 2*z['x']                   # math on a single field
Out[75]: array([2., 0.])

note

There is a np.vectorize function. It takes a function that accepts only (true) scalar arguments, and applies array arguments according to the standard broadcasting rules. So even if your func was implemented with it, you'd still have to use arguments as I did in [70]. Sometimes it's convenient, but it's better to use standard numpy functions and operators where possible - better and much faster.

hpaulj
  • 221,503
  • 14
  • 230
  • 353
0

How do I solve this?

You are after the first-argument vectorized version of func, let's call it vfunc(vfunc is not "vectorization" stricto sensu, since the vectorization job in done internally.)

#               v
def vfunc(scalars, x, vector):
    #           ^
    return numpy.vstack([  # Assuming that's the shape you want.
        scalar * x + vector for scalar in scalars
    ])

print(vfunc(
    scalars = [2],  # No need for array instance actually
    x       = numpy.array([1, 0, 0]),
    vector  = numpy.array([1, 0, 0]),
))
# => [3 0 0], as expected

print(vfunc(
    scalars = [2],
    x       = numpy.array([1, 0, 0]),
    vector  = numpy.array([[1, 0, 0], [0, 1, 0]]),
))
# => [[3 0 0], [2 1 0]], as expected

print(vfunc(
    scalars = [1, 2],
    x       = numpy.array([1, 0, 0]),
    vector  = numpy.array([1, 0, 0]),
))
# => # expected: [[2 0 0], [3 0 0]]

[...] dtypes may be what I'm looking for

No it is not.

Is there some numpy.magic_broadcast_that_understands_what_I_mean()

Yes. It is called numpy.vectorize but it is not worth it.

As it reads in the documentation:

The vectorize function is provided primarily for convenience, not for performance. The implementation is essentially a for loop.

keepAlive
  • 6,369
  • 5
  • 24
  • 39
  • 1
    Thank you, I learned a lot from that answer. I did know about `vectorize`, but indeed, the whole point was to perform math in parallel, not just loop over it. – Niobos Feb 09 '21 at 07:56