odem_rs_util/

random_variable.rs

//! This module contains support for random variables which are useful for
//! aggregating statistical data of measurements taken during a simulation run.
//!
//! This module provides a generic implementation of a random variable collector
//! for statistical analysis of numerical data, measured during a simulation
//! run. It supports the estimation of up to four statistical moments
//! ([arithmetic mean], [variance], [skew], and [kurtosis]) and allows users to
//! select the number of moments to estimate at compile time.
//!
//! The implementation is inspired by SLX's `random_variable` and generalizes
//! the concept to handle numerical data points of arbitrary types. It employs
//! numerically stable online algorithms based on sums of differences to the
//! estimated mean, as described in the work by [Pébaÿ et al.].
//!
//! The user can choose between unweighted and weighted random variables, with
//! support for tracking the development of values over time using a temporal
//! random variable. The module provides methods to estimate various statistical
//! properties, including mean, variance, skew, and kurtosis, along with
//! population variants of these metrics.
//!
//! ## Statistical Moments
//!
//! The four statistical moments supported by this module are:
//! - **Mean**: Location of the distribution, indicating its center.
//! - **Variance**: Size of the distribution, providing a location-independent
//!   distance from the mean.
//! - **Skewness**: Asymmetry of the distribution, offering a location and
//!   variance-independent measure for the balance of values left and right of
//!   the arithmetic mean.
//! - **Kurtosis**: Measure for the number and severity of outliers by weighing
//!   samples higher the further from the mean they get.
//!
//! ## Usage
//!
//! To use this module, create instances of `RandomVariable` with appropriate
//! types and choose the desired number of moments to estimate at compile time
//! (between 0 and 4).
//!
//! ```
//! // create an unweighted random variable estimating all four moments
//! # use odem_rs_util::random_variable::RandomVariable;
//! let rv = RandomVariable::new::<4>();
//!
//! // tabulate data points
//! rv.tabulate(42);
//! rv.tabulate(12);
//! rv.tabulate(15);
//! // ...
//!
//! // retrieve statistical properties
//! let mean = rv.mean();
//! let variance = rv.variance();
//! let skew = rv.skew();
//! let kurtosis = rv.kurtosis();
//! ```
//!
//! For detailed information on each method and trait, refer to the documentation of individual
//! items within the module.
//!
//! [Pébaÿ et al.]: https://link.springer.com/article/10.1007/s00180-015-0637-z
//! [arithmetic mean]: https://en.wikipedia.org/wiki/Arithmetic_mean
//! [variance]: https://en.wikipedia.org/wiki/Variance
//! [skew]: https://en.wikipedia.org/wiki/Skewness
//! [kurtosis]: https://en.wikipedia.org/wiki/Kurtosis

use crate::error::{Causality, NegativeWeight};
use core::{
	cell::Cell,
	cmp::Ordering,
	convert::Infallible,
	fmt,
	num::{NonZero, Wrapping},
	ops::{Add, Sub},
};

/* ********************************************************** Random Variable */

/// A simple collector for statistical data, inspired by SLX's `random_variable`
/// and generalized to numerical data points of an arbitrary type.
///
/// This implementation uses numerically stable online algorithms using sums of
/// differences to the estimated mean described in [Pébaÿ et al.] rather than
/// SLX's naïve approach of using sums directly. It supports the sampling of
/// weighted values and the subsequent estimation of up to four statistical
/// moments ([arithmetic mean], [variance], [skew], and [kurtosis]).
/// The user can select the number of statistical moments they want to estimate
/// by supplying a compile-time constant between `0` (only estimating the
/// minimum and maximum) and `4` (additionally estimating mean, variance, skew,
/// and kurtosis).
///
/// All four moments are useful to describe the kind of distribution:
/// - *mean* tells us the location of the distribution by giving us its center,
/// - *variance* measures the size of the distribution by giving us the
///   location-independent distance from the mean,
/// - *skew* describes the asymmetry of the distribution by giving us a location
///   and variance-independent measure for the balance of values left and right
///   of the arithmetic mean, and
/// - *kurtosis* is a measure for the number and severity of outliers by
///   weighing samples higher the further from the mean they get.
///
/// To sum it up: the former two statistical moments talk about the location and
/// size, and the latter two statistical moments talk about the shape of the
/// distribution. Since we usually cannot sample the whole population of values,
/// efforts have been made to numerically correct the [sampling bias] inherent
/// in sampling only a subset of a population.
///
/// [Pébaÿ et al.]: https://link.springer.com/article/10.1007/s00180-015-0637-z
/// [arithmetic mean]: https://en.wikipedia.org/wiki/Arithmetic_mean
/// [variance]: https://en.wikipedia.org/wiki/Variance
/// [skew]: https://en.wikipedia.org/wiki/Skewness
/// [kurtosis]: https://en.wikipedia.org/wiki/Kurtosis
/// [sampling bias]: https://en.wikipedia.org/wiki/Sampling_bias
pub struct RandomVariable<S: Weighable = f64, const M: usize = 2> {
	/// The estimated central moments of the sampled distribution.
	moments: CentralMoments<M>,
	/// Statistical moments whose types are preserved.
	inner: Cell<Option<Inner<S::Sample, S::Weight>>>,
}

/// The inner state containing the element count, smallest and largest elements
/// in the original type, and an additional field for the type of sampling.
#[derive(Clone, Copy, Eq, PartialEq)]
struct Inner<S, W> {
	/// The number of tabulated observations.
	count: NonZero<usize>,
	/// The smallest tabulated value.
	min: S,
	/// The largest tabulated value.
	max: S,
	/// Weight storage specific to the kind of sampling performed.
	weight: W,
}

/// A structure containing data related to the central statistical moments of
/// the distribution. The size depends on the number of moments requested.
///
/// Only up to four moments (mean, variance, skew, and kurtosis) are supported.
/// See [RandomVariable] for a more detailed description.
#[derive(Clone, PartialEq)]
pub struct CentralMoments<const N: usize>([Cell<f64>; N]);

impl<S: Weighable> RandomVariable<S> {
	/// Creates a new, empty random variable.
	pub const fn new<const M: usize>() -> RandomVariable<S, M> {
		RandomVariable {
			moments: CentralMoments([const { Cell::new(0.0) }; M]),
			inner: Cell::new(None),
		}
	}
}

