Measurement in IAU

[mpusz]

Improve and productize a measurement class and use it in the definition of the IAU system based on https://web.archive.org/web/20131110215339/http://asa.usno.navy.mil/static/files/2014/Astronomical_Constants_2014.pdf.

[mpusz]

It seems that the IAU is not the only library customer for such a utility. #510 mentioned plenty of constants that are also specified with some uncertainty that might change over the years if we get a better way to measure them.

[mpusz]

We should also study what Python did for a similar feature: https://pythonhosted.org/uncertainties.

[mpusz]

https://en.wikipedia.org/wiki/Propagation_of_uncertainty

[RalphSteinhagen]

Some additional resources of interest might be:

• http://www.pennelynn.com/Documents/CUJ/HTML/14.03/MANNING/MANNING.HTM
• https://github.com/martukas/uncertain
• https://github.com/Gattocrucco/uncertainties-cpp

... the issue I see is that this can be very simple for measurement with uncorrelated stdev errors ... up to very complex when handling this thoroughly for correlated errors, n-dimensional systems, or errors that go beyond statistical r.m.s.-style errors.

[RalphSteinhagen]

^{_{xkcd.com CC BY-NC 2.5 DEED}}

Hi @mpusz and @mattkretz, I have been working a bit on the concept of propagation of uncertainties. The implementation is showcased in this compiler-explorer example. The core of this implementation is based on the UncertainValue struct, handling both uncorrelated and correlated uncertainties.

enum class Correlation { UnCorrelated, Correlated };

template <is_arithmetic_or_complex_v T, Correlation correlation>
struct UncertainValue;

template <is_arithmetic_or_complex_v T>
struct UncertainValue<T, Correlation::UnCorrelated> {
    using value_type = T;

    T value;                            /// stores the mean value
    T uncertainty = static_cast<T>(0);  /// stores the standard deviation of the value

   UncertainValue(T val, T uncert) : value(val), uncertainty(uncert) {}
};
template <typename T>
UncertainValue(T, T) -> UncertainValue<T, Correlation::UnCorrelated>;

template <is_arithmetic_or_complex_v T>
struct UncertainValue<T, Correlation::Correlated> {
    using value_type = T;

    T value;                            /// stores the mean value
    T uncertainty = static_cast<T>(0);  /// stores the standard deviation of the value
    std::uint64_t correlationID = 0UZ;  /// correlation ID - fully correlated if ID1==ID2 - issue: redundant 8bytes for each samples

   UncertainValue(T val, T uncert, std::uint64_t correlationID_ = get_unique_id()) : value(val), uncertainty(uncert), correlationID(correlationID_) {}
};
template <typename T>
UncertainValue(T, T, std::uint64_t) -> UncertainValue<T, Correlation::Correlated>;

Notably, I'd would like your input on:

• Operator Implementations: the operator+, '-, *, /implementations that largely follow the examples given [here](https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulae). They should work for fundamental types andstd::complex<>`.There seems to be a significant amount of boilerplate code, and I'm wondering if there's a more streamlined approach. How do other math libraries, particularly those dealing with SIMD, approach this? Is extending/wrapping the entire STL inevitable, or are there more efficient paths?
• Simplified Correlation Handling: I added a simplified propagation for correlated uncertainties ($\rho$ being either '1' or '0'). Is the treatment worth the effort from an 80/20 perspective? Complete and correct handling would always need to compute and propagate the full covariance matrix, which in all its rigour would be pretty heavy.
• SIMD Compatibility: (Specifically for @mattkretz) In our project, we'd be having a lot of containers carrying -- for example -- std::vector<UncertainValue<float>>. Are these structures suitable to simdize them into collections of values and uncertainties? We want to be able to re-use as much of the simd-ified STL library as possible.
• Performance Optimization: are there specific areas in the implementation where the performance cout be significantly improved without compromising accuracy?
• Error Propagation in Complex Functions: How can error propagation be more accurately/efficiently handled in complex arithmetic operations, especially when dealing with non-linear functions?
• Non-linear Functions in general: The Taylor-expansion of some -- for example -- like ln(x) isn't suitable for all ranges but probably would require $\sigma_f = 0.5 * ( |ln(x+\sigma) - ln(x)| + |ln(x-\sigma) - ln(x)|)$... do you know of better alternatives ... would be great to have compile-time differential functions in C++. 😁

To give a flavour of the operators needed for the propagation and to give the flavour of the boilerplate code I mentioned above:

