Update conventions

[Korijn]

We need to add a couple of things to the conventions document as became apparent in #22

• Axis conventions (e.g. in look_at we assume positive z-axis is the view axis, positive y-axis is the up axis, etc)
• Which functions support broadcasting/vectorization, and how (e.g. most functions in vector.py support broadcasting a single transform such as a matrix or quaternion over many vectors)

[FirefoxMetzger]

Are we also worrying about the axis/axes to compute along, i.e., providing axis and axes kwargs?

[Korijn]

Are we also worrying about the axis/axes to compute along, i.e., providing axis and axes kwargs?

I thought about this earlier, but no, it is better to just document our assumptions for the input, I think. Look at the existing functions in vector.py to see the current expectations.

In essence (and we could add this to the conventions document) we expect the last axis/axes of input arrays to have a fixed shape, and there may be an arbitrary amount of axes before them. So a vertices array has shape [..., 3] and a quaternion array has shape [..., 4] and so forth

[FirefoxMetzger]

Sure. For matricies this would then be [..., 4, 4] (if we are talking about affine ones), right? This assumes we are working with row-major arrays which works fine with numpy. What about fortran-contigous / column-major ([3, ...])? Do we convert those to row-major internally?

[Korijn]

Sure. For matricies this would then be [..., 4, 4] (if we are talking about affine ones), right?

That's exactly right.

This assumes we are working with row-major arrays which works fine with numpy.

Yes, that's why. I guess we should document this as well (if I didn't already).

What about fortran-contigous / column-major ([3, ...])? Do we convert those to row-major internally?

Users are expected to convert their data prior to calling our functions.

[FirefoxMetzger]

Users are expected to convert their data prior to calling our functions.

We could also use np.asarray(order="C") to convert arrays to fortran order if they aren't already. Since we are converting to ndarray in the beginning of a function anyway this would only be a small change and could avoid confusing bugs.

[FirefoxMetzger]

Also, while on the topic of updating conventions, what is our definition for np.Inf * 0? Numpy uses np.inf * 0 = np.nan, which is important for some of our edge cases where we, for example, perform a translation using an affine matrix.

For example, the following results in an assertion error:

vecA = np.array((np.inf, np.inf, np.inf))
vecB = np.array((0,0,0))

# affine translation matrix by vecB
matrix = np.eye(4)
matrix[-1, :3] = vecB

# translate vecA using affine transform
vecA_hom = np.array((*vecA, 1))
result = (vecA_hom @ matrix)[:3]

# translate vecA using the simplified approach
# expected = (np.inf, np.inf, np.inf)
expected = vecA + vecB

np.testing.assert_allclose(result, expected, equal_nan=True)
# fails, because result == (np.nan, np.nan, np.nan)