impl<S: Weighable, const M: usize> RandomVariable<S, M> {
	/// Resets the stored statistical moments.
	pub fn reset(&self) {
		for m in &self.moments.0 {
			m.set(0.0);
		}
		self.inner.set(None);
	}

	/// Adds another value to the statistical collection.
	pub fn tabulate(&self, value: S)
	where
		CentralMoments<M>: Univariate,
	{
		assert!(self.try_tabulate(value).is_ok());
	}

	/// Adds a value to the statistical collection and returns an error if it
	/// fails.
	pub fn try_tabulate(&self, value: S) -> Result<(), S::Error>
	where
		CentralMoments<M>: Univariate,
	{
		value.sample(self)
	}

	/// Combines another equally typed random variable with this one, consuming
	/// both and returning a new random variable representing the aggregated
	/// statistical collection.
	pub fn join(self, other: Self) -> Self
	where
		CentralMoments<M>: Univariate,
	{
		let wl = self.weight();
		let wr = other.weight();

		// check if the left random variable contains values
		match self.inner.take() {
			Some(li) => {
				// check if the right random variable contains values
				match other.inner.into_inner() {
					Some(ri) => {
						// values in both left and right random variables
						let count = li
							.count
							.checked_add(ri.count.get())
							.expect("overflow in the number of tabulations detected");

						RandomVariable {
							moments: self.moments.combine(other.moments, wl, wr),
							inner: Cell::new(Some(Inner {
								count,
								weight: S::combine(li.weight, ri.weight),
								min: if li.min < ri.min { li.min } else { ri.min },
								max: if li.max > ri.max { li.max } else { ri.max },
							})),
						}
					}
					_ => {
						// no values in the right, but values in the left
						self.inner.set(Some(li));
						self
					}
				}
			}
			_ => {
				// no values in the left -> return the right
				other
			}
		}
	}

	/// Returns the number of tabulated observations.
	pub fn count(&self) -> usize {
		let inner = self.inner.take();
		let res = inner.as_ref().map_or(0, |inner| inner.count.get());
		self.inner.set(inner);
		res
	}

	/// Returns the smallest tabulated observation or `None` if no observations
	/// have occurred.
	pub fn min(&self) -> Option<S::Sample> {
		let inner = self.inner.take();
		let res = inner.as_ref().map(|inner| inner.min.clone());
		self.inner.set(inner);
		res
	}

	/// Returns the largest tabulated observation or `None` if no observations
	/// have occurred.
	pub fn max(&self) -> Option<S::Sample> {
		let inner = self.inner.take();
		let res = inner.as_ref().map(|inner| inner.max.clone());
		self.inner.set(inner);
		res
	}

	/// Returns the sum of the weights of the observations.
	pub fn weight(&self) -> f64 {
		S::weight(self)
	}

	/// Returns the weighted sum of the observations.
	pub fn sum(&self) -> S::Sample
	where
		CentralMoments<M>: Mean,
	{
		let weight = self.weight();

		if weight == 0.0 {
			return Sample::qualify(0.0);
		}

		Sample::qualify(self.moments.mean() * weight)
	}

	/// Returns the estimated [arithmetic mean] of the observations.
	///
	/// [arithmetic mean]: https://en.wikipedia.org/wiki/Arithmetic_mean
	pub fn mean(&self) -> f64
	where
		CentralMoments<M>: Mean,
	{
		// sample weight 0 is undefined
		if self.weight() == 0.0 {
			return f64::NAN;
		}

		// calculate the mean
		self.moments.mean()
	}

	/// Returns the [unbiased sample variance] of the observations under the
	/// assumption that the samples were taken from an infinite population.
	///
	/// # Note
	/// If you know that you have sampled the whole population, i.e., all
	/// possible samples, use [Self::population_variance] to determine the
	/// real (not estimated) variance.
	///
	/// [unbiased sample variance]: https://en.wikipedia.org/wiki/Variance#Unbiased_sample_variance
	pub fn variance(&self) -> f64
	where
		CentralMoments<M>: Variance,
	{
		let w = self.weight();

		// variance for sample weights under 1.0 is undefined
		if w <= 1.0 {
			return f64::NAN;
		}

		// very small variances are rounded to clean 0.0
		let m2 = self.moments.m2().get();
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return 0.0;
		}

