automatica 1.0.0

//! Rational functions
//!
//! Ratio whose numerator and denominator are polynomials.
//!
//! ```text
//!        b_n*x^n + b_(n-1)*x^(n-1) + ... + b_1*x + b_0
//! f(x) = ---------------------------------------------
//!        a_m*x^m + a_(m-1)*x^(m-1) + ... + a_1*x + a_0
//! ```

use std::{
    fmt,
    fmt::{Debug, Display, Formatter},
    ops::{Add, Div, Mul},
};

use polynomen::Poly;

use crate::{
    wrappers::{PWrapper, PP},
    One, Polynomial, Zero,
};

mod arithmetic;

/// Rational function
#[derive(Clone, Debug, PartialEq)]
pub struct Rf<T> {
    /// Rational function numerator
    num: PP<T>,
    /// Rational function denominator
    den: PP<T>,
}

impl<T> Rf<T>
where
    T: Clone + PartialEq + Zero,
{
    /// Create a new rational function given its numerator and denominator
    ///
    /// # Arguments
    ///
    /// * `num` - Rational function numerator
    /// * `den` - Rational function denominator
    ///
    /// # Example
    /// ```
    /// use automatica::Rf;
    /// let rf = Rf::new([1., 2.], [-4., 6., -2.]);
    /// ```
    #[must_use]
    pub fn new<R, S>(num: R, den: S) -> Self
    where
        R: AsRef<[T]>,
        S: AsRef<[T]>,
    {
        Self {
            num: PP::new_from_coeffs(num.as_ref()),
            den: PP::new_from_coeffs(den.as_ref()),
        }
    }

    #[must_use]
    pub(crate) fn new_pp(num: PP<T>, den: PP<T>) -> Self {
        Self { num, den }
    }

    #[must_use]
    pub(crate) fn new_from_poly(num: Poly<PWrapper<T>>, den: Poly<PWrapper<T>>) -> Self {
        Self {
            num: PP(num),
            den: PP(den),
        }
    }

    /// Extract rational function numerator
    ///
    /// # Example
    /// ```
    /// use automatica::{Polynomial, Rf};
    /// let num = [1., 2.];
    /// let rf = Rf::new(num.clone(), [-4., 6., -2.]);
    /// assert_eq!(num, rf.num().as_slice());
    /// ```
    #[must_use]
    pub fn num(&self) -> impl Polynomial<T> {
        self.num.to_unwrapped_vec()
    }

    #[must_use]
    pub(crate) fn num_int(&self) -> &PP<T> {
        &self.num
    }

    /// Extract rational function denominator
    ///
    /// # Example
    /// ```
    /// use automatica::{Polynomial, Rf};
    /// let den = [-4., 6., -2.];
    /// let rf = Rf::new([1., 2.], den.clone());
    /// assert_eq!(den, rf.den().as_slice());
    /// ```
    #[must_use]
    pub fn den(&self) -> impl Polynomial<T> {
        self.den.to_unwrapped_vec()
    }

    #[must_use]
    pub(crate) fn den_int(&self) -> &PP<T> {
        &self.den
    }
}

impl<T: Clone + PartialEq + Zero> Rf<T> {
    /// Calculate the relative degree between denominator and numerator.
    ///
    /// # Example
    /// ```
    /// use automatica::{Inv, Rf};
    /// let rf = Rf::new([1., 2.], [-4., 6., -2.]);
    /// let expected = rf.relative_degree();
    /// assert_eq!(expected, 1);
    /// assert_eq!(rf.inv().relative_degree(), -1);
    /// ```
    #[must_use]
    #[allow(clippy::cast_possible_wrap, clippy::cast_possible_truncation)]
    pub fn relative_degree(&self) -> i32 {
        match (self.den.0.degree(), self.num.0.degree()) {
            (Some(d), Some(n)) => d as i32 - n as i32,
            (Some(d), None) => d as i32,
            (None, Some(n)) => -(n as i32),
            _ => 0,
        }
    }
}

macro_rules! rf_roots {
    ($ty:ty) => {
        impl Rf<$ty> {
            /// Calculate the poles of the rational function
            #[must_use]
            pub fn real_poles(&self) -> Option<Vec<$ty>> {
                self.den.real_roots()
            }

            /// Calculate the poles of the rational function
            #[must_use]
            pub fn complex_poles(&self) -> Vec<($ty, $ty)> {
                self.den.complex_roots()
            }

            /// Calculate the zeros of the rational function
            #[must_use]
            pub fn real_zeros(&self) -> Option<Vec<$ty>> {
                self.num.real_roots()
            }

            /// Calculate the zeros of the rational function
            #[must_use]
            pub fn complex_zeros(&self) -> Vec<($ty, $ty)> {
                self.num.complex_roots()
            }
        }
    };
}
rf_roots!(f32);
rf_roots!(f64);