template<typename T, typename U>
requires UncertainValueType<T> || UncertainValueType<U>
constexpr auto operator+(const T& lhs, const U& rhs) noexcept {
  using ValueTypeT = UncertainValueType_t<T>;
  using ValueTypeU = UncertainValueType_t<U>;
  if constexpr (UncertainValueType<T> && UncertainValueType<U>) {
      if constexpr (is_complex_v<ValueTypeT> || is_complex_v<ValueTypeU>) {
          // we are dealing with complex numbers        
          if constexpr (CorrelatedUncertainValue<T> && CorrelatedUncertainValue<U>) {
              if (areCorrelated(lhs, rhs)) {
                  // add value and uncertainties as if they are vectors in 2D space because they are correlated.
                  return T{lhs.value + rhs.value, lhs.uncertainty + rhs.uncertainty, lhs.correlationID};
              }
          }
          // values are not both complex and/or not both correlated -> use the standard uncorrelated calculation.          
          ValueTypeT newUncertainty = {
              std::hypot(std::real(lhs.uncertainty), std::real(rhs.uncertainty)),
              std::hypot(std::imag(lhs.uncertainty), std::imag(rhs.uncertainty))
          };
          return T{lhs.value + rhs.value, newUncertainty};
      } else { 
          // both ValueType[T,U] are arithmetic uncertainties
          if constexpr (CorrelatedUncertainValue<T> && CorrelatedUncertainValue<U>) { 
              if (areCorrelated(lhs, rhs)) {
                  return T{lhs.value + rhs.value, lhs.uncertainty + rhs.uncertainty, lhs.correlationID};
              }
          }
          return T{lhs.value + rhs.value, std::hypot(lhs.uncertainty, rhs.uncertainty)};
      }
  } else if constexpr (UncertainValueType<T> && is_arithmetic_or_complex_v<ValueTypeU>) {
      return T{lhs.value + rhs, lhs.uncertainty};
  } else if constexpr (is_arithmetic_or_complex_v<ValueTypeT> && UncertainValueType<U>) {
      return U{lhs + rhs.value, rhs.uncertainty};
  } else {
      static_assert(std::is_arithmetic_v<ValueTypeT> && std::is_arithmetic_v<ValueTypeU>);
      return lhs + rhs; // unlikely to be called due to default '+' definition
  }
}

I appreciate any insights, comments, or suggestions you can provide to improve this implementation. Your expertise in this area is invaluable.

[mpusz]

Hi @RalphSteinhagen, thanks for sharing this nice code! I am on vacation now, so I did not have much time for this review. I will probably come back to you with more details when I will learn more about handling uncertainties.

For now, a really short feedback:

• is_xxx_v is a good name for a type trait but not for a concept
• is_arithmetic_or_complex_v might be too restricting for some use cases
• Your types definitely need constexpr and [[nodiscard]]
• The members should not be publically exposed (the user should not be able to directly mutate only one of the members)
• UncertainValueType concept could do is_specialization_of check instead

[RalphSteinhagen]

• Your types definitely need constexpr and [[nodiscard]]

Absolutely ... 🙈 ... was a bit q-n-d for this concept.

• UncertainValueType concept could do is_specialization_of check instead

I'd prefer to test w.r.t. the functionality rather than the actual implementation. This should help composition and help users to have parallel functionally equivalent implementation of the struct.

• The members should not be publically exposed ...

This is IMO a longer discussion to be had. There are pros/cons for both sides that need to be evaluated.

[mpusz]

ISO 80000-6 provides the following text output for uncertainty:

l = 2,347 82(32) m

Which means 2,34782 m +−0,000 32 m.

[JohelEGP]

That's very similar to the example in Part 1's §7.2.4.

[mpusz]

Resolved by #819.

The measurement class from example/include/measurement.h is now productized as mp_units::utility::uncertain<T> in the mp-units::utility component. The rename follows VIM terminology: a measurement is the process, the value/uncertainty pair with a unit is a
measurement result, and the bare pair is neither — it is just an uncertain number.

The part that took the most design work was where the uncertainty should live. Constants remain units, so they keep their exact symbolic magnitudes, full cancellation, and do not force a representation type on the user. A measured constant additionally carries its
relative standard uncertainty as exact magnitude metadata:

inline constexpr struct newtonian_constant_of_gravitation final :
  named_constant<"G", mag_ratio<667'430, 100'000> * mag_power<10, -11> * cubic(si::metre) /
                      si::kilogram / square(si::second),
                 relative_standard_uncertainty{mag_ratio<22, 10> * mag_power<10, -5>}> {}
  newtonian_constant_of_gravitation;

The uncertainty is relational: expressed in its own unit a measured constant is exactly 1, and it is the conversion factor between units that carries the uncertainty. So 2 * M_SUN stays an exact statement, and a solar mass expressed in Earth masses is exact
too, because the shared G cancels symbolically. The relative uncertainty is materialized only when a quantity leaves units in which the constant cancels:

quantity mass = measurement_of(iau::solar_mass);  // 1 ± 0 M_☉, exact
std::println("{}", mass.in(si::kilogram));        // 1.98892e+30 ± 4.37562e+25 kg

Representation types opt into this through the new fold_conversion_uncertainty customization point, so a custom uncertainty-aware type gets the same behaviour and every other representation is unaffected (the derivation is not even instantiated for them).

On the IAU system specifically: after annotating G from Astronomical Constants 2014 (uᵣ = 2.2 × 10⁻⁵), nothing else in IAU needed annotating. The IAU 2015 Resolution B3 values are nominal — exact by definition precisely so that they can serve as conversion units — so the type system now states which constants are exact and which are measured, a distinction that previously lived only in comments. HEP was audited separately and 30 CODATA 2014/2018/2022 constants were annotated with their published relative uncertainties, including μ₀, which became a measured constant with the 2019 SI redefinition.

Docs:

• Working with measurement uncertainty
• Weighing the Earth example