Operations on vector quantities

mpusz commented 4 days ago

I am in the process of improving the support for vector quantities. We have had many simplifications and errors in this subject so far, and we should try to improve here. The question here is how strict we should be.

Let's see the following example.

We were defining speed and velocity as

inline constexpr speed final : quantity_spec<length / time> {} speed;
inline constexpr velocity final : quantity_spec<speed, position_vector / duration> {} velocity;

which was a simplification, wrong, and inconsistent with ISQ definitions. ISQ defines them in the following way:

inline constexpr velocity final : quantity_spec<displacement / duration> {} velocity;
inline constexpr speed final : quantity_spec<norm(velocity)> speed {};

_Well, velocity is actually defined in terms of delta position_vector, but that is a displacement and should be fixed in our definitions._

Note that speed and velocity are not in the same tree, and the real definition order (dependency) is reversed to what we had before.

Moreover, our flagship example for many years looked like this:

constexpr QuantityOf<isq::speed> auto avg_speed(QuantityOf<isq::length> auto d,
                                                QuantityOf<isq::time> auto t)
{
  return d / t;
}

quantity height = isq::height(20 * m);
quantity time = 5 * s;
quantity speed = avg_speed(height, time);

This is a simplification as length does not necessarily have to describe a displacement. What kind of speed will I get by dividing the height of the building by time? The building does not move or change its height. But this is how we are used to do it in the engineering.

Of course, if I drop a rock from the top of the roof and observe it falling down, I may use the height of this building to measure the displacement of the rock on its way down.

Also, a displacement is a vector quantity and should be described by the vector representation type. However, when we want to get speed and not velocity, we should not care.

To summarize, a proper function template should look like:

constexpr QuantityOf<isq::speed> auto avg_speed(QuantityOf<norm(isq::displacement)> auto d,
                                                QuantityOf<isq::duration> auto t)
{
  return d / t;
}

and we should call it with:

quantity speed = avg_speed(quantity_cast<norm(isq::displacement)>(height), time);

The question is, is that not too much? Does requiring a quantity_cast to a norm(displacement) not sacrifice usability for physical safety/correctness too much? What can we do to improve that?

JohelEGP commented 4 days ago

The first signature still looks good. It's still a valid formula of speed. The onus is then on the author of the function to cast the length to whatever QuantityOf<isq::speed> requires.

mpusz commented 4 days ago

Well, we could argue if that is a valid formula. In the link you provided, it is explicitly stated that d is a "distance covered" which is different from the static height of the building 😉

But yes, we could implement this as well as:

constexpr QuantityOf<isq::speed> auto avg_speed(QuantityOf<isq::length> auto d,
                                                QuantityOf<isq::time> auto t)
{
  return quantity_cast<norm(isq::displacement)>(d) / t;
}

This makes it easy to use but less correct in the physical sense.

mpusz commented 4 days ago

I would like to mention that without introducing norm, vector_product, and scalar_product, it is nearly impossible to make any vector quantities like (velocity, acceleration, force, moment of force) work correctly and resemble their physical properties.

mpusz commented 4 days ago

We also probably should refactor our tree of quantities of kind length to something like:

flowchart TD
    length["<b>length</b><br>[m]"]
    length --- width["<b>width</b> / <b>breadth</b>"]
    length --- height["<b>height</b> / <b>depth</b> / <b>altitude</b>"]
    width --- thickness["<b>thickness</b>"]
    width --- diameter["<b>diameter</b>"]
    width --- radius["<b>radius</b>"]
    radius --- radius_of_curvature["<b>radius_of_curvature</b>"]
    length --- path_length["<b>path_length</b>"]
    path_length --- distance["<b>distance</b>"]
    distance --- displacement["<b>displacement</b><br>{vector}"]
    distance --- radial_distance["<b>radial_distance</b>"]
    radial_distance --- position_vector["<b>position_vector</b><br>{vector}"]
    length --- wavelength["<b>wavelength</b>"]

The primary changes are:

