Unit Testing Numerical Routines

Published

2024-11-12

Reading Time

7 minutes.

How would you unit test this function?

std::tuple<double, double, double> ecef_to_lla(const Eigen::Vector3d& position_ecef_m);

This prototype takes a position expressed as cartesian coordinates in the Earth Centered, Earth Fixed (ECEF) reference frame and converts it into latitude, longitude, and altitude. If you are not familiar with aerospace coordinate systems you will still be able to follow along. Think of it as a fancy cartesian to spherical coordinate transform.

As a prudent programmer, you probably do not consider your code finished until it has passed a rigorous set of unit tests. However, numerical routines such as this one present difficulties for unit testing. How can you be confident the function is correct in general? What if you don't know the expected output? What if the underlying algorithm changes?

Level 1: Magic Numbers

Your first attempt at testing might be to pull out an existing script and copy-paste the inputs and expected results into your unit test. For our routine, you might use this website. You might create something like this:

TEST_CASE("Check ecef_to_lla against web-app result", "[magic numbers]") {
  // Input point, reversed from output using
  // https://www.oc.nps.edu/oc2902w/coord/llhxyz.htm
  const Eigen::Vector3d pos{1318610., -4504041., 4306196.};

  // Do calculation
  const auto [lat, lon, alt] = ecef_to_lla(pos);

  // Check results to 0.001 degrees and 1m
  REQUIRE_THAT(lat, WithinAbs(42.730, 0.001));
  REQUIRE_THAT(lon, WithinAbs(-73.682, 0.001));
  REQUIRE_THAT(alt, WithinAbs(1000.0, 1.));
}

This is certainly better than not having a test. If the script used to generate the test case is trusted and stable, it's a great way to verify you meet requirements. However, if the algorithm is unstable, it may be difficult to maintain. Which version of c_ECEF2geo_final_v2.m did you use to make this unit test? It also doesn't help if your original implementation was wrong.¹

Level 2: Trivial Points

Instead of focusing on input-output pairs, let's make sure the function is returning results that make sense. This requires you to build an intuitive understanding of the system under test.

For example, we know that by definition the x-axis of the ECEF coordinate system must pass through the point at the intersection of the equatorial plane and the prime meridian. We also know that because this point lies on the equator, it must have a distance from the earth's center equal to the ellipsoid's semi-major axis. We can construct such a point, and assert it lands where expected.

TEST_CASE("Point at the intersection of the ellipsoid and the x-axis",
          "[hand calculations]") {
  // Point on the ellipsoid at y=0 and z=0
  const Eigen::Vector3d pos{semi_major_axis_m, 0, 0};

  // Do calculation
  const auto [lat, lon, alt] = ecef_to_lla(pos);

  // Check results are at the intersection of prime meridian and equator.
  REQUIRE_THAT(lat, WithinAbs(0.0, 0.001));
  REQUIRE_THAT(lon, WithinAbs(0.0, 0.001));
  REQUIRE_THAT(alt, WithinAbs(0.0, 1.0));
}

Other points worth testing might be the North Pole and points directly above those other points. The key here is to look for inputs where the output can be trivially derived.² Zeros are a good place to start.

Trivial points are great for developing an understanding, but it is still possible for errors to be lurking in the points you do not test.

Level 3: Property-Based Testing

What if we could automatically check trivial points over a wide variety of inputs? Unit testing frameworks like Catch2 allow this with property-based testing. The idea is to generate random input, and then assert that some property of the output holds, rather than checking for specific output.

For example, the sign of a point's latitude is controlled entirely by the sign of its Z coordinate. Using the GENERATE macro, we can pick random points to test this property on.

TEST_CASE("Latitude has correct sign", "[property based]") {
  constexpr int N = 10;
  const double z =
      GENERATE(take(N, random(-semi_major_axis_m, semi_major_axis_m)));

  // This code is reached 10 times, with z taking on a random value each time
}

const double z = GENERATE(...) tells Catch to create new unit tests with z taking on different values based on the rules inside the macro.

