dsmc 0.18.0

//! System identification: ARX models by least squares / Kalman filter with input delay estimation
//! (`arx`), IV (`iv`), output-error models by the prediction-error method (`oe`), SRIVC (`srivc`,
//! also on multi-inertia plants), Levy and vector fitting (`frequency_response`), the
//! frequency-domain sample maximum likelihood estimator (`sml`), and Gaussian process regression
//! (`gpr`). The validation of the identified models is tested in `validation.rs`.
//!
//! The settings shared by the submodules are here: the tolerances on the identified coefficients
//! and the distance between two transfer functions.

mod arx;
mod frequency_response;
mod gpr;
mod iv;
mod oe;
mod sml;
mod srivc;

use dsmc::TransferFunction;

// Tolerances on normalized coefficients. All data are noise-free, so the measured errors are
// 1e-16 .. 1e-8; these leave a margin while still catching a broken identification.
const TOL_LSM_ARX: f64 = 1e-6;
const TOL_KF_ARX: f64 = 1e-6;
// Continuous-time parameters from the ARX models by the matched z-transform. It is not the inverse
// of the bilinear transform the data come from (frequency warping: relative (ω ts)^2 / 12 ≈ 1e-4
// for ω ≈ 32 rad/s, ts = 1 ms), so this bounds the method difference, not the identification.
const TOL_LSM_POLY: f64 = 1e-8;
const TOL_KF_POLY: f64 = 1e-8;
const TOL_LEVY: f64 = 1e-6;
const TOL_VF: f64 = 1e-8;

/// Max absolute difference between the coefficients of two transfer functions, after
/// normalizing both so that the leading denominator coefficient is 1.
fn tf_distance<D: dsmc::Domain>(a: &TransferFunction<f64, D>, b: &TransferFunction<f64, D>) -> f64 {
    let normalize = |tf: &TransferFunction<f64, D>| {
        let k = tf.denominator[0];
        let n: Vec<f64> = tf.numerator.iter().map(|c| c / k).collect();
        let d: Vec<f64> = tf.denominator.iter().map(|c| c / k).collect();
        (n, d)
    };
    let (na, da) = normalize(a);
    let (nb, db) = normalize(b);
    assert_eq!((na.len(), da.len()), (nb.len(), db.len()), "orders differ: {a:?} vs {b:?}");
    na.iter().zip(&nb).chain(da.iter().zip(&db)).map(|(x, y)| (x - y).abs()).fold(0.0, f64::max)
}