dsmc 0.11.0

//! Validation of an identified model on input / output data (e.g. a separate experiment).
//!
//! `Validation` holds the residual `ε = y - ŷ` of a model (output error or one-step prediction
//! error); the criteria and tests on it are in the submodules:
//!
//! | Module | Methods | Result |
//! | --- | --- | --- |
//! | `information_criteria` | `bic`, `aic`, `aicc` | |
//! | `whiteness` | `autocorrelation` | `WhitenessTest` |
//! | `cross_correlation` | `cross_correlation` | `CrossCorrelation` |
//! | `line_test` | `line_test` (periodic input) | `LineTest` |
//! | `coherence` | `coherence_test` (any input) | `CoherenceTest` |
//! | `frequency_response` | `frequency_response` (diagnostic) | `FrequencyResponseComparison` |
//!
//! (`spectra`: the segment spectra shared by `coherence` and `frequency_response`.)
//!
//! `check` runs several of them at once with one confidence level and pass / fail rules
//! (`Validation::check`, `Check`, `Report`):
//!
//! ```
//! use dsmc::tf;
//! use dsmc::system_identification::validation::{Check, Validation};
//!
//! let g = tf!("1000 / (s^2 + 20 s + 1000)");
//! let u: Vec<f64> = (0..4000).map(|k| ((k * 7919) % 101) as f64 / 50.0 - 1.0).collect();
//! let y0 = Validation::continuous(&g, 1e-3, &u, &vec![0.0; u.len()]).unwrap().simulated;
//! let noise = |k: usize| {
//!     // splitmix64: uniform in [-0.5, 0.5)
//!     let mut z = (k as u64 + 1).wrapping_mul(0x9E37_79B9_7F4A_7C15);
//!     z = (z ^ (z >> 30)).wrapping_mul(0xBF58_476D_1CE4_E5B9);
//!     z = (z ^ (z >> 27)).wrapping_mul(0x94D0_49BB_1331_11EB);
//!     ((z ^ (z >> 31)) >> 11) as f64 / (1u64 << 53) as f64 - 0.5
//! };
//! // Plant output with pseudo-random measurement noise
//! let y: Vec<f64> = y0.iter().enumerate().map(|(k, v)| v + 0.02 * noise(k)).collect();
//!
//! let report = Validation::continuous(&g, 1e-3, &u, &y)
//!     .unwrap()
//!     .check(
//!         &[
//!             Check::InformationCriteria { parameters: 3 },
//!             Check::CrossCorrelation { max_lag: 50 },
//!             Check::Coherence { segment_s: 0.5, excited: 1e-2 },
//!         ],
//!         0.99,
//!     )
//!     .unwrap();
//! println!("{report}");
//! assert!(report.passed());
//! ```

use std::ops::AddAssign;

use nalgebra::{ComplexField, RealField, Scalar};
use num_traits::Float;
use thiserror::Error;

use crate::{Discrete, Polynomial, SimulationError, TransferFunction, TransferFunctionWithDelay};
use crate::sampling::{samples_per_period, whole_samples, DurationError};
use crate::system_identification::arx::Arx;

mod check;
mod coherence;
mod information_criteria;
mod cross_correlation;
mod frequency_response;
mod line_test;
mod spectra;
mod whiteness;

pub use check::{Check, CheckKind, CheckResult, Report};
pub use coherence::CoherenceTest;
pub use cross_correlation::CrossCorrelation;
pub use frequency_response::FrequencyResponseComparison;
pub use line_test::LineTest;
pub use whiteness::WhitenessTest;