These rules are called generators, and they are composable. Here take(N, ...) is a generator that returns the first N values from another generator. We use this because we do not want an infinite number of examples.

Next, random(min, max) generates random numbers between min and max. We choose, somewhat arbitrarily, to generate points that are between the top and bottom of the earth (approximately). The full test follows:

TEST_CASE("Latitude has correct sign", "[property based]") {
  const double z =
      GENERATE(take(10, random(-semi_major_axis_m, semi_major_axis_m)));

  // this code is reached 10 times, with z taking on a random value each time

  const double x = 1318610.;
  const double y = -4504041.;

  // Do computation
  const auto [lat, lon, alt] = ecef_to_lla({x, y, z});

  // Check the property that `lat` and `z` have the same sign
  REQUIRE(std::signbit(lat) == std::signbit(z));
}

Finding these properties requires an understanding of your function, but developing that understanding is essential for writing correct code. For conversion functions like this, it can often be helpful to assert that composition with the inverse (ecef_to_lla(lla_to_ecef(...))) yields the identity.

It is also important for these properties to hold true for all implementations that meet your requirements. They should be based in math or physics, rather than your implementation.

Level 4: Checking Against a Specification

If you are lucky, you can test against a specification. In this case, the WGS84 ellipsoid and the ECEF frame derived from it are defined by the National Geospatial-Intelligence Agency in NGA.STND.0036_1.0.0_WGS84. Using a specification is very similar to the magic number technique, but instead of coming from a poorly change-controlled matlab script, it comes from a peer-reviewed standard.

The WGS84 specification contains a list of reference stations, with their positions given in ECEF and LLA coordinates. We will test the implementation against these values.

The most basic usage of GENERATE is to give it a comma-separated list of values you want your variable to take on. Here, we give it a list of tuples containing the name, ECEF coordinates, latitude, longitude, and ellipsoid height of each station.

TEST_CASE("Test against WGS84 Reference Stations", "[spec]") {
  using Catch::Matchers::WithinAbs;

  // Air Force Reference Station Locations from NGA.STND.0036_1.0.0_WGS84
  const auto [name, x, y, z, lat_spec, lon_spec, alt_spec] = GENERATE(
      std::make_tuple("Colorado Springs", -1248599.695, -4819441.002,
                      3976490.117, 38.80293817, 255.47540411, 1911.778),
      std::make_tuple("Ascension", 6118523.866, -1572350.772, -876463.909,
                      -7.95132931, 345.58786964, 106.281),
      std::make_tuple("Diego Garcia", 1916196.855, 6029998.797, -801737.183,
                      -7.26984216, 72.37092367, -64.371),
      std::make_tuple("Kwajalein", -6160884.028, 1339852.169, 960843.154,
                      8.72250188, 167.73052378, 39.652),
      std::make_tuple("Hawaii", -5511980.264, -2200246.752, 2329481.004,
                      21.56149239, 201.76066695, 425.789),
      std::make_tuple("Cape Canaveral", 918988.062, -5534552.894, 3023721.362,
                      28.48373823, 279.42769502, -24.083));

  // NGA data uses longitude east, even for the western hemisphere. Wrap that to
  // have negative degrees
  const auto lon_spec_wrapped = lon_spec > 180 ? lon_spec - 360 : lon_spec;

  // Compute the latitude, longitude and ellipsoid height
  const auto [lat_computed, lon_computed, alt_computed] =
      ecef_to_lla({x, y, z});

  // Check that the results match the spec values
  REQUIRE_THAT(lat_computed, WithinAbs(lat_spec, 1e-6));
  REQUIRE_THAT(lon_computed, WithinAbs(lon_spec_wrapped, 1e-6));
  REQUIRE_THAT(alt_computed, WithinAbs(alt_spec, 0.1));
}

You are of course lucky if your outputs are so well defined. For most routines, you will have to fall back on the previous techniques. With a good set of sanity checks, you will save time and money when it comes to performing Verification & Validation. In fact, the more of your Verification & Validation plan that can be written up as unit tests, and integrated into your CI/CD pipeline, the faster your subject matter experts and software engineers will be able to iterate.