au 0.10.0

//! # Continuous time linear system
//!
//! The time evolution of the system is performed through ODE (ordinary
//! differential equation) [solvers](../solver/index.html).

use nalgebra::{ComplexField, DVector, RealField};
use num_traits::Float;

use crate::{
    enums::Continuous,
    linear_system::{
        solver::{Order, Radau, Rk, Rkf45},
        Equilibrium, SsGen,
    },
    units::Seconds,
};

/// State-space representation of continuous time linear system
pub type Ss<T> = SsGen<T, Continuous>;

/// Implementation of the methods for the state-space
impl<T: ComplexField> Ss<T> {
    /// Calculate the equilibrium point for continuous time systems,
    /// given the input condition
    /// ```text
    /// x = - A^-1 * B * u
    /// y = - (C * A^-1 * B + D) * u
    /// ```
    ///
    /// # Arguments
    ///
    /// * `u` - Input vector
    ///
    /// # Example
    ///
    /// ```
    /// use au::Ss;
    /// let a = [-1., 1., -1., 0.25];
    /// let b = [1., 0.25];
    /// let c = [0., 1., -1., 1.];
    /// let d = [0., 1.];
    ///
    /// let sys = Ss::new_from_slice(2, 1, 2, &a, &b, &c, &d);
    /// let u = 0.0;
    /// let eq = sys.equilibrium(&[u]).unwrap();
    /// assert_eq!((0., 0.), (eq.x()[0], eq.y()[0]));
    /// ```
    pub fn equilibrium(&self, u: &[T]) -> Option<Equilibrium<T>> {
        assert_eq!(u.len(), self.b.ncols(), "Wrong number of inputs.");
        let u = DVector::from_row_slice(u);
        // 0 = A*x + B*u
        let bu = -&self.b * &u;
        let lu = &self.a.clone().lu();
        // A*x = -B*u
        let x = lu.solve(&bu)?;
        // y = C*x + D*u
        let y = &self.c * &x + &self.d * u;
        Some(Equilibrium::new(x, y))
    }
}

/// Implementation of the methods for the state-space
impl<T: ComplexField + Float + RealField> Ss<T> {
    /// System stability. Checks if all A matrix eigenvalues (poles) are negative.
    ///
    /// # Example
    ///
    /// ```
    /// use au::Ss;
    /// let sys = Ss::new_from_slice(2, 1, 1, &[-2., 0., 3., -7.], &[1., 3.], &[-1., 0.5], &[0.1]);
    /// assert!(sys.is_stable());
    /// ```
    #[must_use]
    pub fn is_stable(&self) -> bool {
        self.poles().iter().all(|p| p.re.is_negative())
    }
}

/// Implementation of the methods for the state-space
impl Ss<f64> {
    /// Time evolution for the given input, using Runge-Kutta second order method
    ///
    /// # Arguments
    ///
    /// * `u` - input function returning a vector (column mayor)
    /// * `x0` - initial state (column mayor)
    /// * `h` - integration time interval
    /// * `n` - integration steps
    pub fn rk2<F>(&self, u: F, x0: &[f64], h: Seconds<f64>, n: usize) -> Rk<F, f64>
    where
        F: Fn(Seconds<f64>) -> Vec<f64>,
    {
        Rk::new(self, u, x0, h, n, Order::Rk2)
    }

    /// Time evolution for the given input, using Runge-Kutta fourth order method
    ///
    /// # Arguments
    ///
    /// * `u` - input function returning a vector (column mayor)
    /// * `x0` - initial state (column mayor)
    /// * `h` - integration time interval
    /// * `n` - integration steps
    pub fn rk4<F>(&self, u: F, x0: &[f64], h: Seconds<f64>, n: usize) -> Rk<F, f64>
    where
        F: Fn(Seconds<f64>) -> Vec<f64>,
    {
        Rk::new(self, u, x0, h, n, Order::Rk4)
    }