		// calculate the variance
		self.moments.sample_variance(w)
	}

	/// Returns the [corrected sample standard deviation] of the observations
	/// under the assumption that the samples were taken from an infinite
	/// population.
	///
	/// # Note
	/// If you know that you have sampled the whole population, i.e., all
	/// possible samples, use [Self::population_std_dev] to determine the
	/// real (not estimated) standard deviation.
	///
	/// [corrected sample standard deviation]: https://en.wikipedia.org/wiki/Standard_deviation#Corrected_sample_standard_deviation
	#[cfg(any(feature = "std", feature = "libm"))]
	pub fn std_dev(&self) -> f64
	where
		CentralMoments<M>: Variance,
	{
		let w = self.weight();

		// a sample weight under 1.0 is undefined
		if w <= 1.0 {
			return f64::NAN;
		}

		// a very small sample variance is rounded to a clean 0.0
		let m2 = self.moments.m2().get();
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return 0.0;
		}

		// calculate the standard deviation
		sqrt(self.moments.sample_variance(w))
	}

	/// Returns the [skewness] of the observations under the assumption
	/// that the samples were taken from an infinite population.
	///
	/// # Note
	/// If you know that you have sampled the whole population, i.e., all
	/// possible samples, use [Self::population_skew] to determine the
	/// real (not estimated) skewness.
	///
	/// [skewness]: https://en.wikipedia.org/wiki/Skewness
	#[cfg(any(feature = "std", feature = "libm"))]
	pub fn skew(&self) -> f64
	where
		CentralMoments<M>: Skewness,
	{
		let w = self.weight();
		let m2 = self.moments.m2().get();

		// a sample weight under two or a very small variance is undefined
		if w <= 2.0 || m2 * m2 <= f64::EPSILON * self.mean() {
			return f64::NAN;
		}

		// calculate the skew
		self.moments.sample_skew(w)
	}

	/// Returns the [sample kurtosis] of the observations under the assumption
	/// that the samples were taken from an infinite population.
	///
	/// # Note
	/// If you know that you have sampled the whole population, i.e., all
	/// possible samples, use [Self::population_kurtosis] to determine the
	/// real (not estimated) kurtosis.
	///
	/// [sample kurtosis]: https://en.wikipedia.org/wiki/Kurtosis#Sample_kurtosis
	pub fn kurtosis(&self) -> f64
	where
		CentralMoments<M>: Kurtosis,
	{
		let w = self.weight();
		let m2 = self.moments.m2().get();

		// a sample weight under three or a very small variance is undefined
		if w <= 3.0 || m2 * m2 <= f64::EPSILON * self.mean() {
			return f64::NAN;
		}

		// calculate the kurtosis
		self.moments.sample_kurtosis(self.weight())
	}

	/// Returns the [excess kurtosis] of the observations under the assumption
	/// that the samples were taken from an infinite population.
	///
	/// # Note
	/// If you know that you have sampled the whole population, i.e., all
	/// possible samples, use [Self::population_excess_kurtosis] to determine
	/// the real (not estimated) kurtosis.
	///
	/// [excess kurtosis]: https://en.wikipedia.org/wiki/Kurtosis#Excess_kurtosis
	pub fn excess_kurtosis(&self) -> f64
	where
		CentralMoments<M>: Kurtosis,
	{
		let w = self.weight();
		let m2 = self.moments.m2().get();

		// a sample weight under three or a very small variance is undefined
		if w <= 3.0 || m2 * m2 <= f64::EPSILON * self.mean() {
			return f64::NAN;
		}

		// calculate the kurtosis
		self.moments.sample_excess_kurtosis(self.weight())
	}

	/// Returns the [population variance] of the observations under the
	/// assumption that the whole population has been sampled.
	///
	/// # Note
	/// This is likely not what you want since you will usually only be able
	/// to sample a limited subset of all possible observations and will
	/// therefore want to use [Self::variance] instead. However, both values
	/// will eventually converge to the same number.
	///
	/// [population variance]: https://en.wikipedia.org/wiki/Variance#Population_variance
	pub fn population_variance(&self) -> f64
	where
		CentralMoments<M>: Variance,
	{
		let w = self.weight();

		// a very small sample variance is rounded to a clean 0.0
		let m2 = self.moments.m2().get();
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return 0.0;
		}

		// calculate the variance
		self.moments.population_variance(w)
	}

	/// Returns the [population standard deviation] of the observations under
	/// the assumption that the whole population has been sampled.
	///
	/// # Note
	/// This is likely not what you want since you will usually only be able
	/// to sample a limited subset of all possible observations and will
	/// therefore want to use [Self::std_dev] instead. However, both values
	/// will eventually converge to the same number.
	///
	/// [population standard deviation]: https://en.wikipedia.org/wiki/Standard_deviation#Definition_of_population_values
	#[cfg(any(feature = "std", feature = "libm"))]
	pub fn population_std_dev(&self) -> f64
	where
		CentralMoments<M>: Variance,
	{
		let w = self.weight();

		// a very small sample variance is rounded to a clean 0.0
		let m2 = self.moments.m2().get();
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return 0.0;
		}

		// calculate the standard deviation
		sqrt(self.moments.population_variance(w))
	}

	/// Returns the [skewness] of the observations under the assumption that
	/// the whole population has been sampled.
	///
	/// # Note
	/// This is likely not what you want since you will usually only be able
	/// to sample a limited subset of all possible observations and will
	/// therefore want to use [Self::skew] instead. However, both values
	/// will eventually converge to the same number.
	///
	/// [skewness]: https://en.wikipedia.org/wiki/Skewness
	#[cfg(any(feature = "std", feature = "libm"))]
	pub fn population_skew(&self) -> f64
	where
		CentralMoments<M>: Skewness,
	{
		let w = self.weight();
		let m2 = self.moments.m2().get();

		// a very small variance is undefined
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return f64::NAN;
		}

		// calculate the skew
		self.moments.population_skew(w)
	}

	/// Returns the [kurtosis] of the observations under the assumption that
	/// the whole population has been sampled.
	///
	/// # Note
	/// This is likely not what you want since you will usually only be able
	/// to sample a limited subset of all possible observations and will
	/// therefore want to use [Self::kurtosis] instead. However, both values
	/// will eventually converge to the same number.
	///
	/// [kurtosis]: https://en.wikipedia.org/wiki/Kurtosis
	pub fn population_kurtosis(&self) -> f64
	where
		CentralMoments<M>: Kurtosis,
	{
		let w = self.weight();
		let m2 = self.moments.m2().get();

		// a very small variance is undefined
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return f64::NAN;
		}

		// calculate the kurtosis
		self.moments.population_kurtosis(w)
	}

	/// Returns the [excess kurtosis] of the observations under the assumption
	/// that the whole population has been sampled.
	///
	/// # Note
	/// The excess kurtosis of a population is just the [population kurtosis]
	/// minus 3. This is usually not what you want since you will usually only
	/// be able to sample a limited subset of all possible observations and will
	/// therefore want to use [Self::excess_kurtosis] instead. However, both
	/// values will eventually converge to the same number.
	///
	/// [excess kurtosis]: https://en.wikipedia.org/wiki/Kurtosis#Excess_kurtosis
	/// [population kurtosis]: Self::population_kurtosis
	pub fn population_excess_kurtosis(&self) -> f64
	where
		CentralMoments<M>: Kurtosis,
	{
		let w = self.weight();
		let m2 = self.moments.m2().get();

		// a very small variance is undefined
		if m2 * m2 <= f64::EPSILON * self.mean() {
			return f64::NAN;
		}

		// calculate the kurtosis
		self.moments.population_excess_kurtosis(w)
	}
}

