CODEWITHSUNDEEP
pythonnumpy
galah92
galah92

Reputation: 3991

Inconsistent cartesian product using NumPy

I've got the following definitions of 3 axes vectors:

arena = np.array([[-0.1, 0.1], [-0.1, 0.1], [0.06, 0.14]], dtype=np.float32)
resolution = np.array([0.01, 0.01, 0.01], dtype=np.float32)

axes_size = np.array((arena[:, 1] - arena[:, 0]) / resolution).astype(np.int32)
axes = [np.linspace(*arena[i], axes_size[i], endpoint=False, dtype=np.float32)
        for i in np.arange(axes_size.shape[0])]

Note that axes_size = array([20, 20, 8], dtype=int32).

I'm trying to vectorize & generalize the following piece of code that create the cartesian product of axes:

points = np.zeros([*axes_size, 3])
points[:, :, :, 0] = axes[0][:, None, None]
points[:, :, :, 1] = axes[1][None, :, None]
points[:, :, :, 2] = axes[2][None, None, :]

Note that points.shape = (20, 20, 8, 3).

I tried the following:

points = np.stack(np.meshgrid(*axes), axis=-1)
points = np.array(np.meshgrid(*axes)).T.reshape(-1, 3).T.reshape(*axes_size, 3)

Both lines seems to fail when comparing the results using np.testing.assert_allclose(), and I can't seem to understand why.

Upvotes: 0

Views: 69

Answers (1)

Han-Kwang Nienhuys
Han-Kwang Nienhuys

Reputation: 3254

The problem is that np.meshgrid has a default setting that is not what most people would expect.

Add indexing='ij' to create sane behavior:

points = np.stack(np.meshgrid(*axes, indexing='ij'), axis=-1)

This will pass the np.allclose test.

The default of np.meshgrid is indexing='xy' which reverses the order of the first two axes. This is supposedly convenient for image data, where a 800x400 image (width 800, height 400) should be a shape (400, 800) array.

>>> np.meshgrid(np.arange(20), np.arange(21), np.arange(3))[0].shape
(21, 20, 3)

>>> np.meshgrid(np.arange(20), np.arange(21), np.arange(3), indexing='ij')[0].shape
(20, 21, 3)

If I deal with multidimensional arrays, I try to have different sizes along different axes, e.g. (20, 21, 8, 3) rather than (20, 20, 8, 3) so that I will catch this type of problem.

Upvotes: 1

Related Questions