mpusz/mp-units

Consider adding expression templates support for `norm`, `dot`, and `cross` products on vector quantities

mpusz opened this issue · 4 comments

mpusz commented

We all expect that length * time will not result in speed because we clearly have an error in calculation.

Dot and cross products in vector quantities domain are really similar operations to * and / in the scalar domain. Using a wrong operation by mistake results in a different quantity than expected. Not having support for it makes the library not safe for projects that work on vector quantities.

As the V2 library has now experimental support for vector and tensor quantities as well, it could be a good idea to add the above to improve safety. Also, it would improve ISQ system specification as we will not have to override the quantity character for derived quantities, i.e.:

QUANTITY_SPEC(power, force * velocity, quantity_character::scalar);
QUANTITY_SPEC(moment_of_force, position_vector * force);  // vector
QUANTITY_SPEC(torque, moment_of_force, quantity_character::scalar);
QUANTITY_SPEC(surface_tension, force / length, quantity_character::scalar);  // TODO what is a correct equation here?

The above could be rewritten as:

QUANTITY_SPEC(power, dot(force, velocity));   // or we can assume that `op*` means a dot product
QUANTITY_SPEC(moment_of_force, cross(position_vector, force));
QUANTITY_SPEC(torque, dot(moment_of_force, unit_vector));  // I do not know what a `unit_vector` here should be
QUANTITY_SPEC(surface_tension, norm(force) / length);

Please do share your ideas or comments.

mpusz commented

How to support things like 42 * (m / s) * vector{2, 1, 0}? The result should compile, not be assignable to speed and be convertible to velocity.

I think my off the cuff reaction is that I would personally find it hard to be confident I was getting this right.

Take this with a grain of salt, since I'm also not at all caught up on how vector and tensor "characters" work. 🙂 But let's take the cross product as an example. The vector cross product is a kludge that only makes sense in 3 dimensions. It's the dual of the wedge product, a more general and better-founded product which works in arbitrary dimensions. (One of my favorite books on this is Geometric Algebra for Computer Science. This gets discussed in chapter 3.)

With the cross product, you end up with something that "kinda" acts like a vector for some operations, but not for others. For example, under inversion, "real" vectors are negated, but cross products aren't. I think the usual way to phrase this is to call the cross product a "pseudo-vector". So, besides the "vector" character, should we have a "pseudo-vector" character? And if so, how about pseudo-scalars and pseudo-tensors?

I do see some appeal here, some benefit to figuring it all out. For example, based on the transformation properties, it seems clear that no operation which adds vectors and pseudo-vectors can be meaningful... right? (Does anyone have a counter-example?) If so, then it'd be great if we could catch this at compile time! But then again, many (most?) users may not even know when they have a pseudo-vector instead of a "real" vector, so this could lead to frustrating compile time errors. (For example, normal vectors are pseudo-vectors, but users are very likely to reach for "regular" vectors for these.)

tl;dr: I know just enough to recognize that there's significant complexity lurking here, but not enough to know how to handle it.

mpusz commented

Well, yes, it seems that is tough indeed, but said that the physical quantities library is easy 😨. I thought about this for a while, and I am not sure if there is any other solution if we want to do it right...

Let's assume a user that tries to do it correctly, and not just to provide slideware code that forgives everything. By this, I mean a user that uses a break linear algebra library to express vector quantities like position_vector and force. This user wants to calculate a moment_of_force. All the existing physical libraries on the market define only multiplication and division operators. Both operators are wrong, but as there is no other way to do it, the user chooses multiplication. All the libraries are fine to create a type of quantity of moment of force from it, but then try to do the multiplication on the underlying vector representation type and crash. As a result, the user has to provide custom overloads for operator* for such quantities that will do a vector product on those for the underlying type. Providing those happens to not be that easy as none of the libraries mark such quantities to have a vector character, so a generic code for them can't be provided.

I might be wrong, but I think that none of the libraries on the market actually has unit tests involving real linear algebra as underlying types for vector and tensor quantities. That is why we do not see bugs for that, and users suffer or do not know they can do any better. Writing unit tests like:

1 * force * (1 * length) == 1 * moment_of_force

is not enough and is simply wrong. But this exactly is what we have unfortunately done for years 😢

mpusz commented

It seems that Pint supports those operations already: https://pint.readthedocs.io/en/stable/user/numpy.html#Function/Method-Support.