mpusz / mp-units

The quantities and units library for C++
https://mpusz.github.io/mp-units/
MIT License
1.04k stars 80 forks source link

Implement simple vector and tensor types #594

Open mpusz opened 1 month ago

mpusz commented 1 month ago

To not depend on external libraries and to provide correct, fast, and simple examples of LA usage with quantities we need simple vector and tensor types (up to 3 dimensions) together with the following operations:

PatrickKa commented 1 month ago

Hi. I have a few questions regarding this issue:

  1. Should tensors of arbitrary rank be supported, or will the rank be capped, since most "normal" physics does not require more than rank-4 tensors? (Even those are rather rare. The only two rank-4 tensors I could think of off the top of my head were the elasticity tensor and the Riemann curvature tensor)
  2. When you speak of inner products between tensors or a tensor and a vector, do you actually mean contractions? So, e.g., if a_ij, and b_jk denote the components of two rank-2 tensors is (a ⋅ b)_ik = Σ_j a_ij b_jk?
  3. What do you mean with "scalar product of two tensors"? I have never heard of such a thing. The way I remember it from my physics studies is that "scalar products" is a synonym for "inner products" which (for non-complex fields) are positive definite, symmetric, bilinear forms. That means, they map two vectors to a positive scalar value, not tensors of other rank. I think I heard people use the term inner product between tensors to talk about what I know as contractions, hence my second question. However, I really don't know what a scalar product of two tensors should then be.
PatrickKa commented 1 month ago

Ad 3. Seems like I was a bit off with my terminology. According to Wikipedia, "scalar product" is a synonym for "dot product". This does not help my confusinon though, since the dot product is just a special case of an inner product.

mpusz commented 1 month ago

Should tensors of arbitrary rank be supported, or will the rank be capped, since most "normal" physics does not require more than rank-4 tensors? (Even those are rather rare. The only two rank-4 tensors I could think of off the top of my head were the elasticity tensor and the Riemann curvature tensor)

Yes, you are right. ISO 80000 part 2 only provides operations on tensors of the second and fourth order.

mpusz commented 1 month ago

When you speak of inner products between tensors or a tensor and a vector, do you actually mean contractions? So, e.g., if a_ij, and b_jk denote the components of two rank-2 tensors is (a ⋅ b)_ik = Σ_j a_ij b_jk?

Yes, it is how ISO 80000-2 defines them. Please contact me via email, and I will send you more details if needed.

mpusz commented 1 month ago

What do you mean with "scalar product of two tensors"?

ISO 80000-2 defines a scalar product of two tensors T and S of the second order as $T : S = \sum_i\sumjT{ij} S_{ij}$ and the result is a scalar quantity.

PatrickKa commented 1 month ago

Thank you for clarifying. It seems like I should have just read ISO 80000-2 😅. I won't have time in the foreseeable future to work on this issue, but if it is still open when I have less on my plate, I'd be happy to help with it.

Also, I didn't know that I can use $\LaTeX$ here. That's cool.

chiphogg commented 1 month ago

I've got a couple of questions/comments on this issue.

To not depend on external libraries and to provide correct, fast, and simple examples of LA usage with quantities we need simple vector and tensor types (up to 3 dimensions) together with the following operations:

Why artificially limit to 3 dimensions?

  • a % b - modulo where both arguments should be of the same quantity kind and character

What does it even mean to take the "modulo" of vectors or tensors?

  • a × b - cross product of two vectors

We could do this, for 3 dimensions only. But the cross product operation is a hacky kludge, mathematically. It's just a disguised wedge product, which would be valid for arbitrary dimensions.

If you're giving different priorities to these operations, then I'd put this one in a lower priority tier.

mpusz commented 1 month ago

Why artificially limit to 3 dimensions?

My understanding is that ISQ talks only about Cartesian (orthonormal) coordinates in ordinary space. This even means that we do not have to provide 1- and 2-dimensional vectors to support ISQ. I proposed to support them because it should be quite easy and would allow us to model specific scenarios for some quantities.

What does it even mean to take the "modulo" of vectors or tensors?

You are right; this operation probably makes sense for scalars only.

We could do this, for 3 dimensions only.

And this is why I said that we need 3-dimensional vectors and matrices/tensors.

But the cross product operation is a hacky kludge, mathematically. It's just a disguised wedge product, which would be valid for arbitrary dimensions.

Agree, but this is how ISQ defines its quantities. And we need to support it as well as all others with the same priority.

mpusz commented 1 month ago

To be clear, I do not want to invent a new full-blown LA library here. We need two simple types working on three components (possibly also less) and a bunch of operations on them. That is it. If someone wants to do more advanced stuff, there are always Eigen, Blaze, and other beautiful libraries out there. We will also provide some interop examples with them, if possible.

mpusz commented 1 month ago

Those types are only meant to be used as representation types of vectors and tensor quantities in this library and will mostly be useful for our unit testing and examples. We may also expose them in the library interface to our users, but we will not force their usage. Alternatives should also be easy to use.

Those types will not be used to model huge vectors or matrices of quantities (even if homogeneous). For that, we already have many LA libraries.

chiphogg commented 1 month ago

If we're going that route, I suggest not even making the dimensionality configurable at all. Just keep it simple and use 3D vectors.

Spammed commented 1 month ago

I'm used to handle 6 degree of freedom (3 distances & 3 orientaton angles) in a 3D carthesian world, hence I'm often use 6D-vectors. E.G. for 'locations' ([m],[m],[m],[rad],[rad],[rad]) and for 'forces' ([N],[N],[N],[Nm],[Nm],[Nm]) etc. At present without any unit control. I realize that this can also be set up differently, e.g. as a struct of 2 3D vectors. But it worries me to have to do this in order to benefit from your great work. Don't let that slow you down, I just wanted to report...

chiphogg commented 1 month ago

@Spammed, a full "linear algebra + quantities" solution would be able to handle this. We would need interfaces to provide a separate unit for each index, and facilities for correctly combining them in matrix multiplication.

And it's known how to do this --- but there just isn't an open source library that provides this yet. Daniel Withopf began presenting one solution publicly back in 2021, but it's a closed source library that's internal to Bosch. Also in 2021, I implemented a proof of concept for wrapping Eigen with this kind of interface, and showed zero runtime performance loss. I presented this idea in this segment of my CppCon talk, but I didn't get a chance to implement actual usable end-user interfaces, so this still isn't implemented in Au yet (see aurora-opensource/au#70 to track progress, but know that it will be a while before we start).

All of which is to say: what you want is possible and desirable, even with a single vector (not a struct of 2 3D vectors), but the stopgap solution in this present issue won't give it to you, and I don't know of any open source solutions that are planned with a known timeframe (this is planned for Au, but there is no timeframe, and we surely won't start on it in 2024).

Spammed commented 1 month ago

Thank you very much for your response and for addressing and summarizing the problem. I greatly appreciate the work of all of you and know that your joint deliberations are leading step by step to a better future.