norm(displacement) is a distance, but not every distance is a norm(displacement). Thanks to this, if we will have a distance quantity, we will be able to call avg_speed(norm(isq::displacement)(distance), time) with no need for quantity_cast.
norm(position_vector) could be considered as radial_distance. Please let me know if you think this is incorrect. We can move position_vector directly under distance, but that will require a quantity_cast to convert norm(position_vector) to radial_distance or vice versa.
Not that long ago, I thought that position_vector should be a leaf from displacement to express that position_vector can be modeled as a displacement from the origin. However, position_vector does not imply relocation, which is inherent to displacement.

Spammed commented 4 days ago

Perhaps each vectorial quantity should be named '..._vector'. So 'displacement_vector' instead of 'displacement' and 'velocity_vector' instead of 'velocity'? 'speed = norm(velocity_vector)' is perhaps more natural than 'speed = norm(velocity)'.

mpusz commented 3 days ago

We can't rename ISQ quantities. Their names are standardized.

A user can always provide an alternative name with:

constexpr auto velocity_vector = velocity;

Also, a helper for norm(isq::displacement) can be provided:

constexpr auto displ_norm = norm(isq::displacement);

With the above the avg_speed could be called as:

avg_speed(displ_norm(distance), time);

mpusz commented 3 days ago

Looking at the tree above and discussing the problem with automotive experts in the ISO C++ Committee who know mp-units really well already, we decided to redefine speed to:

inline constexpr velocity final : quantity_spec<displacement / duration> {} velocity;
inline constexpr speed final : quantity_spec<path_length / duration> {} speed;

It is inconsistent with ISQ but is more correct and useful. velocity is always a temporary speed vector attached to the path we travel. The magnitude of it gives us a temporary speed. However, we often need to calculate the speed of the entire path. This is why I believe that path_length / duration is a better definition here.

As a result, we still have to quantity_cast from height, but that is exactly why we have a quantity_cast in the library. If someone wants to calculate speed, then a path_length or distance should be provided.

To summarize, our new avg_speed definition may be:

constexpr QuantityOf<isq::speed> auto avg_speed(QuantityOf<isq::path_length> auto d,
                                                QuantityOf<isq::duration> auto t)
{
  return d / t;
}

norm(isq::displacement) can still be provided if needed as it is implicitly convertible to isq::path_length.

mpusz commented 3 days ago

Now the question is, does passing a norm(isq::position_vector) make any physical sense? It will compile, but it probably is not a correct behavior.

Despite radial_distance being defined as :

distance (item 3‑1.8), where one point is located on an axis or within a closed non self-intersecting curve or surface

Maybe it should not be in the same but parallel branch to distance. Maybe radial_distance is actually not a distance? Or maybe it is, but not every distance and path_length make sense as an argument to speed? distance does not always imply a relocation of the object, which is needed to talk about speed. I am lost again now 😢

Any thoughts...

Spammed commented 3 days ago

What exactly does the arrow in the tree represent? If an arrow means 'The lower is a specialized form of the upper'? A specialization, I think, does not change the character and therefore scalars and vectors must be kept in separate trees. Oh dear!

mpusz commented 3 days ago

Yes, it means exactly that.

Separate tree is not an option. Displacement and position vector are still quantities of length and they are measured in metres.

I considered putting them in a separate branch (this is exactly what we have until 2.4.0) but I think it's wrong as well. There are at least two reasons for that:

I think that saying that a 3D length is a more specialized version of length is a perfectly fine statement. This applies to position_vector being a 3D version of radial_distance, and displacement being a 3D version of distance travelled. Of course vectors can also have other dimensionality than 3D.
magnitude(displacement) is distance_traveled which is a more specialized version of distance. Even though in most cases we can't directly convert a 3D displacement to distance because of representation type incompatibility, the magnitude of such displacement can be perfectly converted to distance. Which is exactly what we want to achieve here.

mpusz commented 3 days ago

How about:

constexpr auto distance_traveled = norm(isq::displacement);

constexpr QuantityOf<isq::speed> auto avg_speed(QuantityOf<distance_traveled> auto d, QuantityOf<isq::time> auto t)
{
  return d / t;
}