(I've been writing some unit-tests for the objects API we've been discussing while waiting for #22 to finalize.)

[Korijn]

My initial thought is to stick to numpy defaults with regards to np.inf and np.nan as much as possible even if it is less convenient.

We could also use np.asarray(order="C") to convert arrays to fortran order if they aren't already. Since we are converting to ndarray in the beginning of a function anyway this would only be a small change and could avoid confusing bugs.

Well yes, but that only works if the ndarray metadata is accurate (right?). It is very easy to create fortran style arrays that are marked as C contiguous. I'm not keen on relying on that metadata...

[FirefoxMetzger]

It is very easy to create fortran style arrays that are marked as C contiguous.

Do you have an example of this? Outside of bugs in numpy, I've never seen numpy's metadata be wrong and - unless I am mistaken - a lot of internals rely on this being correct, i.e., if this metadata doesn't match the actual data layout a lot of numpy functions will stop working as expected.

[almarklein]

Note that on the GPU, constructing a matrix out of elements or vectors follows column major:

mat2x2<f32>(e1,e2,e3,e4)
mat2x2<f32>(c1, c2)

This means that if we use row-major, we must transpose the matrix before flattening it and putting it in the uniform.

This stuff always confuses me, but don't we currently use column major? We don't transpose matrices before uploading, plus our transformations look like `matrix * vec*, which suggests column major, right?

[Korijn]

It is very easy to create fortran style arrays that are marked as C contiguous.

Do you have an example of this? Outside of bugs in numpy, I've never seen numpy's metadata be wrong and - unless I am mistaken - a lot of internals rely on this being correct, i.e., if this metadata doesn't match the actual data layout a lot of numpy functions will stop working as expected.

>>> row_major = np.array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
>>> row_major.flags
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False

>>> col_major = np.array([[0, 3, 6], [1, 4, 7], [2, 5, 8]])
>>> col_major.flags
  C_CONTIGUOUS : True
  F_CONTIGUOUS : False

>>> np.asfortranarray(row_major)
array([[0, 1, 2],
       [3, 4, 5],
       [6, 7, 8]])
>>> np.asfortranarray(row_major).flags
  C_CONTIGUOUS : False
  F_CONTIGUOUS : True

What I'm trying to show here is that the flag maintained by numpy doesn't actually influence the indexing pattern/axes, it only determines the layout of the underlying buffer in memory.

Additionally, the two arrays created above are laid out in row vs column major format in terms of indexing patterns/axes, but the flag is c contiguous for both (of course). It's only set when you explicitly set it, or when you convert to fortran explicitly.

It could be that I'm just ignorant of some mechanism in numpy that can prove useful to us here. Please point it out to me, if so . :)

To sum up; I think we should just document that we expect input arrays to be c contiguous, and to have row major indexing pattern/axes.

[Korijn]

Note that on the GPU, constructing a matrix out of elements or vectors follows column major:

mat2x2<f32>(e1,e2,e3,e4)
mat2x2<f32>(c1, c2)

This means that if we use row-major, we must transpose the matrix before flattening it and putting it in the uniform.

This stuff always confuses me, but don't we currently use column major? We don't transpose matrices before uploading, plus our transformations look like `matrix * vec*, which suggests column major, right?

Our conventions are set in stone in vector_apply_matrix where vectors and matrices come together:

https://github.com/pygfx/pylinalg/blob/main/pylinalg/func/vector.py#L104-L108

In that function, we can slice off the last row of the matrix, because it holds no useful information. The last column holds the translation component.

Indeed, when you perform matrix * vec, you need the vector to be a column, as in this reflection transform example from wikipedia:

However, our vector array is row major (every row is a vector), so we transpose the matrix after slicing off the last row, to be able to broadcast the transformation.

So actually, our matrices are laid out column major, and the vectors are row major! 👀

(And yet they are all C_CONTIGUOUS)

[FirefoxMetzger]

This means that if we use row-major, we must transpose the matrix before flattening it and putting it in the uniform.

That depends on how wgpu is implemented. row vs column major goes back to C vs fortran in that C uses row-major (also called c-contiguous) and fortran uses colum-major (also called f-contiguous). Everything newer than C tends to use row-major, and I would have expected that wgpu handles this conversion internally given how uncommon column-major is on the CPU.

From what I have seen so far, we are consistently using row-major, so if we are not transposing at the moment we hand over to the GPU, we are probably fine :)