// differentially weighted random variable
impl<T: Sample, S: Sample, const M: usize> RandomVariable<Utilized<T, S>, M>
where
	T: Add<Output = T> + Sub<Output = T>,
{
	/// Returns the currently tracked value within the most recent interval.
	pub fn tracked(&self) -> Option<S> {
		let inner = self.inner.take();
		let res = inner.as_ref().map(|inner| inner.weight.2.clone());
		self.inner.set(inner);
		res
	}

	/// Updates the temporal random variable's state to the current time without
	/// changing the tracked value and panics if the given time is before the
	/// previous tabulation.
	pub fn update(&self, time: T)
	where
		CentralMoments<M>: Univariate,
	{
		assert!(
			self.try_update(time).is_ok(),
			"updating the time to before the previous tabulation is not permitted"
		)
	}

	/// Updates the temporal random variable's state to the current time without
	/// changing the tracked value and returns an error if the given time is
	/// before the previous tabulation.
	pub fn try_update(&self, time: T) -> Result<(), Causality<T>>
	where
		CentralMoments<M>: Univariate,
	{
		let mut result = Ok(());

		if let Some(mut inner) = self.inner.take() {
			match time.partial_cmp(&inner.weight.1) {
				None | Some(Ordering::Less) => {
					result = Err(Causality {
						cause: inner.weight.1.clone(),
						effect: time,
					});
				}
				Some(Ordering::Equal) => {}
				Some(Ordering::Greater) => {
					// update the central moments
					self.moments.sample_w(
						inner.weight.2.quantify(),
						(time.clone() - inner.weight.1).quantify(),
						(time.clone() - inner.weight.0.clone()).quantify(),
					);

					// update the maximum
					inner.max = time.clone();

					// update the reference time
					inner.weight.1 = time;
				}
			}

			// restore the inner value
			self.inner.set(Some(inner));
		}
		// no current interval to update -> no update

		result
	}
}

/* *********************************** Weighted and Unweighted Tag Structures */

/// Helper-Trait for the kinds of value accepted for
/// [RandomVariables](RandomVariable).
pub trait Sample: Clone + PartialOrd {
	/// Converts the value into a floating point number, potentially stripping
	/// units of measurements.
	///
	/// The units of measurement may be recovered (albeit with a potential loss
	/// in precision) through the use of [qualify](Self::qualify).
	fn quantify(&self) -> f64;

	/// Converts a floating point number back into the sample type.
	///
	/// This method is allowed to produce values with reduced precision since
	/// its results are not used cumulatively. This method is used to
	/// recover units of measurement that were stripped earlier during the
	/// conversation via [quantify](Self::quantify).
	fn qualify(val: f64) -> Self;
}

/// Helper-trait unifying calculations for the weighted and unweighted case.
///
/// Since variance and join are computed differently in those cases, this trait
/// allows us to unify the interface for the user.
pub trait Weighable: Clone + Sized {
	/// The value-type that is being sampled.
	type Sample: Sample;
	/// The weight-type of the samples.
	type Weight: Clone;
	/// The error-type in case sampling fails.
	type Error;

	/// Adds this value to the collection of values represented by the random
	/// variable.
	fn sample<const M: usize>(self, rv: &RandomVariable<Self, M>) -> Result<(), Self::Error>
	where
		CentralMoments<M>: Univariate;

	/// Returns the weight of the random variable as a floating point number.
	fn weight<const M: usize>(rv: &RandomVariable<Self, M>) -> f64;

	/// Combines the weight of two random variables.
	fn combine(weight: Self::Weight, other: Self::Weight) -> Self::Weight;
}

impl<S: Sample> Weighable for S {
	type Sample = S;
	type Weight = ();
	type Error = Infallible;

	fn sample<const M: usize>(self, rv: &RandomVariable<Self, M>) -> Result<(), Self::Error>
	where
		CentralMoments<M>: Univariate,
	{
		rv.inner.set(Some(match rv.inner.replace(None) {
			Some(mut inner) => {
				// update the counter
				inner.count = inner
					.count
					.checked_add(1)
					.expect("overflow in the number of tabulations detected");

				// update minimum and maximum
				if inner.min > self {
					inner.min = self.clone();
				}
				if inner.max < self {
					inner.max = self.clone();
				}

				// update the central moments
				rv.moments.sample_u(self.quantify(), inner.count.get());

				inner
			}
			_ => {
				// initialize the sum of values and the inner structure
				rv.moments.sample_u(self.quantify(), 1);
				Inner {
					count: NonZero::new(1).unwrap(),
					min: self.clone(),
					max: self,
					weight: (),
				}
			}
		}));

		Ok(())
	}

	fn weight<const M: usize>(rv: &RandomVariable<Self, M>) -> f64 {
		rv.count() as f64
	}

	fn combine(_: Self::Weight, _: Self::Weight) -> Self::Weight {}
}

/// Represents a weighted sample.
#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
pub struct Weighted<S, W = f64>(
	/// The value.
	pub S,
	/// The weight.
	pub W,
);

impl<S: Sample, W: Sample + Add<Output = W>> Weighable for Weighted<S, W> {
	type Sample = S;
	type Weight = W;
	type Error = NegativeWeight<W>;

	fn sample<const M: usize>(self, rv: &RandomVariable<Self, M>) -> Result<(), Self::Error>
	where
		CentralMoments<M>: Univariate,
	{
		let Self(value, weight) = self;

		// quantify the weight
		let w = weight.quantify();

		// check the sign of the weight and return an error if it's negative
		if w.partial_cmp(&0.0).is_none_or(Ordering::is_lt) {
			return Err(NegativeWeight { weight });
		}

		rv.inner.set(Some(match rv.inner.replace(None) {
			Some(mut inner) => {
				// update count and weight
				inner.count = inner
					.count
					.checked_add(1)
					.expect("overflow in the number of tabulations detected");
				inner.weight = inner.weight + weight;

				// update minimum and maximum
				if inner.min > value {
					inner.min = value.clone();
				}
				if inner.max < value {
					inner.max = value.clone();
				}

				// update the central moments
				rv.moments
					.sample_w(value.quantify(), w, inner.weight.quantify());

				inner
			}
			_ => {
				// insert the first value
				rv.moments.sample_w(value.quantify(), w, w);
				Inner {
					count: const { NonZero::new(1).unwrap() },
					min: value.clone(),
					max: value.clone(),
					weight,
				}
			}
		}));

		Ok(())
	}

	fn weight<const M: usize>(rv: &RandomVariable<Self, M>) -> f64 {
		let inner = rv.inner.take();
		let w = inner.as_ref().map_or(0.0, |inner| inner.weight.quantify());
		rv.inner.set(inner);
		w
	}