/// Errors of `Validation`.
#[derive(Clone, Debug, Error, PartialEq)]
pub enum ValidationError {
    #[error("input and output lengths differ: {input} vs {output}")]
    LengthMismatch { input: usize, output: usize },
    #[error("improper model: numerator degree {numerator} > denominator degree {denominator}")]
    Improper { numerator: usize, denominator: usize },
    #[error("zero denominator")]
    ZeroDenominator,
    #[error("{periods} whole period(s) evaluated; at least 2 are needed")]
    TooFewPeriods { periods: usize },
    #[error("line {line} is not between 1 and period / 2 = {half}")]
    InvalidLine { line: usize, half: usize },
    #[error("{segments} segment(s) of {segment_len} samples evaluated; at least 2 are needed")]
    TooFewSegments { segments: usize, segment_len: usize },
    #[error("duration {duration} s is negative or not finite, or sampling period {ts} s not positive and finite")]
    InvalidDuration { duration: f64, ts: f64 },
    #[error("duration {duration} s is not a whole number of sampling periods {ts} s")]
    FractionalDuration { duration: f64, ts: f64 },
    #[error("fundamental frequency {fundamental_frequency} Hz is not positive and finite")]
    InvalidFundamental { fundamental_frequency: f64, ts: f64 },
    #[error("the period of the fundamental frequency {fundamental_frequency} Hz is not a whole number of sampling periods {ts} s")]
    FractionalPeriod { fundamental_frequency: f64, ts: f64 },
    #[error(transparent)]
    Simulation(#[from] SimulationError),
}

/// Validation of a model on measured data by the residual `ε[k] = y[k] - ŷ[k]`, with `ŷ` either
///
/// - the output of the model simulated from rest with the measured input `u` (output error:
///   `continuous`, `discrete`), which validates the input-output model `G`; or
/// - the one-step prediction from the measured past outputs and inputs (prediction error:
///   `one_step_prediction`), which validates a model together with its noise model (ARX).
///
/// `bic`, `aic`, `autocorrelation` and `cross_correlation` apply to both. `line_test`,
/// `coherence_test` and `frequency_response` look at the input-output relation of `ŷ` and are
/// meant for the output error.
///
/// The residual is evaluated from sample `start` on (0 unless set by `evaluated_from`): to validate
/// on the second half of one experiment whose first half was used for the identification, simulate
/// the whole record (so that the state at the split is reproduced) and evaluate the second half.
///
/// Durations (`evaluated_from`, the segments of `coherence_test` and `frequency_response`) are
/// given in seconds and the fundamental frequency of `line_test` in Hz; they must be whole numbers
/// of the sampling period `ts` (the period `1 / f0` for the fundamental; relative tolerance
/// `max(1e-9, 4 eps)`). The lags of `autocorrelation` and `cross_correlation` are in samples.
#[derive(Clone, Debug)]
pub struct Validation<T> {
    /// Measured input `u[k]` over the whole record.
    pub input: Vec<T>,
    /// Model output `ŷ[k]` over the whole record (simulated, or one-step predicted).
    pub simulated: Vec<T>,
    /// Residual `ε[k] = y[k] - ŷ[k]` over the whole record.
    pub residual: Vec<T>,
    /// First evaluated sample.
    pub start: usize,
    /// Sampling period \[s\].
    pub ts: T,
}

impl<T: Float> Validation<T> {
    /// Continuous-time model `e^(-delay s) G(s)` (a `TransferFunctionWithDelay`, or a
    /// `TransferFunction` without delay) driven by `u` through a zero-order hold, sampled with
    /// period `ts` (`TransferFunctionWithDelay::simulate`: exact; the delay must be a whole number
    /// of sampling periods).
    pub fn continuous(model: impl Into<TransferFunctionWithDelay<T>>, ts: T, u: &[T], y: &[T]) -> Result<Self, ValidationError>
    where
        T: AddAssign + ComplexField + RealField,
    {
        if u.len() != y.len() {
            return Err(ValidationError::LengthMismatch { input: u.len(), output: y.len() });
        }
        let simulated = model.into().simulate(ts, u)?;
        Self::new(ts, u, simulated, y)
    }

    /// One-step prediction error of an ARX model (its parameters `arx.parameter`, orders and input
    /// delay; its history is not used): from rest,
    ///
    /// ```text
    /// ŷ[k|k-1] = Σ_i a_i y[k-i] + Σ_j b_j u[k-nk-j] = φ[k]ᵀ θ
    /// ```
    ///
    /// with the measured past outputs, as in the identification (`Arx::push`, `Arx::predict`).
    /// The ARX model `A(z) y = B(z) u + e` assumes a white equation error `e`, i.e. noise
    /// `e / A(z)` on the output: the model (with this noise model) is right only if the prediction
    /// error `ε = A(z) y - B(z) u` is white (`autocorrelation`) and independent of the input
    /// (`cross_correlation`). With white noise on the output instead (`y = G u + v`), even the true
    /// `A`, `B` leave the colored error `A(z) v`: the ARX noise model does not fit, and the IV
    /// method (validated by the output error, `discrete`) is the one to use.
    ///
    /// `ts` \[s\] is the sampling period of the data (for the durations of the tests).
    pub fn one_step_prediction(arx: &Arx<T>, ts: T, u: &[T], y: &[T]) -> Result<Self, ValidationError>
    where
        T: AddAssign + Scalar,
    {
        if u.len() != y.len() {
            return Err(ValidationError::LengthMismatch { input: u.len(), output: y.len() });
        }
        let mut model = arx.clone();
        model.clear_history();
        let predicted = (0..u.len())
            .map(|k| {
                model.push(u[k], if k > 0 { y[k - 1] } else { T::zero() });
                model.predict()
            })
            .collect();
        Self::new(ts, u, predicted, y)
    }