    /// Runge-Kutta-Fehlberg 45 with adaptive step for time evolution.
    ///
    /// # Arguments
    ///
    /// * `u` - input function returning a vector (column vector)
    /// * `x0` - initial state (column vector)
    /// * `h` - integration time interval
    /// * `limit` - time evaluation limit
    /// * `tol` - error tolerance
    pub fn rkf45<F>(
        &self,
        u: F,
        x0: &[f64],
        h: Seconds<f64>,
        limit: Seconds<f64>,
        tol: f64,
    ) -> Rkf45<F, f64>
    where
        F: Fn(Seconds<f64>) -> Vec<f64>,
    {
        Rkf45::new(self, u, x0, h, limit, tol)
    }

    /// Radau of order 3 with 2 steps method for time evolution.
    ///
    /// # Arguments
    ///
    /// * `u` - input function returning a vector (column vector)
    /// * `x0` - initial state (column vector)
    /// * `h` - integration time interval
    /// * `n` - integration steps
    /// * `tol` - error tolerance
    pub fn radau<F>(&self, u: F, x0: &[f64], h: Seconds<f64>, n: usize, tol: f64) -> Radau<F, f64>
    where
        F: Fn(Seconds<f64>) -> Vec<f64>,
    {
        Radau::new(self, u, x0, h, n, tol)
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[allow(clippy::many_single_char_names)]
    #[test]
    fn equilibrium() {
        let a = [-1., 1., -1., 0.25];
        let b = [1., 0.25];
        let c = [0., 1., -1., 1.];
        let d = [0., 1.];

        let sys = Ss::new_from_slice(2, 1, 2, &a, &b, &c, &d);
        let u = 0.0;
        let eq = sys.equilibrium(&[u]).unwrap();
        assert_eq!((0., 0.), (eq.x()[0], eq.y()[0]));
        assert!(!format!("{}", eq).is_empty());
    }

    #[test]
    fn stability() {
        let eig1 = -2.;
        let eig2 = -7.;
        let sys = Ss::new_from_slice(
            2,
            1,
            1,
            &[eig1, 0., 3., eig2],
            &[1., 3.],
            &[-1., 0.5],
            &[0.1],
        );
        assert!(sys.is_stable())
    }

    #[test]
    fn new_rk2() {
        let a = [-1., 1., -1., 0.25];
        let b = [1., 0.25];
        let c = [0., 1.];
        let d = [0.];
        let sys = Ss::new_from_slice(2, 1, 1, &a, &b, &c, &d);
        let iter = sys.rk2(|_| vec![1.], &[0., 0.], Seconds(0.1), 30);
        assert_eq!(31, iter.count());
    }

    #[test]
    fn new_rk4() {
        let a = [-1., 1., -1., 0.25];
        let b = [1., 0.25];
        let c = [0., 1.];
        let d = [0.];
        let sys = Ss::new_from_slice(2, 1, 1, &a, &b, &c, &d);
        let iter = sys.rk4(|_| vec![1.], &[0., 0.], Seconds(0.1), 30);
        assert_eq!(31, iter.count());
    }

    #[test]
    fn new_rkf45() {
        let a = [-1., 1., -1., 0.25];
        let b = [1., 0.25];
        let c = [0., 1.];
        let d = [0.];
        let sys = Ss::new_from_slice(2, 1, 1, &a, &b, &c, &d);
        let iter = sys.rkf45(|_| vec![1.], &[0., 0.], Seconds(0.1), Seconds(2.), 1e-5);
        assert_relative_eq!(2., iter.last().unwrap().time().0, max_relative = 0.01);
    }

    #[test]
    fn new_radau() {
        let a = [-1., 1., -1., 0.25];
        let b = [1., 0.25];
        let c = [0., 1.];
        let d = [0.];
        let sys = Ss::new_from_slice(2, 1, 1, &a, &b, &c, &d);
        let iter = sys.radau(|_| vec![1.], &[0., 0.], Seconds(0.1), 30, 1e-5);
        assert_eq!(31, iter.count());
    }
}