	fn combine(weight: Self::Weight, other: Self::Weight) -> Self::Weight {
		weight + other
	}
}

/// Represents a time-weighted sample that can effectively be used to calculate
/// the utilization rate of model elements by tabulating (time-index, state)
/// pairs for every state change.
#[derive(Copy, Clone, Debug, Default, Eq, PartialEq)]
pub struct Utilized<T, S = bool>(
	/// The time index.
	pub T,
	/// The value.
	pub S,
);

impl<T: Sample, S: Sample> Weighable for Utilized<T, S>
where
	T: Add<Output = T> + Sub<Output = T>,
{
	type Sample = T;
	type Weight = (T, T, S);
	type Error = Causality<T>;

	fn sample<const M: usize>(self, rv: &RandomVariable<Self, M>) -> Result<(), Self::Error>
	where
		CentralMoments<M>: Univariate,
	{
		use core::mem::replace;
		let Self(time, value) = self;
		let mut result = Ok(());

		let inner = match rv.inner.take() {
			Some(mut inner) => {
				match time.partial_cmp(&inner.weight.1) {
					// ensure monotonicity in the time sequence
					None | Some(Ordering::Less) => {
						result = Err(Causality {
							cause: inner.weight.1.clone(),
							effect: time,
						});
					}
					Some(Ordering::Equal) => {
						// update the inner weight
						inner.weight.2 = value;
					}
					Some(Ordering::Greater) => {
						// update count and reference time
						inner.count = inner
							.count
							.checked_add(1)
							.expect("overflow in the number of tabulations detected");

						// update the maximum
						inner.max = time.clone();

						// update the central moments
						rv.moments.sample_w(
							replace(&mut inner.weight.2, value).quantify(),
							(time.clone() - inner.weight.1.clone()).quantify(),
							(time.clone() - inner.weight.0.clone()).quantify(),
						);

						// update the current time
						inner.weight.1 = time;
					}
				}

				inner
			}
			_ => {
				// open the first interval
				Inner {
					count: const { NonZero::new(1).unwrap() },
					min: time.clone(),
					max: time.clone(),
					weight: (time.clone(), time, value),
				}
			}
		};

		rv.inner.set(Some(inner));
		result
	}

	fn weight<const M: usize>(rv: &RandomVariable<Self, M>) -> f64 {
		let inner = rv.inner.take();
		let w = inner.as_ref().map_or(0.0, |inner| {
			(inner.weight.1.clone() - inner.weight.0.clone()).quantify()
		});
		rv.inner.set(inner);
		w
	}

	fn combine(weight: Self::Weight, other: Self::Weight) -> Self::Weight {
		(weight.0, weight.1 + (other.1 - other.0), other.2)
	}
}

/* **************************************************** Trait Implementations */

impl<S: Weighable, const M: usize> Default for RandomVariable<S, M> {
	fn default() -> Self {
		RandomVariable::<S>::new::<M>()
	}
}

impl<S: Weighable, const M: usize> Clone for RandomVariable<S, M> {
	fn clone(&self) -> Self {
		let inner = self.inner.take();
		self.inner.set(inner.clone());

		RandomVariable {
			moments: self.moments.clone(),
			inner: Cell::new(inner),
		}
	}
}

impl<S1, S2> PartialEq<RandomVariable<S2>> for RandomVariable<S1>
where
	S1: Weighable,
	S2: Weighable<Sample = S1::Sample>,
{
	fn eq(&self, other: &RandomVariable<S2>) -> bool {
		self.moments.eq(&other.moments)
	}
}

impl<S: Weighable> fmt::Debug for RandomVariable<S, 0>
where
	S::Sample: fmt::Debug,
{
	fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
		let mut dbs = f.debug_struct("RandomVariable");
		let len = self.count();

		dbs.field("n", &len);

		if let Some(inner) = self.inner.take() {
			dbs.field("min", &inner.min).field("max", &inner.max);
			self.inner.set(Some(inner));
		}

		dbs.finish()
	}
}

impl<S: Weighable> fmt::Debug for RandomVariable<S, 1>
where
	S::Sample: fmt::Debug,
{
	fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
		let mut dbs = f.debug_struct("RandomVariable");
		let len = self.count();

		dbs.field("n", &len);

		if let Some(inner) = self.inner.take() {
			dbs.field("min", &inner.min).field("max", &inner.max);
			self.inner.set(Some(inner));
			dbs.field("mean", &self.mean());
		}

		dbs.finish()
	}
}

impl<S: Weighable> fmt::Debug for RandomVariable<S, 2>
where
	S::Sample: fmt::Debug,
{
	fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
		let mut dbs = f.debug_struct("RandomVariable");
		let len = self.count();

		dbs.field("n", &len);

		if let Some(inner) = self.inner.take() {
			dbs.field("min", &inner.min).field("max", &inner.max);
			self.inner.set(Some(inner));
			dbs.field("mean", &self.mean());

			#[cfg(any(feature = "std", feature = "libm"))]
			dbs.field("sdev", &self.std_dev());

			#[cfg(all(not(feature = "std"), not(feature = "libm")))]
			dbs.field("variance", &self.variance());
		}

		dbs.finish()
	}
}

impl<S: Weighable> fmt::Debug for RandomVariable<S, 3>
where
	S::Sample: fmt::Debug,
{
	fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
		let mut dbs = f.debug_struct("RandomVariable");
		let len = self.count();

		dbs.field("n", &len);

		if let Some(inner) = self.inner.take() {
			dbs.field("min", &inner.min).field("max", &inner.max);
			self.inner.set(Some(inner));
			dbs.field("mean", &self.mean());

			#[cfg(any(feature = "std", feature = "libm"))]
			dbs.field("sdev", &self.std_dev())
				.field("skew", &self.skew());

			#[cfg(all(not(feature = "std"), not(feature = "libm")))]
			dbs.field("variance", &self.variance());
		}

		dbs.finish()
	}
}