?

We can also consider adding distance_traveled to the length hierarchy between distance and displacement.

Spammed commented 3 days ago

The term 'distance_traveled' reminds me a lot of 'path_length'. Why is 'distance' a specialization of 'path_length'? Or what is 'path_length' good for?

Can we delete the old 'path_length' and introduce a new 'path_length' as a specialization of 'displacement'? (ISQ probably contradicts this)

I think I remain confused by "Specialization can change the 'character'".

mpusz commented 3 days ago

Let's describe the elements in the branch here:

path_length - length of a rectifiable curve between two of its points
distance - shortest path length between two points in a metric space (does not imply repositioning). It does not have to be equivalent to the magnitude of the displacement as you reposition from New York to Tokyo and certainly you did not travel the shortest straight line through the globe of Earth :wink:
radial_distance - distance, where one point is located on an axis or within a closed non self-intersecting curve or surface
displacement - vector describing repositioning the object between two points in space (perfect for velocity)
position_vector - vector from the origin of the coordinate system to the point in space (not suitable for the velocity)

A magnitude/norm of position_vector is a radial_distance. Every radial_distance is distance. I claim that speed should not be obtained from a radial_distance as it does not make much sense (unless you actually traveled from the center of the circle to its border). This means that not every distance implies repositioning a body over this distance. When I stand in the middle of the forest, there are certain distances from me to nearby trees. If I stand there for 10 seconds, does it mean that I have multiple speeds at the same time? 😉

path_length is even more generic than distance. It does not imply the shortest route between the points.

dwith-ts commented 3 days ago

Just a side note as it has been mentioned twice already: taking the norm of a position vector should / can not be a valid operation. You can only take the norm of a displacement vector not of a single point / position vector.

mpusz commented 3 days ago

@dwith-ts, thanks for feedback. I was not aware of this. This complicates things a bit, but I think I know how to make it work.

Do you also suggest that we limit position_vector quantity_spec type to quantity_point only, on which we can't do any arithmetics anyway? Said otherwise, quantity<isq::position_vector[m]> should be a compile-time error?

Spammed commented 3 days ago

Many thanks for the patient discussion. With '_Every radialdistance [vector] is [a scalar] distance.' (or 'Every [position_vector is a scalar] _radialdistance is [a scalar] distance.' ?) I have to get along somehow.

mpusz commented 3 days ago

Just a side note as it has been mentioned twice already: taking the norm of a position vector should / can not be a valid operation. You can only take the norm of a displacement vector not of a single point / position vector.

@dwith-ts, is this true even for polar or spherical coordinate systems? I thought that the norm of a position vector would be a radial distance in such cases.

mpusz commented 3 days ago

With 'Every radial_distance [vector] is [a scalar] distance.' (or 'Every [position_vector is a scalar] radial_distance is [a scalar] distance.' ?) I have to get along somehow.

radial_distance is defined as a scalar quantity. Every radial_distance is a distance. position_vector is a vector quantity.

dwith-ts commented 3 days ago

Just a side note as it has been mentioned twice already: taking the norm of a position vector should / can not be a valid operation. You can only take the norm of a displacement vector not of a single point / position vector.

@dwith-ts, is this true even for polar or spherical coordinate systems? I thought that the norm of a position vector would be a radial distance in such cases.

Well, with both Cartesian and polar / spherical coordinate systems you could decide to use a position vector interchangeably with its displacement from the origin, then allowing the norm operation could be justified.

However, in math terms a norm is always defined on a vector space and the definition of such a vector space contains (mathematical) points and vectors. The norm operation can then be applied only to such vectors and is independent of the chosen origin.

I don’t know all the library details so I can’t say much about whether a quantity<isq::position_vector[m]> should be a compile time error. I think there is a more general issue that there are now 3 places that can be used to encode the affine space properties / distinction between points and displacement vectors:

quantity vs quantity_point
position_vector[m] vs displacement_vector[m]
The dimension also encodes some sort of information about the affine space as [m/s] is always an indication that we have a quantity and not a quantity_point; => position_vector[m/s] and displacement_vector[m/s] are interchangeable me as the velocity information inside the vector is a quantity and not a quantity_point (velocities are always relative to something)