impl<T> Rf<T>
where
    T: Clone + Div<Output = T> + One + PartialEq + Zero,
{
    /// Normalization of rational function. If the denominator is zero the same
    /// rational function is returned.
    ///
    /// from:
    /// ```text
    ///        b_n*z^n + b_(n-1)*z^(n-1) + ... + b_1*z + b_0
    /// G(z) = ---------------------------------------------
    ///        a_n*z^n + a_(n-1)*z^(n-1) + ... + a_1*z + a_0
    /// ```
    /// to:
    /// ```text
    ///        b'_n*z^n + b'_(n-1)*z^(n-1) + ... + b'_1*z + b'_0
    /// G(z) = -------------------------------------------------
    ///          z^n + a'_(n-1)*z^(n-1) + ... + a'_1*z + a'_0
    /// ```
    ///
    /// # Example
    /// ```
    /// use automatica::Rf;
    /// let rf = Rf::new([1., 2.], [-4., 6., -2.]);
    /// let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
    /// assert_eq!(expected, rf.normalize());
    /// ```
    #[must_use]
    pub fn normalize(&self) -> Self {
        if self.den.0.is_zero() {
            return self.clone();
        }
        let (den, an) = self.den.monic();
        let num = PP(&self.num.0 / an);
        Self::new_pp(num, PP(den))
    }

    /// In place normalization of rational function. If the denominator is zero
    /// no operation is done.
    ///
    /// from:
    /// ```text
    ///        b_n*z^n + b_(n-1)*z^(n-1) + ... + b_1*z + b_0
    /// G(z) = ---------------------------------------------
    ///        a_n*z^n + a_(n-1)*z^(n-1) + ... + a_1*z + a_0
    /// ```
    /// to:
    /// ```text
    ///        b'_n*z^n + b'_(n-1)*z^(n-1) + ... + b'_1*z + b'_0
    /// G(z) = -------------------------------------------------
    ///          z^n + a'_(n-1)*z^(n-1) + ... + a'_1*z + a'_0
    /// ```
    ///
    /// # Example
    /// ```
    /// use automatica::Rf;
    /// let mut rf = Rf::new([1., 2.], [-4., 6., -2.]);
    /// rf.normalize_mut();
    /// let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
    /// assert_eq!(expected, rf);
    /// ```
    pub fn normalize_mut(&mut self) {
        if self.den.0.is_zero() {
            return;
        }
        let an = self.den.0.monic_mut();
        self.num.0.div_mut(&an);
    }
}

impl<T> Rf<T>
where
    T: Clone + PartialEq,
{
    /// Evaluate the rational function.
    ///
    /// # Arguments
    ///
    /// * `s` - Value at which the rational function is evaluated.
    ///
    /// # Example
    /// ```
    /// use automatica::{Complex as C, Rf};
    /// let rf = Rf::new([1., 2., 3.], [-4., -3., 1.]);
    /// assert_eq!(-8.5, rf.eval_by_val(3.));
    /// assert_eq!(C::new(0.64, -0.98), rf.eval_by_val(C::new(0., 2.0)));
    /// ```
    pub fn eval_by_val<N>(&self, s: N) -> N
    where
        N: Add<T, Output = N> + Clone + Div<Output = N> + Mul<Output = N> + Zero,
    {
        let res = self.num.0.eval_by_val(PWrapper(s.clone())) / self.den.0.eval_by_val(PWrapper(s));
        res.0
    }
}

impl<T> Rf<T>
where
    T: Clone + PartialEq,
{
    /// Evaluate the rational function.
    ///
    /// # Arguments
    ///
    /// * `s` - Value at which the rational function is evaluated.
    ///
    /// # Example
    /// ```
    /// use automatica::{Complex as C, Rf};
    /// let rf = Rf::new([1., 2., 3.], [-4., -3., 1.]);
    /// assert_eq!(-8.5, rf.eval(&3.));
    /// assert_eq!(C::new(0.64, -0.98), rf.eval(&C::new(0., 2.0)));
    /// ```
    pub fn eval<N>(&self, s: &N) -> N
    where
        N: Add<T, Output = N> + Clone + Div<Output = N> + Mul<N, Output = N> + Zero,
    {
        // let x = PWrapper(s.clone());
        // let res = self.num.0.eval(&x) / self.den.0.eval(&x);
        let res = self.num.0.eval_by_val(PWrapper(s.clone()))
            / self.den.0.eval_by_val(PWrapper(s.clone()));
        res.0
    }
}

