Use consistent epsilon values for testing inertia and reporting errors.

[mitiguy]

The epsilon used to test for an invalid inertia differs from the epsilon = 0 used to report an error (in the error message).

---

This change is 

[mitiguy]

Feature review +@sherm1

Reviewable status: LGTM missing from assignee sherm1(platform), needs at least two assigned reviewers

[sherm1]

Feature with minor comments
+@SeanCurtis-TRI for (Thursday) platform review per rotation, please

Reviewed 2 of 2 files at r1, all commit messages.
Reviewable status: 2 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform)

---

multibody/tree/rotational_inertia.h line 944 at r1 (raw file):

  }

  // Tests whether each moment of inertia is non-negative (to within `epsilon`)

minor: this comment and the note below still refer to `epsilon` (in several places) as though it was still a parameter. If you want to discuss the tolerance here, consider introducing a mathematical ε instead.

---

multibody/tree/rotational_inertia.h line 948 at r1 (raw file):

  // triangle-inequality test requires `epsilon` when the sum of two moments are
  // nearly equal to the third one. Example: Ixx = Iyy = 50, Izz = 100.00000001,
  // or Ixx = -0.0001 (negative),  Ixx = 49.9999,  Iyy = 50.

BTW this Ixx = -0.0001, Ixx = 49.9999 comment doesn't make sense to me

[sherm1]

Reviewed 1 of 1 files at r2, all commit messages.
Reviewable status: 2 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @mitiguy)

[mitiguy]

Reviewable status: LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits)

---

multibody/tree/rotational_inertia.h line 944 at r1 (raw file):

Previously, sherm1 (Michael Sherman) wrote…

minor: this comment and the note below still refer to `epsilon` (in several places) as though it was still a parameter. If you want to discuss the tolerance here, consider introducing a mathematical ε instead.

Done .

---

multibody/tree/rotational_inertia.h line 948 at r1 (raw file):

Previously, sherm1 (Michael Sherman) wrote…

BTW this Ixx = -0.0001, Ixx = 49.9999 comment doesn't make sense to me

Done. Removed.

Note: I consolidated, simplified, and shortened the comments about epsilon.

[SeanCurtis-TRI]

Good morning. First pass done.

Unfortunately, while addressing the inconsistent epsilons, this has revealed a number of other inconsistencies in this corner of the code.

Reviewed 1 of 2 files at r1, 1 of 1 files at r2, all commit messages.
Reviewable status: 12 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @mitiguy)

---

a discussion (no related file):
No tests?

---

multibody/tree/rotational_inertia.h line 950 at r2 (raw file):
BTW "is typically much smaller than the largest possible principal moment of inertia"

I have two issues with this:

• The word "typically" hides a lot of sins. It doesn't say "always" or simply say that epsilon scales with the largest principal moment. This document implies that describing epsilon is important here, but then hiding the "atypical" case seems counter to that purpose. What is being hidden here is that otherwise valid inertias whose principal moments of inertia are simply "too small" also are deemed invalid. This is a typical Drake-ism knowing that we're targeting human-scale dynamical systems, but it might be worth being more clear about it.
• When I first read this phrase, my brain went in the wrong direction and I misinterpreted it. I believe a slight rephrasing would help people like me and not offend your sensibilities (although, if it does, you can skip the suggested text):

Because the re-expression error scales with the largest principal moment, the epsilon also scales (but relative to an upper-bound estimate of the largest principal moment). We also exclude rotational inertias whose overall scale are considered too small relative to Drake's expected robot workspace scale.

---

multibody/tree/rotational_inertia.h line 592 at r2 (raw file):

    // Calculate principal moments of inertia p and then test these principal
    // moments to be mostly non-negative and also satisfy triangle inequality.
    const Vector3<double> p = CalcPrincipalMomentsOfInertia();

nit: Again, if you keep the |tr/2| proxy, you shouldn't be computing p, you should simply be passing the generic Ixx, Iyy, and Izz into AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality() and not waste the work of computing the principal moments for this test.

---

multibody/tree/rotational_inertia.h line 948 at r2 (raw file):

  // The triangle-inequality test requires ε when the sum of two moments are
  // nearly equal to the third one. Example: Ixx = Iyy = 50, Izz = 100.00000001.
  // The positive (non-zero) ε accounts for round-off errors, e.g., from