    /// Discrete-time model `G(z)` driven by `u`, sampled with period `ts` \[s\] (for the durations
    /// of the tests).
    pub fn discrete(model: &TransferFunction<T, Discrete>, ts: T, u: &[T], y: &[T]) -> Result<Self, ValidationError> {
        let (numer, denom) = monic(&model.numerator, &model.denominator)?;
        // G = z^-r (Σ n_i z^-i) / (1 + Σ d_i z^-i), r = deg D - deg N
        let r = denom.len() - numer.len();
        let mut simulated: Vec<T> = Vec::with_capacity(u.len());
        for k in 0..u.len() {
            let mut yk = T::zero();
            for (i, &c) in numer.iter().enumerate() {
                if k >= r + i {
                    yk = yk + c * u[k - r - i];
                }
            }
            for (i, &c) in denom.iter().enumerate().skip(1) {
                if k >= i {
                    yk = yk - c * simulated[k - i];
                }
            }
            simulated.push(yk);
        }
        Self::new(ts, u, simulated, y)
    }

    fn new(ts: T, u: &[T], simulated: Vec<T>, y: &[T]) -> Result<Self, ValidationError> {
        if u.len() != y.len() {
            return Err(ValidationError::LengthMismatch { input: u.len(), output: y.len() });
        }
        let residual = y.iter().zip(&simulated).map(|(&a, &b)| a - b).collect();
        // Checked here, so that the durations of the tests can only be wrong themselves
        whole_samples(T::zero(), ts).map_err(|_| ValidationError::InvalidDuration { duration: 0.0, ts: ts.to_f64().unwrap_or(f64::NAN) })?;
        Ok(Self { input: u.to_vec(), simulated, residual, start: 0, ts })
    }

    /// Evaluate the residual from the time `from_s` \[s\] on (sample `from_s / ts`), e.g. after the
    /// transient from rest (one period of a periodic input).
    pub fn evaluated_from(mut self, from_s: T) -> Result<Self, ValidationError> {
        self.start = self.duration_samples(from_s)?;
        Ok(self)
    }

    /// Evaluate the residual from sample `start` on (for the options in samples of SRIVC).
    pub(crate) fn evaluated_from_sample(mut self, start: usize) -> Self {
        self.start = start;
        self
    }

    /// Period `1 / (fundamental_frequency ts)` as a whole number of samples.
    pub(crate) fn period_samples(&self, fundamental_frequency: T) -> Result<usize, ValidationError> {
        samples_per_period(fundamental_frequency, self.ts).map_err(|kind| {
            let (fundamental_frequency, ts) = (fundamental_frequency.to_f64().unwrap_or(f64::NAN), self.ts.to_f64().unwrap_or(f64::NAN));
            match kind {
                DurationError::Invalid => ValidationError::InvalidFundamental { fundamental_frequency, ts },
                DurationError::Fractional => ValidationError::FractionalPeriod { fundamental_frequency, ts },
            }
        })
    }

    /// `duration / ts` as a whole number of samples.
    pub(crate) fn duration_samples(&self, duration: T) -> Result<usize, ValidationError> {
        whole_samples(duration, self.ts).map_err(|kind| {
            let (duration, ts) = (duration.to_f64().unwrap_or(f64::NAN), self.ts.to_f64().unwrap_or(f64::NAN));
            match kind {
                DurationError::Invalid => ValidationError::InvalidDuration { duration, ts },
                DurationError::Fractional => ValidationError::FractionalDuration { duration, ts },
            }
        })
    }

    /// Evaluated residual `ε[start..]`.
    pub fn evaluated(&self) -> &[T] {
        &self.residual[self.start.min(self.residual.len())..]
    }

    /// Number of evaluated samples `N`.
    pub fn samples(&self) -> usize {
        self.evaluated().len()
    }

    /// Mean squared residual `V = Σ ε^2 / N`.
    pub fn mse(&self) -> T {
        let e = self.evaluated();
        e.iter().fold(T::zero(), |acc, &v| acc + v * v) / T::from(e.len()).unwrap()
    }
}

/// Numerator and denominator without leading zeros, divided by the leading denominator coefficient.
fn monic<T: Float>(numerator: &Polynomial<T>, denominator: &Polynomial<T>) -> Result<(Vec<T>, Vec<T>), ValidationError> {
    let trim = |p: &Polynomial<T>| p.iter().copied().skip_while(|c| c.is_zero()).collect::<Vec<T>>();
    let (numer, denom) = (trim(numerator), trim(denominator));
    let Some(&lead) = denom.first() else {
        return Err(ValidationError::ZeroDenominator);
    };
    if numer.len() > denom.len() {
        return Err(ValidationError::Improper { numerator: numer.len() - 1, denominator: denom.len() - 1 });
    }
    let numer = if numer.is_empty() { vec![T::zero()] } else { numer };
    Ok((numer.iter().map(|&c| c / lead).collect(), denom.iter().map(|&c| c / lead).collect()))
}