/// Implementation of rational function printing
impl<T> Display for Rf<T>
where
    T: Clone + Display + One + PartialEq + PartialOrd + Zero,
{
    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
        let (s_num, s_den) = if let Some(precision) = f.precision() {
            let num = format!("{poly:.prec$}", poly = self.num.0, prec = precision);
            let den = format!("{poly:.prec$}", poly = self.den.0, prec = precision);
            (num, den)
        } else {
            let num = format!("{}", self.num.0);
            let den = format!("{}", self.den.0);
            (num, den)
        };
        let length = s_num.len().max(s_den.len());
        let dash = "\u{2500}".repeat(length);
        write!(f, "{}\n{}\n{}", s_num, dash, s_den)
    }
}

#[cfg(test)]
mod tests {
    use crate::{Complex, Inv};

    use super::*;

    #[test]
    #[allow(clippy::float_cmp)]
    fn rational_function_creation() {
        let num = [1., 2., 3.];
        let den = [-4.2, -3.12, 0.0012];
        let rf = Rf::new(num, den);
        assert_eq!(num, rf.num().as_slice());
        assert_eq!(den, rf.den().as_slice());
    }

    #[test]
    fn relative_degree() {
        let rf = Rf::new([1., 2.], [-4., 6., -2.]);
        let expected = rf.relative_degree();
        assert_eq!(expected, 1);
        assert_eq!(rf.inv().relative_degree(), -1);
        assert_eq!(-1, Rf::new([1., 1.], Poly::zero()).relative_degree());
        assert_eq!(1, Rf::new(Poly::zero(), [1., 1.]).relative_degree());
        assert_eq!(
            0,
            Rf::<f32>::new(Poly::zero(), Poly::zero()).relative_degree()
        );
    }

    #[test]
    fn evaluation() {
        let rf = Rf::new([-0.75, 0.25], [0.75, 0.75, 1.]);
        let res = rf.eval(&Complex::new(0., 0.9));
        assert_abs_diff_eq!(0.429, res.re, epsilon = 0.001);
        assert_abs_diff_eq!(1.073, res.im, epsilon = 0.001);
    }

    #[test]
    fn evaluation_by_value() {
        let rf = Rf::new([-0.75, 0.25], [0.75, 0.75, 1.]);
        let res1 = rf.eval(&Complex::new(0., 0.9));
        let res2 = rf.eval_by_val(Complex::new(0., 0.9));
        assert_eq!(res1, res2);
    }

    #[test]
    fn poles() {
        let rf = Rf::new([1.0_f64], [6., -5., 1.]);
        assert_eq!(Some(vec![2., 3.]), rf.real_poles());
    }

    #[test]
    fn complex_poles() {
        let rf = Rf::new([1.0_f32], [10., -6., 1.]);
        assert_eq!(
            vec![Complex::new(3., -1.), Complex::new(3., 1.)],
            rf.complex_poles()
        );
    }

    #[test]
    fn zeros() {
        let rf = Rf::new([1.0_f64], [6., -5., 1.]);
        assert_eq!(None, rf.real_zeros());
    }

    #[test]
    fn complex_zeros() {
        let rf = Rf::new([3.25_f32, 3., 1.], [10., -3., 1.]);
        assert_eq!(
            vec![Complex::new(-1.5, -1.), Complex::new(-1.5, 1.)],
            rf.complex_zeros()
        );
    }

    #[test]
    fn print() {
        let rf = Rf::new(Poly::one(), Poly::new_from_roots(&[-1.]));
        assert_eq!(
            "1\n\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\n1 +1s",
            format!("{}", rf)
        );

        let rf2 = Rf::new([1.123], [0.987, -1.321]);
        assert_eq!(
            "1.12\n\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\u{2500}\n0.99 -1.32s",
            format!("{:.2}", rf2)
        );
    }

    #[test]
    fn normalization() {
        let rf = Rf::new([1., 2.], [-4., 6., -2.]);
        let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
        assert_eq!(expected, rf.normalize());

        let rf2 = Rf::new([1.], [0.]);
        assert_eq!(rf2, rf2.normalize());
    }

    #[test]
    fn rf_normalization_mutable() {
        let mut rf = Rf::new([1., 2.], [-4., 6., -2.]);
        rf.normalize_mut();
        let expected = Rf::new([-0.5, -1.], [2., -3., 1.]);
        assert_eq!(expected, rf);

        let mut rf2 = Rf::new([1.], [0.]);
        let rf3 = rf2.clone();
        rf2.normalize_mut();
        assert_eq!(rf2, rf3);
    }

    #[test]
    fn eval_rational_function() {
        let s_num = [-1., 1.];
        let s_den = [0., 1.];
        let s = Rf::new(s_num, s_den);
        let p = Poly::new_from_coeffs(&[1., 2., 3.]);
        let r = p.eval(&s);
        let expected = Rf::new([3., -8., 6.], [0., 0., 1.]);
        assert_eq!(expected, r);
    }
}