This stuff always confuses me, but don't we currently use column major? We don't transpose matrices before uploading, plus our transformations look like `matrix * vec*, which suggests column major, right?

Which way we do the multiplication doesn't matter in regards to row/column-major. The reason I suggest to left-multiply vectors is that we use the convention that batches are stacked along the first dimension, i.e., a batch of n-vertices has shape (N, 3). Because of this, users either have to transpose the transformation matrix or the batch of vectors when multiplying: vectors @ matrix.T or matrix @ vectors.T. vectors @ matrix.T is a bit nicer though, because the WorldObject controls the matrix. This means we can directly construct matrix.T and call it transformation_matrix and get vectors @ transformation_matrix as the user-facing API, which is neat.

Row/column-major relates to the memory layout of the array, i.e., the order in which elements appear in the (flat) buffer that underlies the array. Here is a quick 2x2 example:

>>> foo = np.array([[1, 2], [3, 4]], order="C")
>>> foo
array([[1, 2],
       [3, 4]])
>>> foo.ravel(order="K")  # order=K keeps current memory layout
array([1, 2, 3, 4])
>>> bar = np.array([[1, 2], [3, 4]], order="F")
>>> bar
array([[1, 2],
       [3, 4]])
>>> bar.ravel(order="K")
array([1, 3, 2, 4])

You can work with either order just the same in numpy, completely ignore row/column-major and be absolutely fine (numpy handles the rest). The only time this matters is when you need to do (1) performance optimization [some ops are faster in column-major, most are not], (2) serialization [when you deserialize you need to know what format the incoming data is in], or (3) you work on compatibility between low-level libraries and wish to pass arrays by reference.

[FirefoxMetzger]

So actually, our matrices are laid out column major, and the vectors are row major! 👀 (And yet they are all C_CONTIGUOUS)

Oh! Now I get it, we are talking about different things @Korijn. You are talking about multi-indexing (always row-major) whereas I am talking about memory layout (C/F-contiguous). The latter is relevant for performance in pylinalg because we do a lot of matmul with small matrices.

For large and arbitrary matrices (Nx111x200 @ Nx200x77), doing a @ b column-major is faster:

For smaller matrices (Nx3x4 @ Nx4x1), i.e. a batch of transformation matrices times a batch of vectors/vertices, row-major is faster on modern architecture:

import numpy as np
import matplotlib.pyplot as plt
from timeit import timeit

print(np.show_config())

rng = np.random.default_rng(42)
matrix_size = np.linspace(10, int(1e5), 200).astype(int)


n_repeats = 100

c_times = []
f_times = []

for idx in range(len(matrix_size)):
    c_times.append(
        timeit(
            "(matrix @ vectors)",
            setup=(
                "matrix=rng.random((idx, 3, 4));"
                "vectors=rng.random((idx, 4, 1));"
            ),
            globals=globals(),
            number=n_repeats,
        )
        / n_repeats
    )
    f_times.append(
        timeit(
            "(matrix @ vectors)",
            setup=(
                "matrix = np.asarray(rng.random((idx, 3, 4)), order='F');"
                "vectors = np.asarray(rng.random((idx, 4, 1)), order='F')"
            ),
            globals=globals(),
            number=n_repeats,
        )
        / n_repeats
    )

plt.plot(matrix_size, c_times, label="Column-Major")
plt.plot(matrix_size, f_times, label="Row-Major")
plt.gca().set_title("Performance (log-log)")
plt.gca().set_xlabel("Number of Elements")
plt.gca().set_ylabel("Time per run")
plt.xscale("log")
plt.yscale("log")
plt.legend(title="Method")
plt.tight_layout()
plt.savefig("distance_timings.png")

For a single 4x3 matrix broadcasted along a batch of vectors (Nx4 @ (3x4).T) performance is the same::

import numpy as np
import matplotlib.pyplot as plt
from timeit import timeit

print(np.show_config())

rng = np.random.default_rng(42)
matrix_size = np.logspace(1, 5, 200).astype(int)


n_repeats = 100

c_times = []
f_times = []

for idx in range(len(matrix_size)):
    c_times.append(
        timeit(
            "(vectors @ matrix.T)",
            setup=(
                "matrix=rng.random((3, 4));"
                "vectors=rng.random((idx, 4));"
            ),
            globals=globals(),
            number=n_repeats,
        )
        / n_repeats
    )
    f_times.append(
        timeit(
            "(vectors @ matrix.T)",
            setup=(
                "matrix = np.asarray(rng.random((3, 4)), order='F');"
                "vectors = np.asarray(rng.random((idx, 4)), order='F')"
            ),
            globals=globals(),
            number=n_repeats,
        )
        / n_repeats
    )

plt.plot(matrix_size, c_times, label="Column-Major")
plt.plot(matrix_size, f_times, label="Row-Major")
plt.gca().set_title("Performance (log-log)")
plt.gca().set_xlabel("Number of Elements")
plt.gca().set_ylabel("Time per run")
plt.xscale("log")
plt.yscale("log")
plt.legend(title="Method")
plt.tight_layout()
plt.savefig("distance_timings.png")

Bottomline: Some people may prefer to operate on column-major arrays, and some on row-major arrays depending on the application scenario. Since we are trying to build something that can be re-used I am wondering to what extent we want to support these user groups vs optimize for pygfx.

[almarklein]

Attempt to clarify the row-major thing.

Row-major can mean two things

• The memory layout: whether rows are contiguous in memory, or the columns.
• Whether vectors are considered columns or rows. E.g. whether m[i] gets a row or a column from the matrix.

TL;DR;

Numpy is row-major for both interpretations, wgsl is column-major. And this is fine.

We could consider transposing the data where the uniform is set, so that in the shaders we can multiply in the same way as we do on the CPU.

Numpy is row-major for both (by default)

>>> m = np.array([[1, 2],[3, 4]])
>>> m
array([[1, 2],
       [3, 4]])
>>> m[0]
array([1, 2])

WGPU / WGSL is column-major for both

The docs clearly state that it's column-major for the latter interpretation of row/col-major:

So what about the memory layout? I did a quick test, putting color values in a matrix like this:

This produces a matrix:

r  g  b  0
0  0  0  0 
0  0  0  0
0  0  0  0

and then to get the correct result, the shader has to sample the colors like this:

From the docs we know that m[0] is column from the matrix, so the matrix layout on the GPU is this:

r  0  0  0
g  0  0  0 
b  0  0  0
0  0  0  0

This expected, because because the raw data is passed unchanged, first interpreted as row-major, now interpreted as column-major. We've effectively transposed the matrix as we uploaded it.

From the code in Vector4.apply_matrix4() it indeed looks like we multiply matrices like v * M. And indeed, in our shaders we consistently do M * v.

[almarklein]

Bottom line of the above: just do whatever works best on the CPU - what happens on the GPU is irrelevant for the current discussion :)