BTW Every positive number is non-zero, making "(non-zero)" a redundant elaboration.

Code quote:

 (non-zero)

---

multibody/tree/rotational_inertia.h line 950 at r2 (raw file):

  // The positive (non-zero) ε accounts for round-off errors, e.g., from
  // re-expressing inertia in another frame and is typically much smaller than
  // the largest possible principal moment of inertia in the rotational inertia.

BTW I'm not blocking on this, but I'll reiterate the following:

The justification of using epsilon (i.e., the original inertia satisfied the constraint, but the re-expressed values don't because of rounding error) is perfectly valid. This increases the space of "valid" inertia values to include those cases. Yay!

Remember, for every such inertia, there is another inertia that didn't satisfy the triangle inequality constraint but, due to re-expression's rounding error, it now falls within the epsilon, thus it gets incorrectly classified as valid.

It might be worth acknowledging that fact. That we'd rather include some ever-so-slightly invalid inertias rather than exclude the ones that were definitely valid but were penalized by finite precision during valid mathematical operations.

I just hate the idea of thinking of "epsilon" as some amazing panacea, especially when used like this where you're partitioning the value space into two regions instead of three.

---

multibody/tree/rotational_inertia.h line 955 at r2 (raw file):

  static boolean<T> AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality(
      const T& Ixx, const T& Iyy, const T& Izz) {
    // Denoting Imin and Imax as the smallest and largest possible moments of

nit: You've defined Imax as the largest "possible" moment of inertia. And then you've asserted that tr/3 ≤ Imax ≤ tr/2.

This is an error. The true statement is that Imax = tr/2.

In other words, you've blurred whether Imax refers to the actual largest moment, or the possible largest moment. Your English says "possible", your math says "actual".

Now, if you'd defined Imax as the actual largest moment, then your inequality would be true.

See my follow up note on this documentation.

Code quote:

possible

---

multibody/tree/rotational_inertia.h line 963 at r2 (raw file):

    // related to machine precision multiplied by Imax = tr/2.
    // Tangent: The product of inertia Iᵢⱼ with the largest absolute value is at
    // most half the largest moment of inertia, i.e., |Iᵢⱼ| ≤ Imax/2.

BTW You've flooded this comment with pedagogical trivia. As I read the comment and the code, you've made many assertions that don't directly reflect what is in the code or even necessary to understand the code (you've even explicitly tagged on a "Tangent", knowing it's not relevant but finding it very interesting). Here's what I see as important to understanding the code:

    // Defining tr = Ixx + Iyy + Izz (the trace of the rotational inertia), one
    // can prove that the magnitude of the largest moment of inertia is upper
    // bound by |tr/2|. We want to scale epsilon by that maximum moment value;
    // finding it is too expensive, so we use |tr/2| as an acceptable proxy.

Nothing else you've noted (Imed, Imin, special case of Imin=0, lower bound on Imax, etc.) contributes to understanding the actual code here.

Is this really the place to teach about the linear algebra implications of the inertia matrix?

Code quote:

    // Denoting Imin and Imax as the smallest and largest possible moments of
    // inertia in a rotational inertia, denoting Imed as the intermediate moment
    // inertia, and denoting tr = Ixx + Iyy + Izz as the trace of the rotational
    // inertia, one can prove that if Imin == 0, then Imed == Imax == tr / 2 and
    // in all cases  0 ≤ Imin ≤ tr/3,  tr/3 ≤ Imed ≤ tr/2,  tr/3 ≤ Imax ≤ tr/2.
    // Hence: to check the validity of `this` rotational inertia, use an ε value
    // related to machine precision multiplied by Imax = tr/2.
    // Tangent: The product of inertia Iᵢⱼ with the largest absolute value is at
    // most half the largest moment of inertia, i.e., |Iᵢⱼ| ≤ Imax/2.

---

multibody/tree/rotational_inertia.h line 964 at r2 (raw file):

    // Tangent: The product of inertia Iᵢⱼ with the largest absolute value is at
    // most half the largest moment of inertia, i.e., |Iᵢⱼ| ≤ Imax/2.
    using std::abs;

nit: You have the function CalcMaximumPossibleMomentOfInertia(). Why have you duplicated that function and its documentation here?

---

multibody/tree/rotational_inertia.h line 966 at r2 (raw file):

    using std::abs;
    const double precision = 16 * std::numeric_limits<double>::epsilon();
    const T max_possible_inertia_moment = 0.5 * abs(Ixx + Iyy + Izz);

nit The only mention of |tr| comes in your tangent. But you've defined your upper bound by |tr| and not tr. (I accounted for that in my alternative text proposal above.) If you reject the proposal, you should, at least, consider tweaking your comments so that they better align with your code.

---

multibody/tree/rotational_inertia.h line 966 at r2 (raw file):

    using std::abs;
    const double precision = 16 * std::numeric_limits<double>::epsilon();
    const T max_possible_inertia_moment = 0.5 * abs(Ixx + Iyy + Izz);

nit: This function only gets called twice. Once in CouldBePhysicallyValid() and once in ThrowNotPhysicallyValid(). In both cases the values of Ixx, Iyy, and Izz are the actual principal moments of inertia. We shouldn't have to compute this proxy, we actually have the maximum moment.

It seems silly to preface each call with the cost of computing the principal moments and then to ignore that value by creating a bound that could be as much as 50% larger than we'd want.

It seems that the actual usage and the implementation could be better aligned.

Options:

1. Document them as being the principal moments (Imin, Imed, Imax) and throw out the proxy (leaving the call sites intact).
2. Double down on the generic Ixx, Iyy, and Izz and stop unnecessarily pre-computing the traces at the call sites.
3. Throw out the parameters entirely, and allow this function to decide whether it uses the principal moments or not. It shouldn't be subject to whether a call site has decided to use generic diagonal values or actual principal moments.

---

multibody/tree/rotational_inertia.h line 975 at r2 (raw file):

    const auto are_moments_near_positive =
        AreMomentsOfInertiaNearPositive(Ixx, Iyy, Izz, epsilon);
    const auto is_triangle_inequality_satisified = Ixx + Iyy + epsilon >= Izz &&

nit: typo

(Predating this PR, but as long as we're here...)

Suggestion:

satisfied 

---

multibody/tree/rotational_inertia.cc line 197 at r2 (raw file):

      const Vector3<double> p = CalcPrincipalMomentsOfInertia();
      if (!AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality(
              p(0), p(1), p(2))) {

nit: If you keep AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality() using the |tr/2| proxy value, then you should invoke it here with the actual "generic" Ixx, Iyy, and Izz and only calculate the actual principal moments if that test fails and you need them for the error message. But right now, you're paying the cost whether the inertia is obviously valid or not.

[mitiguy]

Good -- thank you for the scrutiny.

Reviewable status: 4 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits)

---

a discussion (no related file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

No tests?

There are no new tests (nor do I think new ones are needed).
This PR was to address the inconsistent epsilons.

There are tests that exercise the functions in question (indirectly -- as these are not public functions).
These tests are in the file rotational_inertia_test.cc e.g., see
GTEST_TEST(RotationalInertia, MakeFromMomentsAndProductsOfInertia).

---

multibody/tree/rotational_inertia.h line 592 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: Again, if you keep the |tr/2| proxy, you shouldn't be computing p, you should simply be passing the generic Ixx, Iyy, and Izz into AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality() and not waste the work of computing the principal moments for this test.

Although it is worthwhile to test Ixx, Iyy, Izz, the best litmus test is the principal moments. I added an example to clarify this.

Note: Subsequent PRs will make additional fixes. For example, as of before this PR (and is still), this function does NOT actually make all the checks alleged in the documentation above, i.e., as highlighted in the bold below:

/// - Ixx, Iyy, Izz and principal moments are all non-negative.
/// - Ixx, Iyy Izz and principal moments satisfy the triangle inequality

---

multibody/tree/rotational_inertia.h line 948 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

BTW Every positive number is non-zero, making "(non-zero)" a redundant elaboration.

Done.

---

multibody/tree/rotational_inertia.h line 950 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

BTW I'm not blocking on this, but I'll reiterate the following:

The justification of using epsilon (i.e., the original inertia satisfied the constraint, but the re-expressed values don't because of rounding error) is perfectly valid. This increases the space of "valid" inertia values to include those cases. Yay!

Remember, for every such inertia, there is another inertia that didn't satisfy the triangle inequality constraint but, due to re-expression's rounding error, it now falls within the epsilon, thus it gets incorrectly classified as valid.

It might be worth acknowledging that fact. That we'd rather include some ever-so-slightly invalid inertias rather than exclude the ones that were definitely valid but were penalized by finite precision during valid mathematical operations.

I just hate the idea of thinking of "epsilon" as some amazing panacea, especially when used like this where you're partitioning the value space into two regions instead of three.

Done. I added this note:

// Note: A side effect of ε is that some inertias are incorrectly classified
// as valid. We prefer to include some ever-so-slightly invalid inertias 
// rather than exclude ones that were definitely valid but were penalized by
// finite precision during valid mathematical operations. 

I am unsure why this needs to be reiterated.
I do not regard epsilon as an "amazing" panacea.
It is an artifact of finite precision mathematics.
I guess one might call it a "pragmatic solution" to a known finite precision problem.
If there is a relatively simple better way to separate the wheat and the chaff (in the context of this PR), please let me know.

---

multibody/tree/rotational_inertia.h line 955 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: You've defined Imax as the largest "possible" moment of inertia. And then you've asserted that tr/3 ≤ Imax ≤ tr/2.

This is an error. The true statement is that Imax = tr/2.

In other words, you've blurred whether Imax refers to the actual largest moment, or the possible largest moment. Your English says "possible", your math says "actual".

Now, if you'd defined Imax as the actual largest moment, then your inequality would be true.

See my follow up note on this documentation.

Done. Simplified significantly.

---

multibody/tree/rotational_inertia.h line 963 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

BTW You've flooded this comment with pedagogical trivia. As I read the comment and the code, you've made many assertions that don't directly reflect what is in the code or even necessary to understand the code (you've even explicitly tagged on a "Tangent", knowing it's not relevant but finding it very interesting). Here's what I see as important to understanding the code:

    // Defining tr = Ixx + Iyy + Izz (the trace of the rotational inertia), one
    // can prove that the magnitude of the largest moment of inertia is upper
    // bound by |tr/2|. We want to scale epsilon by that maximum moment value;
    // finding it is too expensive, so we use |tr/2| as an acceptable proxy.

Nothing else you've noted (Imed, Imin, special case of Imin=0, lower bound on Imax, etc.) contributes to understanding the actual code here.

Is this really the place to teach about the linear algebra implications of the inertia matrix?

Done. I agree. Good suggestions. Comment is significantly shortened

---

multibody/tree/rotational_inertia.h line 964 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: You have the function CalcMaximumPossibleMomentOfInertia(). Why have you duplicated that function and its documentation here?

I like your suggestion #3 below, so check out the changes.

To answer your question (no longer relevant after your suggested changes): This was a static function so there is no "this" rotational inertia. Hence no call to CalcMaximumPossibleMomentOfInertia().

---

multibody/tree/rotational_inertia.h line 966 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit The only mention of |tr| comes in your tangent. But you've defined your upper bound by |tr| and not tr. (I accounted for that in my alternative text proposal above.) If you reject the proposal, you should, at least, consider tweaking your comments so that they better align with your code.

Done. PTAL.

---

multibody/tree/rotational_inertia.h line 966 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: This function only gets called twice. Once in CouldBePhysicallyValid() and once in ThrowNotPhysicallyValid(). In both cases the values of Ixx, Iyy, and Izz are the actual principal moments of inertia. We shouldn't have to compute this proxy, we actually have the maximum moment.

It seems silly to preface each call with the cost of computing the principal moments and then to ignore that value by creating a bound that could be as much as 50% larger than we'd want.

It seems that the actual usage and the implementation could be better aligned.

Options:

1. Document them as being the principal moments (Imin, Imed, Imax) and throw out the proxy (leaving the call sites intact).
2. Double down on the generic Ixx, Iyy, and Izz and stop unnecessarily pre-computing the traces at the call sites.
3. Throw out the parameters entirely, and allow this function to decide whether it uses the principal moments or not. It shouldn't be subject to whether a call site has decided to use generic diagonal values or actual principal moments.

Done. I really like suggestion #3.
This change will be useful for subsequent PRs.

---

multibody/tree/rotational_inertia.h line 975 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: typo

(Predating this PR, but as long as we're here...)

Done. Good catch. Thanks Sean.

---

multibody/tree/rotational_inertia.cc line 197 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: If you keep AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality() using the |tr/2| proxy value, then you should invoke it here with the actual "generic" Ixx, Iyy, and Izz and only calculate the actual principal moments if that test fails and you need them for the error message. But right now, you're paying the cost whether the inertia is obviously valid or not.

Done. (I think).

[jwnimmer-tri]

Reviewable status: 5 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @mitiguy)

---

multibody/tree/rotational_inertia.h line 958 at r3 (raw file):

  // performed for `this` rotational inertia's principal moments of inertia or
  // for `this` rotational inertia's diagonal moments of inertia.
  boolean<T> AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality(

BTW it seems to me like this function does not meet our current standards for being inline: https://drake.mit.edu/styleguide/cppguide.html#Inline_Functions

[SeanCurtis-TRI]

A core source of confusion, flagged below.

Reviewed 2 of 2 files at r3, all commit messages.
Reviewable status: 6 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @mitiguy)

---

a discussion (no related file):

Previously, mitiguy (Mitiguy) wrote…

There are no new tests (nor do I think new ones are needed).
This PR was to address the inconsistent epsilons.

There are tests that exercise the functions in question (indirectly -- as these are not public functions).
These tests are in the file rotational_inertia_test.cc e.g., see
GTEST_TEST(RotationalInertia, MakeFromMomentsAndProductsOfInertia).

Technically, you did change the sensitivity of one call site (having gone from zero to non-zero epsilon). The fact that no test got angry about that means that the tests don't care about the epsilon value. Is the use of epsilon insufficiently important to justify a test?

---

multibody/tree/rotational_inertia.h line 592 at r2 (raw file):

Previously, mitiguy (Mitiguy) wrote…

Although it is worthwhile to test Ixx, Iyy, Izz, the best litmus test is the principal moments. I added an example to clarify this.

Note: Subsequent PRs will make additional fixes. For example, as of before this PR (and is still), this function does NOT actually make all the checks alleged in the documentation above, i.e., as highlighted in the bold below:

/// - Ixx, Iyy, Izz and principal moments are all non-negative.
/// - Ixx, Iyy Izz and principal moments satisfy the triangle inequality

The example you give here has persuaded me that the formulation you have of the triangle inequality test only works on the principal moments. Therefore, if that logic only works with principal moments, why give the caller the option to not use the principal moments? Doesn't that give the caller the power to make AreMomentsOfInertia....()`'s result invalid? Why do that?

Plain and simple AreMomentsOfInerita...() should do what it needs to do to always give a correct answer (for a given value of correct). Its correctness should not be dependent on the user turning on the "please be correct" switch.

---

multibody/tree/rotational_inertia.h line 948 at r2 (raw file):

Previously, mitiguy (Mitiguy) wrote…

Done.

Near zero is not generally correct and quite imprecise. It scales with the largest moment. However, I'm not going to block because of my pedantry. There's obviously something here that is very important for you to express. My inability to see it is not sufficient basis to further worry.

---

multibody/tree/rotational_inertia.h line 950 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

BTW "is typically much smaller than the largest possible principal moment of inertia"

I have two issues with this:

• The word "typically" hides a lot of sins. It doesn't say "always" or simply say that epsilon scales with the largest principal moment. This document implies that describing epsilon is important here, but then hiding the "atypical" case seems counter to that purpose. What is being hidden here is that otherwise valid inertias whose principal moments of inertia are simply "too small" also are deemed invalid. This is a typical Drake-ism knowing that we're targeting human-scale dynamical systems, but it might be worth being more clear about it.
• When I first read this phrase, my brain went in the wrong direction and I misinterpreted it. I believe a slight rephrasing would help people like me and not offend your sensibilities (although, if it does, you can skip the suggested text):

Because the re-expression error scales with the largest principal moment, the epsilon also scales (but relative to an upper-bound estimate of the largest principal moment). We also exclude rotational inertias whose overall scale are considered too small relative to Drake's expected robot workspace scale.

I guess the change in documentation simply made this concern go away; the word "typically" has left the building. I will still note that the question of "otherwise valid but too small for Drake" behavior still isn't mentioned anywhere.

---

multibody/tree/rotational_inertia.h line 975 at r2 (raw file):

Previously, mitiguy (Mitiguy) wrote…

Done. Good catch. Thanks Sean.

That was just my vs code yelling at me. :)

---

multibody/tree/rotational_inertia.h line 958 at r3 (raw file):

Previously, jwnimmer-tri (Jeremy Nimmer) wrote…

BTW it seems to me like this function does not meet our current standards for being inline: https://drake.mit.edu/styleguide/cppguide.html#Inline_Functions

Agreed.

---

multibody/tree/rotational_inertia.cc line 197 at r2 (raw file):

Previously, mitiguy (Mitiguy) wrote…

Done. (I think).

Done enough. (It doesn't have the savings I hoped for, but you've persuaded me that the triangle inequality test requires the principal moments, period; the savings were illusory from the beginning.)

---

multibody/tree/rotational_inertia.h line 959 at r3 (raw file):

  // for `this` rotational inertia's diagonal moments of inertia.
  boolean<T> AreMomentsOfInertiaNearPositiveAndSatisfyTriangleInequality(
      bool is_test_principal_moments_of_inertia) const {

BTW You like verbose parameters. If you keep this, perhaps you can be more concise?

use_principal_moments, perhaps?

It's hard to imagine that any reader would think you mean the principal moments of something else. And the imperative "use" sounds like the directive you mean this parameter to be, rather than the sometimes boolean predicate question that "is_test" is frequently used as.

To be clear, I assume this parameter will simply go away.

---

multibody/tree/rotational_inertia.h line 970 at r3 (raw file):

    using std::max;
    const double precision = 16 * std::numeric_limits<double>::epsilon();
    const T max_possible_inertia_moment = CalcMaximumPossibleMomentOfInertia();

nit: The maximum possible inertia moment has no value if you have the actual moments (including the maximum moment). We don't need to approximate the biggest moment if we have the biggest moment.

Implication:

1. If I've persuaded you that this should always compute the principal moment, then all code and documentation relating to "maximum possible" moment disappears completely from this function and precision is simply scaled by the largest principal moment.
2. If, it turns out, there is a reason to do the triangle inequality test on the diagonal of an arbitrary rotation inertia, then you should still only use the "maximum possible" moment if the user has requested you don't use the principal moments. E.g.,

    const Vector3<double> moments = is_test_principal_moments_of_inertia
                                        ? CalcPrincipalMomentsOfInertia()
                                        : ExtractDoubleOrThrow(get_moments());
    const double scale = is_test_principal_moments_of_inertia
                             ? moments.z()
                             : CalcMaximumPossibleMomentOfInertia();
    const double precision = 16 * std::numeric_limits<double>::epsilon();
    const T epsilon = precision * max(1.0, scale);

[mitiguy]

Ideally, confusion has been addressed by new comments and the reviewable comments about the sequence of proper triangle-inequality tests.

Reviewable status: 4 unresolved discussions, LGTM missing from assignee SeanCurtis-TRI(platform), commits need curation (https://drake.mit.edu/reviewable.html#curated-commits) (waiting on @mitiguy)

---

a discussion (no related file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

Technically, you did change the sensitivity of one call site (having gone from zero to non-zero epsilon). The fact that no test got angry about that means that the tests don't care about the epsilon value. Is the use of epsilon insufficiently important to justify a test?

New tests added which care about epsilon.

---

multibody/tree/rotational_inertia.h line 592 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

The example you give here has persuaded me that the formulation you have of the triangle inequality test only works on the principal moments. Therefore, if that logic only works with principal moments, why give the caller the option to not use the principal moments? Doesn't that give the caller the power to make AreMomentsOfInertia....()`'s result invalid? Why do that?

Plain and simple AreMomentsOfInerita...() should do what it needs to do to always give a correct answer (for a given value of correct). Its correctness should not be dependent on the user turning on the "please be correct" switch.

I disavow persuading you on this.
The triangle inequality test is valid for all inertia matrices, not just inertia matrices whose off-diagonal elements (products of inertia) are zero.

Necessary, but not sufficient: Although testing a rotational inertia's central moments of inertia is more conclusive than testing a generic rotational inertia, it is still an insufficient test. Failing a triangle inequality test is a definite test of invalidity. Passing a triangle inequality test does not validate the rotational inertia because a rotational inertia is always associated with an "about-point", and this RoationalInertia class does not have that information.

A proper sequence for the triangle-inequality test:
1 a. Test the user's rotational inertia (only need rotational inertia).
If it fails, the error is just due to the rotational inertia.

1 b. Test the user's principal moments of inertia (only need rotational inertia).
If it fails, the error is due to a combination of the rotational inertia and principal moments of inertia.

2 a. Shift the rotational inertia to center of mass (requires spatial inertia class).
If it fails, the error is due to a combination of the rotational inertia and center of mass position.

2 b. Test central principal moments of inertia (requires spatial inertia class).
If it fails, the error is due to a combination of the rotational inertia, center of mass position, and principal moments of inertia.

Note: A good error message only contains the pertinent information on why there is failure. It is confusing to complicate an error message with eigenvalues and principal moments of inertia when the user's matrix immediately fails the simplest triangle inequality test.
Improved error message construction will be addressed in a subsequent PR.

---

multibody/tree/rotational_inertia.h line 950 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

I guess the change in documentation simply made this concern go away; the word "typically" has left the building. I will still note that the question of "otherwise valid but too small for Drake" behavior still isn't mentioned anywhere.

No. The epsilon uses a max to be more permissive.
const T epsilon = precision * max(1.0, max_possible_inertia_moment);

I am perplexed as to why this causes you the concerns you stated.

What is being hidden here is that otherwise valid inertias whose principal moments of inertia are simply "too small" _also_ are deemed invalid.'' and We also exclude rotational inertias whose overall scale are considered too small relative to Drake's expected robot workspace scale.''

Note: Tests were added to ensure:

1. A zero rotational inertia is valid.
2. A rotational inertia with tiny negative moments of inertia does NOT throw an exception.
3. A rotational inertia with small negative moments of inertia DOES throw an exception.

---

multibody/tree/rotational_inertia.h line 958 at r3 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

Agreed.

Done.

---

multibody/tree/rotational_inertia.h line 959 at r3 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

BTW You like verbose parameters. If you keep this, perhaps you can be more concise?

use_principal_moments, perhaps?

It's hard to imagine that any reader would think you mean the principal moments of something else. And the imperative "use" sounds like the directive you mean this parameter to be, rather than the sometimes boolean predicate question that "is_test" is frequently used as.

To be clear, I assume this parameter will simply go away.

Done. Yes -- good suggestion. Thank you.
Changed parameter is_test_principal_moments_of_inertia to use_principal_moments
To be clear, this parameter should not go away.

---

multibody/tree/rotational_inertia.h line 970 at r3 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

nit: The maximum possible inertia moment has no value if you have the actual moments (including the maximum moment). We don't need to approximate the biggest moment if we have the biggest moment.

Implication:

1. If I've persuaded you that this should always compute the principal moment, then all code and documentation relating to "maximum possible" moment disappears completely from this function and precision is simply scaled by the largest principal moment.
2. If, it turns out, there is a reason to do the triangle inequality test on the diagonal of an arbitrary rotation inertia, then you should still only use the "maximum possible" moment if the user has requested you don't use the principal moments. E.g.,

    const Vector3<double> moments = is_test_principal_moments_of_inertia
                                        ? CalcPrincipalMomentsOfInertia()
                                        : ExtractDoubleOrThrow(get_moments());
    const double scale = is_test_principal_moments_of_inertia
                             ? moments.z()
                             : CalcMaximumPossibleMomentOfInertia();
    const double precision = 16 * std::numeric_limits<double>::epsilon();
    const T epsilon = precision * max(1.0, scale);

1. We do not always need to compute the principal moments of inertia.

2. I do not think the scale factor should change due to the argument passed.
   Example: the two nearly identical rotational matrices below are related. The 2nd matrix just contains the principal moments of inertia of the 1st. I think their scale factors should also be nearly identical.

[ 2          0        0.000010
  0          2        0
  0.00001    0        2 ]
scale = CalcMaximumPossibleMomentOfInertia() = 3;

[ 1.99999     0       0
  0           2       0 
  0           0       2.00001 ]
scale = moments.z() = 2.00001;

One might argue that the name of the function CalcMaximumPossibleMomentOfInertia() should be changed to something like CalcMaximumPossibleMomentOfInertiaIgnoringProductsOfInertia().
If so, I AGREE as it is more EXPLICIT.

---

multibody/tree/rotational_inertia.cc line 197 at r2 (raw file):

Previously, SeanCurtis-TRI (Sean Curtis) wrote…

Done enough. (It doesn't have the savings I hoped for, but you've persuaded me that the triangle inequality test requires the principal moments, period; the savings were illusory from the beginning.)

I was NOT trying to persuade you that for a generic RotationalInertia (which may contain non-zero products of inertia) that the triangle inequality test was invalid or useless.