impl<S: Weighable> fmt::Debug for RandomVariable<S, 4>
where
	S::Sample: fmt::Debug,
{
	fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
		let mut dbs = f.debug_struct("RandomVariable");
		let len = self.count();

		dbs.field("n", &len);

		if let Some(inner) = self.inner.take() {
			dbs.field("min", &inner.min).field("max", &inner.max);
			self.inner.set(Some(inner));
			dbs.field("mean", &self.mean());

			#[cfg(any(feature = "std", feature = "libm"))]
			dbs.field("sdev", &self.std_dev())
				.field("skew", &self.skew());

			#[cfg(all(not(feature = "std"), not(feature = "libm")))]
			dbs.field("variance", &self.variance());

			dbs.field("kurt", &self.kurtosis());
		}

		dbs.finish()
	}
}

macro_rules! impl_primitive_value {
	($($T:ident),*) => {$(
		impl Sample for $T {
			fn quantify(&self) -> f64 {
				*self as f64
			}

			fn qualify(val: f64) -> Self {
				val as $T
			}
		}
	)*}
}

macro_rules! impl_wrapped_value {
	($($T:ident < $V:ty >),*) => {$(
		impl Sample for $T<$V> {
			fn quantify(&self) -> f64 {
				self.0 as f64
			}

			fn qualify(val: f64) -> Self {
				$T(val as $V)
			}
		}
	)*}
}

impl_primitive_value!(
	i8, i16, i32, i64, i128, isize, u8, u16, u32, u64, u128, usize, f32, f64
);

impl_wrapped_value!(
	Wrapping<i8>,
	Wrapping<i16>,
	Wrapping<i32>,
	Wrapping<i64>,
	Wrapping<isize>,
	Wrapping<u8>,
	Wrapping<u16>,
	Wrapping<u32>,
	Wrapping<u64>,
	Wrapping<usize>
);

impl Sample for bool {
	fn quantify(&self) -> f64 {
		*self as u8 as f64
	}

	fn qualify(val: f64) -> Self {
		val != 0.0
	}
}

/* **************************************************** Central Moment Traits */

/// Helper-trait that unifies sampling of single values for different orders
/// of central moments.
#[doc(hidden)]
pub trait Univariate {
	/// Samples an unweighted value.
	fn sample_u(&self, x: f64, n: usize);

	/// Samples a weighted value.
	fn sample_w(&self, x: f64, w: f64, s: f64);

	/// Combines two central moments.
	fn combine(self, other: Self, wl: f64, wr: f64) -> Self;
}

/// Trait that marks support for the first central moment.
#[doc(hidden)]
pub trait Mean {
	/// Access to the statistical moment.
	fn m1(&self) -> &Cell<f64>;

	/// Calculates the mean based on the first moment.
	fn mean(&self) -> f64 {
		self.m1().get()
	}
}

/// Trait that marks support for the second central moment.
#[doc(hidden)]
pub trait Variance: Mean {
	/// Access to the statistical moment.
	fn m2(&self) -> &Cell<f64>;

	/// Calculates the variance under the assumption of a finite sample size out
	/// of an infinite population.
	fn sample_variance(&self, w: f64) -> f64 {
		let m2 = self.m2().get();
		m2 / (w - 1.0)
	}

	/// Calculates the variance under the assumption of a complete sample of the
	/// whole population.
	fn population_variance(&self, w: f64) -> f64 {
		let m2 = self.m2().get();
		m2 / w
	}
}

/// Trait that marks support for the third central moment.
#[doc(hidden)]
pub trait Skewness: Variance {
	/// Access to the statistical moment.
	fn m3(&self) -> &Cell<f64>;

	/// Calculates the skewness under the assumption of a finite sample size
	/// out of an infinite population.
	#[cfg(any(feature = "std", feature = "libm"))]
	fn sample_skew(&self, w: f64) -> f64 {
		let variance = self.sample_variance(w);
		let stddev = sqrt(variance);
		w / ((w - 1.0) * (w - 2.0)) * self.m3().get() / (variance * stddev)
	}

	/// Calculates the skewness under the assumption of a complete sample of the
	/// whole population.
	#[cfg(any(feature = "std", feature = "libm"))]
	fn population_skew(&self, w: f64) -> f64 {
		let variance = self.population_variance(w);
		let stddev = sqrt(variance);
		self.m3().get() / (w * variance * stddev)
	}
}

/// Trait that marks support for the fourth central moment.
#[doc(hidden)]
pub trait Kurtosis: Variance {
	/// Access to the statistical moment.
	fn m4(&self) -> &Cell<f64>;

	/// Calculates the kurtosis under the assumption of a finite sample size
	/// out of an infinite population.
	fn sample_kurtosis(&self, w: f64) -> f64 {
		let v = self.sample_variance(w);
		(w / (w - 1.0)) * ((w + 1.0) / (w - 2.0)) / (w - 3.0) * self.m4().get() / (v * v)
	}

	/// Calculates the excess kurtosis that returns 0.0 for distributions
	/// that spread their outliers like the normal distribution does,
	/// negative values for more localized and positive values for more spread
	/// out samples, relative to the normal distribution, under the assumption
	/// of a finite sample size out of an infinite population.
	fn sample_excess_kurtosis(&self, w: f64) -> f64 {
		self.sample_kurtosis(w) - 3.0 * ((w - 1.0) / (w - 2.0)) * ((w - 1.0) / (w - 3.0))
	}

	/// Calculates the kurtosis under the assumption of a complete sample of a
	/// finite population.
	fn population_kurtosis(&self, w: f64) -> f64 {
		let v = self.population_variance(w);
		self.m4().get() / w / (v * v)
	}

	/// Calculates the excess kurtosis that returns 0.0 for distributions
	/// that spread their outliers like the normal distribution does,
	/// negative values for more localized and positive values for more
	/// spread-out samples, relative to the normal distribution, under the
	/// assumption of a complete sample of the whole population.
	fn population_excess_kurtosis(&self, w: f64) -> f64 {
		self.population_kurtosis(w) - 3.0
	}
}