And now the question is which of the combinations should be allowed / are meaningful:

quantity<position_vector[m]>
quantity_point<position_vector[m]>
quantity<displacement_vector[m]>
quantity_point<displacement_vector[m]>
quantity<position_vector[m/s]>
quantity_point<position_vector[m/s]>
quantity<displacement_vector[m/s]>
quantity_point<displacement_vector[m/s]>

I hope I got the library-related terminology right, if not feel free to clarify.

JohelEGP commented 3 days ago

displacement - vector describing repositioning the object between two points in space (perfect for velocity)

Is this definition from an ISO draft?

mpusz commented 2 days ago

Is this definition from an ISO draft?

No, I just improvised a bit :wink:

The actual definitions are:

`displacement`

Definition: vector quantity between any two points in space
Remarks: Displacement vectors are so-called free vectors, i.e. their magnitude and direction do not depend on a particular coordinate system. The magnitude of this vector is also called displacement.
Usage examples: shear_strain

`position_vector`

Definition: vector quantity from the origin of a coordinate system to a point in space
Remarks: Position vectors are so-called bounded vectors, i.e. their magnitude and direction depend on the particular coordinate system used
Usage examples: moment_of_force, angular_momentum

Additionally:

speed and electric_dipole_moment use subtraction of two position_vectors in their results. Subtracting two position vectors gives us a displacement
mechanical_work, magnetic_tension, magnetomotive_force, and linked_flux integrate over position_vector. I am pretty sure that it can also be replaced with a displacement for our needs?

mpusz commented 2 days ago

@dwith-ts, I am afraid that, according to the above definitions, we can't disallow taking the magnitude of a position vector.

mpusz commented 2 days ago

The dimension also encodes some sort of information about the affine space as [m/s] is always an indication that we have a quantity and not a quantity_point; => position_vector[m/s] and displacement_vector[m/s] are interchangeable me as the velocity information inside the vector is a quantity and not a quantity_point (velocities are always relative to something)

I am unsure if I agree with that. If I drive a car that has a speed-o-meter and some G-sensors, I can make a measurement/readout from those instruments. I think that it is OK to model this with a quantity_point<isq::velocity[km/h]>. Adding two such results does not make sense. If I do the same after a few seconds, I can subtract two of them to get quantity<isq::velocity[km/h]> to find out how my velocity changes in time. Does it make sense?

quantity<position_vector[m/s]> quantity_point<position_vector[m/s]> quantity<displacement_vector[m/s]> quantity_point<displacement_vector[m/s]>

According to the ISQ, position_vector and displacement are quantities of length, so we can't express them in m/s.

dwith-ts commented 2 days ago

I am unsure if I agree with that. If I drive a car that has a speed-o-meter and some G-sensors, I can make a measurement/readout from those instruments. I think that it is OK to model this with a quantity_point<isq::velocity[km/h]>. Adding two such results does not make sense. If I do the same after a few seconds, I can subtract two of them to get quantity<isq::velocity[km/h]> to find out how my velocity changes in time. Does it make sense?

