Skip to content
Permalink

Comparing changes

Choose two branches to see what’s changed or to start a new pull request. If you need to, you can also or learn more about diff comparisons.

Open a pull request

Create a new pull request by comparing changes across two branches. If you need to, you can also . Learn more about diff comparisons here.
base repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected base ref is valid, then try again.
Loading
base: 5d762b2fea2a
Choose a base ref
...
head repository: mockingbirdnest/Principia
Failed to load repositories. Confirm that selected head ref is valid, then try again.
Loading
compare: f17f9fee970a
Choose a head ref
  • 10 commits
  • 5 files changed
  • 1 contributor

Commits on Jun 24, 2020

  1. A failure in Diagonalize.

    pleroy committed Jun 24, 2020
    Copy the full SHA
    fb6486d View commit details

Commits on Jul 3, 2020

  1. Copy the full SHA
    f81cea0 View commit details

Commits on Jul 4, 2020

  1. Copy the full SHA
    e0b1381 View commit details
  2. Copy the full SHA
    9ae88e5 View commit details
  3. A more thorough test.

    pleroy committed Jul 4, 2020
    Copy the full SHA
    44d2e2b View commit details
  4. Copy the full SHA
    b103097 View commit details
  5. Cleanup of Diagonalize.

    pleroy committed Jul 4, 2020
    Copy the full SHA
    a68cb16 View commit details
  6. Readying for review.

    pleroy committed Jul 4, 2020
    Copy the full SHA
    7acfb03 View commit details
  7. Typo.

    pleroy committed Jul 4, 2020
    Copy the full SHA
    8a2f021 View commit details
  8. Merge pull request #2626 from pleroy/2611

    Properly deal with two very close eigenvalues in Diagonalize
    pleroy authored Jul 4, 2020
    Copy the full SHA
    f17f9fe View commit details
Showing with 127 additions and 17 deletions.
  1. +3 −3 geometry/symmetric_bilinear_form.hpp
  2. +34 −12 geometry/symmetric_bilinear_form_body.hpp
  3. +87 −1 geometry/symmetric_bilinear_form_test.cpp
  4. +1 −1 journal/player_test.cpp
  5. +2 −0 journal/profiles.cpp
6 changes: 3 additions & 3 deletions geometry/symmetric_bilinear_form.hpp
Original file line number Diff line number Diff line change
@@ -88,10 +88,10 @@ class SymmetricBilinearForm {

private:
// Given a matrix that contains in columns eigenvectors for a form, picks the
// column with the largest norm and return its normalized value. This is
// useful to extract eigenvectors when eigenvalues are known.
// column with the largest norm and return its value. This is useful to
// extract eigenvectors when eigenvalues are known.
template<typename S>
static R3Element<double> PickEigenvector(R3x3Matrix<S> const& matrix);
static R3Element<S> PickEigenvector(R3x3Matrix<S> const& matrix);

// All the operations on this class must ensure that this matrix remains
// symmetric.
46 changes: 34 additions & 12 deletions geometry/symmetric_bilinear_form_body.hpp
Original file line number Diff line number Diff line change
@@ -21,6 +21,7 @@ using quantities::Angle;
using quantities::ArcCos;
using quantities::Cos;
using quantities::Sqrt;
using quantities::Square;
using quantities::si::Radian;

template<typename Scalar,
@@ -149,6 +150,7 @@ typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
return {form, Rotation<Eigenframe, Frame>::Identity()};
}

// When det_B is close to -2 or 2, two of the eigenvalues are close.
R3x3Matrix<double> const B = A_minus_qI / p;
double const det_B = B.Determinant();
Angle const θ = ArcCos(std::clamp(det_B, -2.0, 2.0) * 0.5);
@@ -170,21 +172,38 @@ typename SymmetricBilinearForm<Scalar, Frame, Multivector>::
{zero, zero, α₂}));

// Use the Cayley-Hamilton theorem to efficiently find the eigenvectors. The
// m matrices contain, in columns, eigenvectors for the corresponding α.
// However it's possible for a column to be identically 0. To deal with this
// the call to PickEigenvector extracts the column with the largest norm, and
// we make sure that the eigenvectors remain orthonormal.
// mᵢ matrices contain, in columns, eigenvectors for the corresponding αᵢ.
R3x3Matrix<Scalar> const A_minus_α₀I = A - α₀ * I;
R3x3Matrix<Scalar> const A_minus_α₁I = A - α₁ * I;
R3x3Matrix<Scalar> const A_minus_α₂I = A - α₂ * I;
auto const m₀ = A_minus_α₁I * A_minus_α₂I;