macro_rules! central_moments {
	(Mean for $CM:ident <$($N:literal),*>) => {
		$(impl Mean for $CM<$N> {
			#[inline]
			fn m1(&self) -> &Cell<f64> {
				&self.0[0]
			}
		})*
	};
	(Variance for $CM:ident <$($N:literal),*>) => {
		$(impl Variance for $CM<$N> {
			#[inline]
			fn m2(&self) -> &Cell<f64> {
				&self.0[1]
			}
		})*
	};
	(Skewness for $CM:ident <$($N:literal),*>) => {
		$(impl Skewness for $CM<$N> {
			#[inline]
			fn m3(&self) -> &Cell<f64> {
				&self.0[2]
			}
		})*
	};
	(Kurtosis for $CM:ident <$($N:literal),*>) => {
		$(impl Kurtosis for $CM<$N> {
			#[inline]
			fn m4(&self) -> &Cell<f64> {
				&self.0[3]
			}
		})*
	};
}

central_moments!(Mean for CentralMoments<1, 2, 3, 4>);
central_moments!(Variance for CentralMoments<2, 3, 4>);
central_moments!(Skewness for CentralMoments<3, 4>);
central_moments!(Kurtosis for CentralMoments<4>);

impl Univariate for CentralMoments<0> {
	fn sample_u(&self, _x: f64, _n: usize) {}

	fn sample_w(&self, _x: f64, _w: f64, _s: f64) {}

	fn combine(self, _other: Self, _wl: f64, _wr: f64) -> Self {
		self
	}
}

impl Univariate for CentralMoments<1> {
	fn sample_u(&self, x: f64, n: usize) {
		let m1 = self.m1();
		m1.set(m1.get() + (x - m1.get()) / n as f64);
	}

	fn sample_w(&self, x: f64, w: f64, s: f64) {
		let m1 = self.m1();
		m1.set(m1.get() + w * (x - m1.get()) / s);
	}

	fn combine(self, other: Self, wl: f64, wr: f64) -> Self {
		let sm1 = self.m1();
		let om1 = other.m1();

		sm1.set((wl * sm1.get() + wr * om1.get()) / (wl + wr));

		self
	}
}

impl Univariate for CentralMoments<2> {
	fn sample_u(&self, x: f64, n: usize) {
		let m1 = self.m1();
		let m2 = self.m2();
		let d = x - m1.get();

		m1.set(m1.get() + d / n as f64);
		m2.set(m2.get() + d * (x - m1.get()));
	}

	fn sample_w(&self, x: f64, w: f64, s: f64) {
		let m1 = self.m1();
		let m2 = self.m2();
		let d = w * (x - m1.get());

		m1.set(m1.get() + d / s);
		m2.set(m2.get() + d * (x - m1.get()));
	}

	fn combine(self, other: Self, wl: f64, wr: f64) -> Self {
		let sm1 = self.m1();
		let om1 = other.m1();
		let sm2 = self.m2();
		let om2 = other.m2();

		let m1d = om1.get() - sm1.get();
		let w = wl + wr;

		sm2.set(sm2.get() + om2.get() + m1d * m1d * wl * wr / w);
		sm1.set((wl * sm1.get() + wr * om1.get()) / w);

		self
	}
}

impl Univariate for CentralMoments<3> {
	fn sample_u(&self, x: f64, n: usize) {
		let m1 = self.m1();
		let m2 = self.m2();
		let m3 = self.m3();
		let d = x - m1.get();
		let s = n as f64;

		let a = d / s;
		m1.set(m1.get() + a);
		let b = x - m1.get();
		m3.set(m3.get() + a * (d * b * (s - 2.0) - 3.0 * m2.get()));
		m2.set(m2.get() + d * b);
	}

	fn sample_w(&self, x: f64, w: f64, s: f64) {
		let m1 = self.m1();
		let m2 = self.m2();
		let m3 = self.m3();
		let d = x - m1.get();

		let a = w * d / s;
		m1.set(m1.get() + a);
		let b = x - m1.get();
		m3.set(m3.get() + a * (d * b * (s - 2.0 * w) - 3.0 * m2.get()));
		m2.set(m2.get() + w * d * b);
	}

	fn combine(self, other: Self, wl: f64, wr: f64) -> Self {
		let sm1 = self.m1();
		let om1 = other.m1();
		let sm2 = self.m2();
		let om2 = other.m2();
		let sm3 = self.m3();
		let om3 = other.m3();

		let m1d = om1.get() - sm1.get();
		let m1d2 = m1d * m1d;
		let w = wl + wr;

		sm3.set(
			sm3.get()
				+ om3.get() + 3.0 * m1d * (wl * om2.get() - wr * sm2.get()) / w
				+ m1d2 * m1d * wl * wr * (wl - wr) / (w * w),
		);

		sm2.set(sm2.get() + om2.get() + m1d2 * wl * wr / w);
		sm1.set((wl * sm1.get() + wr * om1.get()) / w);

		self
	}
}

impl Univariate for CentralMoments<4> {
	fn sample_u(&self, x: f64, n: usize) {
		let m1 = self.m1();
		let m2 = self.m2();
		let m3 = self.m3();
		let m4 = self.m4();
		let d = x - m1.get();
		let s = n as f64;

		let a = d / s;
		m4.set(
			m4.get()
				+ a * (6.0 * a * m2.get() - 4.0 * m3.get()
					+ a * a * d * (s - 1.0) * (s * (s - 3.0) + 3.0)),
		);
		m1.set(m1.get() + a);
		let b = x - m1.get();
		m3.set(m3.get() + a * (d * b * (s - 2.0) - 3.0 * m2.get()));
		m2.set(m2.get() + d * b);
	}

	fn sample_w(&self, x: f64, w: f64, s: f64) {
		let m1 = self.m1();
		let m2 = self.m2();
		let m3 = self.m3();
		let m4 = self.m4();
		let d = x - m1.get();

		let a = d / s;
		let t = w * a;
		m4.set(
			m4.get()
				+ t * (6.0 * t * m2.get() - 4.0 * m3.get()
					+ a * a * d * (s - w) * (s * (s - 3.0 * w) + 3.0 * w * w)),
		);
		m1.set(m1.get() + t);
		let b = x - m1.get();
		m3.set(m3.get() + t * (d * b * (s - 2.0 * w) - 3.0 * m2.get()));
		m2.set(m2.get() + w * d * b);
	}