Yes, you raise an excellent point here: In the most general form, we would need what is called a 3-frame notation to describe the full problem (https://rpg.ifi.uzh.ch/docs/teaching/2024/FurgaleTutorial.pdf, not sure if I already linked that). The 3 frames indicate

what we describe (e.g. "position of coordinate frame C", "velocity of coordinate frame C"
relative to what and (e.g. "relative to frame B")
in which frame / using which basis vectors this relation is expressed (e.g. "in the basis vectors of frame A")

Now your example correctly points out that we cannot add the velocity of frame C (your car in the example) relative to the earth at time t_1 to the velocity of the same frame at time t_2, even if they are expressed in the same basis vectors. However, we can add the velocity of frame D (a second car) relative to frame C (the first car) to the velocity of frame C relative to the earth which would give us the velocity of frame D relative to the earth (this only works if both use the same basis vectors for the representation).

Now we could make the same point for displacement vectors: it is also not very meaningful to add the displacement between points / frames C and D to the displacement between A and B. This addition is only meaningful if we would add the displacement of C relative to B and the displacement of B relative to A which gives us the displacement of C relative to A (again assuming that the unmentioned 3rd frame describing the basis vectors is identical).

So both velocities and displacement have a similar issue here that they cannot be added in the general case, only in special cases. However, it seems that this distinction is more obvious and better known for velocities than for displacements.

So what does this help us? I think that we most likely will have to make certain approximations / simplifications of the full problem space, the only question is which simplification exactly.

What I have done in my library:

there is a full distinction between position and displacement vectors with the usual rules (e.g. 2 positions cannot be added)
for velocity vectors I also distinguish between "points" and "delta-vectors" with the same rules, but the general recommendation is to just use a delta-vector if the vector contains only velocities

dwith-ts commented 2 days ago

According to the ISQ, position_vector and displacement are quantities of length, so we can't express them in m/s.

Is ISO here maybe putting the distinction between quantity and quantity_point into the names position_vector and displacement as ISO does not have the distinction between quantity and quantity_point? It feels like we should have something like quantity<3d-length> and quantity_point<3d-length> instead of having both quantity / quantity_point and position_vector / displacement which gives us the four options

quantity<position_vector> *
quantity<displacement>
quantity_point<position_vector>
quantity_point<displacement> *

where the ones marked with * don't seem meaningful to me because they are somehow self-contradictory.

mpusz commented 2 days ago

In the entire ISO/IEC 80000 series that defines ISQ, you will not find any mention of the affine space abstractions, points, or deltas. I mentioned it even in the blog post. Also, as mentioned in #647, even if we can sense some point/delta properties they are they by accident and are inconsistent.

It feels like we should have something like quantity<3d-length> and quantity_point<3d-length> instead of having both quantity / quantity_point and position_vector / displacement which gives us the four options

This is exactly why I ask if we should limit:

position_vector, altitude, time to quantity_point
displacement, height, duration to quantity

Then, it should be clearer.

I am right now working to allow the following happen automatically:

quantity_point<isq::position_vector[m]> qp1(100 * m); 
quantity_point<isq::position_vector[m]> qp2(50 * m);
auto q = qp1 - qp2;   // results with quantity<isq::displacement[m]>

The same would work for altitude and time if we define them the same way in our ISQ definition.

Please note that as of today, also the following would work:

quantity<isq::position_vector[m]> q1(100 * m); 
quantity<isq::position_vector[m]> q2(50 * m);
auto q = q1 - q2;   // results with quantity<isq::displacement[m]>

but as you wrote already, it contradicts itself. Also, using quantity_point<isq::displacement[m]> might not have sense.

However, as I wrote in #646:

However, we can't do arithmetic on quantity_point. If there is an operation that requires a point (e.g., usage of position_vector in moment_of_force), we would need to do qp.quantity_from_zero(). However, this should also return a quantity of displacement (and not position_vector), which will violate the requirements of the equation of moment_of_force.

dwith-ts commented 2 days ago

It feels like we should have something like quantity<3d-length> and quantity_point<3d-length> instead of having both quantity / quantity_point and position_vector / displacement which gives us the four options

This is exactly why I ask if we should limit:
* `position_vector`, `altitude`, `time` to `quantity_point`

* `displacement`, `height`, `duration` to `quantity`

The limiting suggestion you made is still a bit different from what I'd imagine ideally.

For a clean-sheet design I'd go with either something like inversing the relation by using vector<quantity[m]> and vector<quantity_point[m]> which would also be future-proof for non-uniform vectors. But I presume that this would be a big switch away from the current design so I'll focus on the 2nd option that is closer to the current design:

We could use aliases for both position_vector and displacement

using position_vector = quantity_point<vector[m]>
using displacement = quantity<vector[m]>

This also has the advantage that we can do the exact same thing for velocity vectors with

using velocity_point = quantity_point<vector[m/s]>
using velocity = quantity<vector[m/s]>

Is this something that has some fundamental contradiction to ISO or could this be used?

If all of the above is not possible, then I'd say the option you outlined (restricting what can be put inside quantity / quantity_point) is something worth considering as it is slightly better that what is currently implemented in the library. On the other hand, this could also be a source of confusion for users as it makes usage a bit more difficult.

mpusz commented 2 days ago

For a clean-sheet design I'd go with either something like inversing the relation by using vector<quantity[m]> and vector<quantity_point[m]> which would also be future-proof for non-uniform vectors. But I presume that this would be a big switch away from the current design so I'll focus on the 2nd option that is closer to the current design:

With such a design, we have several issues:

mp-units will not know if the outer LA library is performing dot or cross product, so it will be unable to determine and validate the resulting quantity and conversions
allowing the above would also require the isq::position_vector(42 * m) * isq::force(10 * N) to work which is bad,
requires us to type more e.g., vector{1 * m, 2 * m, 3 * m} instead of vector{1, 2, 3} * m
raises plenty of questions about what to do in the case of things like vector{isq::height(1 * m), 2.13 * m, 3 * mm} or vector{1 * m, 1 * s, 1 * N}

mpusz commented 2 days ago

We could use aliases for both position_vector and displacement

using position_vector = quantity_point<vector[m]> using displacement = quantity<vector[m]>

Yes, if ISO/IEC 80000 would be aware of affine spaces as a primary design construct, we could use:

using displacement = quantity<m, vector<3, double>>;
using position_vector = quantity_point<m, default_point_origin(kind_of<isq::length>), vector<3, double>>;

or

using displacement = quantity<isq::displacement[m], vector<3, double>>;
using position_vector = quantity_point<isq::position_vector[m], default_point_origin(isq::position_vector), vector<3, double>>;

Yes, I know, the second line in both cases is inconvenient to type for a custom type 😢

mpusz commented 2 days ago

This also has the advantage that we can do the exact same thing for velocity vectors with

using velocity_point = quantity_point<vector[m/s]> using velocity = quantity<vector[m/s]>

We do it now as follows:

using velocity = quantity<m / s, vector<3, double>>;
using position_vector = quantity_point<m / s, default_point_origin(kind_of<isq::length / isq::time>), vector<3, double>>;

or

using velocity = quantity<isq::velocity[m/s], vector<3, double>>;
using position_vector = quantity_point<isq::velocity[m/s], default_point_origin(isq::velocity), vector<3, double>>;

So this is also compatible in both cases.

mpusz commented 2 days ago

Is this something that has some fundamental contradiction to ISO or could this be used?

No contradictions. It is exactly how it may be used now. My question is if we should prevent creating clearly "point-like" ISQ quantity types (e.g., position_vector, altitude) with quantity?

rothmichaels commented 2 days ago

Not the main point of this thread (still reviewing the rest) but one question I had while looking at the first example in the original post: why are we spelling this function norm vs abs or magnitude?

mpusz commented 2 days ago

That is a very good question. For now, norm is just a proposal, and I am open to alternatives. I discussed this a lot with some LA experts in ISO lately:

abs() is probably the least correct, but it allows generic programming.
magnitude() conflicts with unit magnitude, but in order to prepare grounds for such usage, I renamed magnitude to unit_magnitude and Magnitude to UnitMagnitude in the last days. Note that we still have mag, mag_ratio, and mag_power helpers. We can also rename them to unit_XXX but I did not do it. Thanks to this, we are now able to introduce it.
norm(), as I've heard, is the most correct term for vectors. It is also less to type than magnitude.

Does anyone have a good overview of LA libraries on the market? Which one is the most common for vectors? We want to use whatever is the most popular to enable generic programming.

Anyway, our CPO accounts for the most common options already:

https://github.com/mpusz/mp-units/blob/bda112418a6bff3fd58e435e4e5fafefbcbf7a58/src/core/include/mp-units/framework/customization_points.h#L115-L129

mpusz / mp-units

Operations on vector quantities #646

`displacement`

`position_vector`