// If two eigenvalues are very close, we want to be as precise as possible for
// the eigenvector that's not close to the other two. That happens for the
// inertia of an object that is a disc or a needle. So we locate the third
// eigenvalue and read its eigenvector directly from the matrix mᵢ. The other
// eigenvectors are orthogonal, but they may be far from the truth because of
// the singularity.
auto const m₁ = A_minus_α₂I * A_minus_α₀I;
auto const v₀ = Vector<double, Frame>(PickEigenvector(m₀));
auto const v₁ = Vector<double, Frame>(PickEigenvector(m₁)).
OrthogonalizationAgainst(v₀);
std::unique_ptr<Rotation<Eigenframe, Frame> const> rotation;
if (α₁ - α₀ < α₂ - α₁) {
auto const m₂ = A_minus_α₀I * A_minus_α₁I;
auto const v₂ =
Normalize(Vector<Square<Scalar>, Frame>(PickEigenvector(m₂)));
auto const v₁ = Normalize(Vector<Square<Scalar>, Frame>(PickEigenvector(m₁))
.OrthogonalizationAgainst(v₂));
rotation =
std::make_unique<Rotation<Eigenframe, Frame>>(Wedge(v₁, v₂), v₁, v₂);
} else {
auto const m₀ = A_minus_α₁I * A_minus_α₂I;
auto const v₀ =
Normalize(Vector<Square<Scalar>, Frame>(PickEigenvector(m₀)));
auto const v₁ = Normalize(Vector<Square<Scalar>, Frame>(PickEigenvector(m₁))
.OrthogonalizationAgainst(v₀));
rotation =
std::make_unique<Rotation<Eigenframe, Frame>>(v₀, v₁, Wedge(v₀, v₁));
}

Rotation<Eigenframe, Frame> const rotation{v₀, v₁, Wedge(v₀, v₁)};
return {form, rotation};
return {form, *rotation};
}

template<typename Scalar,
@@ -212,19 +231,22 @@ template<typename Scalar,
typename Frame,
template<typename, typename> typename Multivector>
template<typename S>
R3Element<double>
R3Element<S>
SymmetricBilinearForm<Scalar, Frame, Multivector>::PickEigenvector(
R3x3Matrix<S> const& matrix) {
static R3Element<double> const e₀{1, 0, 0};
static R3Element<double> const e₁{0, 1, 0};
static R3Element<double> const e₂{0, 0, 1};
std::array<R3Element<S>, 3> vs = {matrix * e₀, matrix * e₁, matrix * e₂};

// It's possible for a column to be identically 0. To deal with this we
// extract the column with the largest norm.
std::sort(vs.begin(),
vs.end(),
[](R3Element<S> const& left, R3Element<S> const& right) {
return left.Norm²() < right.Norm²();
});
return NormalizeOrZero(vs.back());
return vs.back();
}

template<typename Frame, template<typename, typename> typename Multivector>
88 changes: 87 additions & 1 deletion geometry/symmetric_bilinear_form_test.cpp
Original file line number Diff line number Diff line change
@@ -14,6 +14,8 @@
#include "serialization/geometry.pb.h"
#include "testing_utilities/almost_equals.hpp"
#include "testing_utilities/componentwise.hpp"
#include "testing_utilities/is_near.hpp"
#include "testing_utilities/numerics_matchers.hpp"
#include "testing_utilities/matchers.hpp"
#include "testing_utilities/vanishes_before.hpp"

@@ -25,10 +27,13 @@ using quantities::Length;
using quantities::Pow;
using quantities::Square;
using quantities::si::Metre;
using testing_utilities::AbsoluteErrorFrom;
using testing_utilities::AlmostEquals;
using testing_utilities::Componentwise;
using testing_utilities::EqualsProto;
using testing_utilities::IsNear;
using testing_utilities::VanishesBefore;
using testing_utilities::operator""_⑴;
using ::testing::Eq;

class SymmetricBilinearFormTest : public ::testing::Test {
@@ -355,7 +360,7 @@ TEST_F(SymmetricBilinearFormTest, Diagonalize) {
AlmostEquals(1, 0),
VanishesBefore(1, 0)));
EXPECT_THAT(f_eigensystem.rotation.Inverse()(w₂),
Componentwise(VanishesBefore(1, 1),
Componentwise(VanishesBefore(1, 2),
VanishesBefore(1, 2),
AlmostEquals(1, 0)));
}
@@ -422,8 +427,89 @@ TEST_F(SymmetricBilinearFormTest, Diagonalize) {
+9.56737022360950050e+02}}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

EXPECT_THAT(f_eigensystem.rotation.quaternion().Norm(),
AlmostEquals(1.0, 1)) << f;
}