	fn combine(self, other: Self, wl: f64, wr: f64) -> Self {
		let sm1 = self.m1();
		let om1 = other.m1();
		let sm2 = self.m2();
		let om2 = other.m2();
		let sm3 = self.m3();
		let om3 = other.m3();
		let sm4 = self.m4();
		let om4 = other.m4();

		let m1d = om1.get() - sm1.get();
		let m1d2 = m1d * m1d;
		let w = wl + wr;
		let w2 = w * w;

		sm4.set(
			sm4.get()
				+ om4.get() + 4.0 * m1d * (wl * om3.get() - wr * sm3.get()) / w
				+ 6.0 * m1d2 * (wl * wl * om2.get() / w2 + wr * wr * sm2.get() / w2)
				+ m1d2 * m1d2 * wl * wr / w * (wl * wl - wl * wr + wr * wr) / w2,
		);

		sm3.set(
			sm3.get()
				+ om3.get() + 3.0 * m1d * (wl * om2.get() - wr * sm2.get()) / w
				+ m1d2 * m1d * wl * wr * (wl - wr) / w2,
		);

		sm2.set(sm2.get() + om2.get() + m1d2 * wl * wr / w);
		sm1.set((wl * sm1.get() + wr * om1.get()) / w);

		self
	}
}

/// Private `sqrt` implementation that falls back to `libm` in `no_std`
/// environments.
#[cfg(any(feature = "std", feature = "libm"))]
fn sqrt(val: f64) -> f64 {
	#[cfg(feature = "std")]
	{
		val.sqrt()
	}
	#[cfg(all(feature = "libm", not(feature = "std")))]
	{
		libm::sqrt(val)
	}
}

#[cfg(test)]
mod tests {
	use float_eq::float_eq;
	use prop::collection::vec;
	use proptest::prelude::*;

	use super::*;

	#[test]
	fn new_rv_is_empty() {
		let rv1 = RandomVariable::<i32>::default();
		let rv2 = RandomVariable::<Weighted<i32>>::default();
		let rv3 = RandomVariable::<Utilized<f32, i32>>::default();

		assert!(
			rv1.count() == 0
				&& rv1.min().is_none()
				&& rv1.max().is_none()
				&& rv1.mean().is_nan()
				&& rv1.variance().is_nan()
				&& rv1.std_dev().is_nan()
		);

		assert!(
			rv2.count() == 0
				&& rv2.min().is_none()
				&& rv2.max().is_none()
				&& rv2.mean().is_nan()
				&& rv2.variance().is_nan()
				&& rv2.std_dev().is_nan()
		);

		assert!(
			rv3.count() == 0
				&& rv3.min().is_none()
				&& rv3.max().is_none()
				&& rv3.mean().is_nan()
				&& rv3.variance().is_nan()
				&& rv3.std_dev().is_nan()
		);
	}

	#[test]
	#[should_panic]
	fn negative_weight_panics() {
		let rv = RandomVariable::new::<4>();

		for e in 0..100 {
			rv.tabulate(Weighted(e, 1.5));
		}

		rv.tabulate(Weighted(0, -1.0));
	}

	#[test]
	#[should_panic]
	fn negative_time_delta_panics() {
		let rv = RandomVariable::new::<4>();
		let mut time = 0;

		for e in 0..100 {
			rv.tabulate(Utilized(time, e));
			time += 1;
		}

		rv.tabulate(Utilized(time - 2, 0));
	}

	proptest! {
		#[test]
		#[cfg_attr(miri, ignore)]
		fn mean_calculated_correctly(elems in vec(0..1000, 2..100)) {
			let rv1 = RandomVariable::new::<4>();
			let n = elems.len();
			let mut s = 0;

			for elem in elems {
				rv1.tabulate(elem);
				s += elem;
			}

			let mean = s as f64 / n as f64;
			prop_assert!(
				float_eq!(
					rv1.mean(),
					mean,
					r2nd <= ((n as f64) * 0.5) * f64::EPSILON
				),
				"{:?} != {:?}", rv1.mean(), mean
			);
		}

		#[test]
		#[cfg_attr(miri, ignore)]
		fn join_with_empty_rv_unweighted(elems in vec(-100.0f64..100.0, 2..100)) {
			let rv1 = RandomVariable::default();

			for elem in elems {
				rv1.tabulate(elem);
			}

			prop_assert_eq!(&rv1, &rv1.clone().join(RandomVariable::new()));
			prop_assert_eq!(&rv1, &RandomVariable::new().join(rv1.clone()));
		}

		#[test]
		#[cfg_attr(miri, ignore)]
		fn join_with_empty_rv_weighted(elems in vec((0..100, 0.0f64..1000.0), 2..100)) {
			let rv1 = RandomVariable::default();

			for elem in elems {
				rv1.tabulate(Weighted(elem.0, elem.1));
			}

			prop_assert_eq!(&rv1, &rv1.clone().join(RandomVariable::default()));
			prop_assert_eq!(&rv1, &RandomVariable::default().join(rv1.clone()));
		}

		#[test]
		#[cfg_attr(miri, ignore)]
		fn join_with_empty_rv_temporal(elems in vec((0.0f64..1000.0, 0..100), 2..100)) {
			let rv1 = RandomVariable::default();
			let mut time = 0.0;

			for elem in elems {
				time += elem.0;
				rv1.tabulate(Utilized(time, elem.1));
			}

			prop_assert_eq!(&rv1, &rv1.clone().join(RandomVariable::default()));
			prop_assert_eq!(&rv1, &RandomVariable::default().join(rv1.clone()));
		}

		#[test]
		#[cfg_attr(miri, ignore)]
		fn equivalence_between_unweighted_and_weighted(elems in vec(-100.0..100.0, 3..100)) {
			let rv1 = RandomVariable::default();
			let rv2 = RandomVariable::default();

			for elem in elems {
				rv1.tabulate(elem);
				rv2.tabulate(Weighted(elem, 1));
			}

			prop_assert_eq!(rv1, rv2);
		}

		#[test]
		#[cfg_attr(miri, ignore)]
		fn equivalence_between_weighted_and_temporal(elems in vec((1..100, 1..100), 3..100)) {
			let rv1 = RandomVariable::default();
			let rv2 = RandomVariable::default();
			let mut t = 0;

			for (w,e) in elems {
				rv1.tabulate(Weighted(e, w));
				rv2.tabulate(Utilized(t, e));
				t += w;
			}
			rv2.update(t);

			prop_assert_eq!(rv1, rv2);
		}
	}
}