// A case with two eigenvalues that are very close (exact relative error
// 2e-15) that used to yield a non-unit quaternion because of a missing
// normalization. Found in game in #2611.
{
auto const f = MakeSymmetricBilinearForm<World>(
R3x3Matrix<double>({{+6.25360854308065672e+00,
+2.24243333089292812e-01,
+1.68316543009972008e-02},
{+2.24243333089292812e-01,
+1.09414207983843497e+01,
+3.52669451554594282e-01},
{+1.68316543009972008e-02,
+3.52669451554594282e-01,
+6.26937749903824937e+00}}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

EXPECT_THAT(f_eigensystem.rotation.quaternion().Norm(),
AlmostEquals(1.0, 0)) << f;
Vector<double, Eigenworld> const e₀({1, 0, 0});
Vector<double, Eigenworld> const e₁({0, 1, 0});
Vector<double, Eigenworld> const e₂({0, 0, 1});

// Real eigenvectors obtained with Mathematica. Note how bad the first two
// eigenvectors are: that's because the object is very nearly a disc, and
// our computation yields eigenvectors that are roughly 45° from the real
// ones. But at least they are in the right plane and determine a rotation
// that correctly aligns the isolated eigenvector.
EXPECT_THAT(f_eigensystem.rotation(e₀),
Componentwise(
AbsoluteErrorFrom(-0.71267592684216601514, IsNear(0.7_⑴)),
AbsoluteErrorFrom(+0.086267747252993862323, IsNear(0.01_⑴)),
AbsoluteErrorFrom(-0.69616872888945048034, IsNear(0.3_⑴))));
EXPECT_THAT(f_eigensystem.rotation(e₁),
Componentwise(
AbsoluteErrorFrom(-0.69988076924532898781, IsNear(0.3_⑴)),
AbsoluteErrorFrom(-0.02018794239465783510, IsNear(0.07_⑴)),
AbsoluteErrorFrom(+0.71397433835008141314, IsNear(0.7_⑴))));
EXPECT_THAT(f_eigensystem.rotation(e₂),
Componentwise(AlmostEquals(+0.04753874357012595212, 10),
AlmostEquals(+0.99606742882485799451, 1),
AlmostEquals(+0.07476459786563599011, 4)));
}

// A case similar to the previous one, but constructed so that the two largest
// eigenvalues are very close (a needle).
{
auto const f = MakeSymmetricBilinearForm<World>(
R3x3Matrix<double>({{+3.02958892130780040082,
+0.34454179629510833170,
-3.74544524094010311734},
{+0.34454179629510833170,
+9.98296957696546981827,
+0.18513433665111470179},
{-3.74544524094010311734,
+0.18513433665111470179,
+7.98744150172682978089}}));
auto const f_eigensystem = f.Diagonalize<Eigenworld>();

EXPECT_THAT(f_eigensystem.rotation.quaternion().Norm(),
AlmostEquals(1.0, 0)) << f;
Vector<double, Eigenworld> const e₀({1, 0, 0});
Vector<double, Eigenworld> const e₁({0, 1, 0});
Vector<double, Eigenworld> const e₂({0, 0, 1});

// Same comment as above regarding the last two eigenvectors.
EXPECT_THAT(f_eigensystem.rotation(e₀),
Componentwise(AlmostEquals(+0.88005120297326585654, 1),
AlmostEquals(-0.04350022098854655144, 1),
AlmostEquals(+0.47288223789782508101, 2)));
EXPECT_THAT(f_eigensystem.rotation(e₁),
Componentwise(
AbsoluteErrorFrom(-0.42509236497268917818, IsNear(0.07_⑴)),
AbsoluteErrorFrom(+0.37170808088366222608, IsNear(1.1_⑴)),
AbsoluteErrorFrom(+0.82530575173550731715, IsNear(0.2_⑴))));
EXPECT_THAT(f_eigensystem.rotation(e₂),
Componentwise(
AbsoluteErrorFrom(+0.21167513171658507292, IsNear(0.1_⑴)),
AbsoluteErrorFrom(+0.92732994849715299942, IsNear(1.6_⑴)),
AbsoluteErrorFrom(-0.30863053191969508327, IsNear(0.3_⑴))));
}
}

2 changes: 1 addition & 1 deletion journal/player_test.cpp
Original file line number Diff line number Diff line change
@@ -89,7 +89,7 @@ TEST_F(PlayerTest, DISABLED_SECULAR_Debug) {
// An example of how journaling may be used for debugging. You must set
// |path| and fill the |method_in| and |method_out_return| protocol buffers.
std::string path =
R"(P:\Public Mockingbird\Principia\Crashes\2530\JOURNAL.20200416-093807)"; // NOLINT
R"(P:\Public Mockingbird\Principia\Crashes\2611\JOURNAL.20200623-194313)"; // NOLINT
Player player(path);
int count = 0;
while (player.Play(count)) {
2 changes: 2 additions & 0 deletions journal/profiles.cpp
Original file line number Diff line number Diff line change
@@ -76,6 +76,8 @@ std::uint64_t SerializePointer(T* t) {

} // namespace

// To remove the check, define this macro to be:
// auto aa = (a); auto bb = (b);
#define PRINCIPIA_CHECK_EQ(a, b) CHECK((a) == (b))
#define PRINCIPIA_SET_VERBOSE_LOGGING 1