automatica 1.0.0

//! # Transfer function and matrices of transfer functions
//!
//! This module contains the generic methods for transfer functions
//! * calculation of zeros and poles (real and complex)
//! * arithmetic operations (addition, subtraction, multiplication, division,
//!   negation, inversion)
//! * positive and negative feedback
//! * conversion from a generic state-space representation of a single input
//!   single output system
//! * evaluation of the transfer function at the given complex number
//!
//! [continuous](continuous/index.html) module contains the specialized
//! structs and methods for continuous systems.
//!
//! [discrete](discrete/index.html) module contains the specialized structs and
//! methods for discrete systems.
//!

use std::{
    fmt,
    fmt::{Debug, Display, Formatter},
    marker::PhantomData,
    ops::{Add, AddAssign, Div, Mul, Neg, Sub},
};

use polynomen::Poly;

use crate::{
    enums::Time,
    error::{Error, ErrorKind},
    linear_system::{self, SsGen},
    polynomial_matrix::{MatrixOfPoly, PolyMatrix},
    rational_function::Rf,
    wrappers::{PWrapper, PP},
    Complex, Inv, NumCast, One, Polynomial, Zero,
};

pub mod continuous;
pub mod discrete;
pub mod discretization;
pub mod matrix;

/// Transfer function representation of a linear system
#[derive(Clone, Debug, PartialEq)]
pub struct TfGen<T, U: Time> {
    /// Rational function
    rf: Rf<T>,
    /// Tag to disambiguate continuous and discrete
    time: PhantomData<U>,
}

impl<T, U> TfGen<T, U>
where
    T: Clone + PartialEq + Zero,
    U: Time,
{
    /// Create a new transfer function given its numerator and denominator
    ///
    /// # Arguments
    ///
    /// * `num` - Transfer function numerator
    /// * `den` - Transfer function denominator
    ///
    /// # Example
    /// ```
    /// use automatica::Tfz;
    /// let tfz = Tfz::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 {
            rf: Rf::new(num, den),
            time: PhantomData::<U>,
        }
    }

    pub(crate) fn new_from_poly(num: Poly<PWrapper<T>>, den: Poly<PWrapper<T>>) -> Self {
        Self {
            rf: Rf::new_from_poly(num, den),
            time: PhantomData::<U>,
        }
    }
}

impl<T, U> TfGen<T, U>
where
    T: Clone + PartialEq + Zero,
    U: Time,
{
    /// Calculate the relative degree between denominator and numerator.
    ///
    /// # Example
    /// ```
    /// use automatica::{Inv, Tfz};
    /// let tfz = Tfz::new([1., 2.], [-4., 6., -2.]);
    /// let expected = tfz.relative_degree();
    /// assert_eq!(expected, 1);
    /// assert_eq!(tfz.inv().relative_degree(), -1);
    /// ```
    #[must_use]
    pub fn relative_degree(&self) -> i32 {
        self.rf.relative_degree()
    }
}

impl<T, U> TfGen<T, U>
where
    T: Clone + PartialEq + Zero,
    U: Time,
{
    /// Extract transfer function numerator
    ///
    /// # Example
    /// ```
    /// use automatica::{Polynomial, Tfz};
    /// let num = [1., 2.];
    /// let tfz = Tfz::new(num, [-4., 6., -2.]);
    /// assert_eq!(num, tfz.num().as_slice());
    /// ```
    #[must_use]
    pub fn num(&self) -> impl Polynomial<T> {
        self.rf.num()
    }

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

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

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

impl<T, U> TfGen<T, U>
where
    T: Clone + PartialEq,
    U: Time,
{
    /// Compute the reciprocal of a transfer function in place.
    pub fn inv_mut(&mut self) {
        self.rf.inv_mut();
    }
}

impl<T, U> Inv for &TfGen<T, U>
where
    T: Clone + PartialEq + Zero,
    U: Time,
{
    type Output = TfGen<T, U>;

    /// Compute the reciprocal of a transfer function.
    fn inv(self) -> Self::Output {
        Self::Output {
            rf: Inv::inv(&self.rf),
            time: PhantomData,
        }
    }
}

impl<T, U> Inv for TfGen<T, U>
where
    T: Clone + PartialEq,
    U: Time,
{
    type Output = Self;

    /// Compute the reciprocal of a transfer function.
    fn inv(mut self) -> Self::Output {
        self.inv_mut();
        self
    }
}

macro_rules! tf_roots {
    ($ty:ty) => {
        impl<U: Time> TfGen<$ty, U> {
            /// Calculate the poles of the transfer function
            #[must_use]
            pub fn real_poles(&self) -> Option<Vec<$ty>> {
                self.rf.real_poles()
            }

            /// Calculate the poles of the transfer function
            #[must_use]
            pub fn complex_poles(&self) -> Vec<Complex<$ty>> {
                self.rf
                    .complex_poles()
                    .into_iter()
                    .map(|(r, i)| Complex::new(r, i))
                    .collect()
            }

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

            /// Calculate the zeros of the transfer function
            #[must_use]
            pub fn complex_zeros(&self) -> Vec<Complex<$ty>> {
                self.rf
                    .complex_zeros()
                    .into_iter()
                    .map(|(r, i)| Complex::new(r, i))
                    .collect()
            }
        }
    };
}
tf_roots!(f32);
tf_roots!(f64);

impl<T, U> TfGen<T, U>
where
    T: Add<Output = T> + Clone + Div<Output = T> + One + PartialEq + Sub<Output = T> + Zero,
    U: Time,
{
    /// Negative feedback.
    ///
    /// ```text
    ///           L(s)
    /// G(s) = ----------
    ///         1 + L(s)
    /// ```
    /// where `self = L(s)`
    #[must_use]
    pub fn feedback_n(&self) -> Self {
        Self {
            rf: Rf::new_pp(
                PP(self.rf.num_int().0.clone()),
                PP(self.den_int().0.clone() + self.num_int().0.clone()),
            ),
            time: PhantomData,
        }
    }

    /// Positive feedback
    ///
    /// ```text
    ///           L(s)
    /// G(s) = ----------
    ///         1 - L(s)
    /// ```
    /// where `self = L(s)`
    #[must_use]
    pub fn feedback_p(&self) -> Self {
        Self {
            rf: Rf::new_pp(
                PP(self.rf.num_int().0.clone()),
                PP(self.den_int().0.clone() - self.num_int().0.clone()),
            ),
            time: PhantomData,
        }
    }

    /// Normalization of transfer function. If the denominator is zero the same
    /// transfer 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::Tfz;
    /// let tfz = Tfz::new([1., 2.], [-4., 6., -2.]);
    /// let expected = Tfz::new([-0.5, -1.], [2., -3., 1.]);
    /// assert_eq!(expected, tfz.normalize());
    /// ```
    #[must_use]
    pub fn normalize(&self) -> Self {
        Self {
            rf: self.rf.normalize(),
            time: PhantomData,
        }
    }

    /// In place normalization of transfer 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::Tfz;
    /// let mut tfz = Tfz::new([1., 2.], [-4., 6., -2.]);
    /// tfz.normalize_mut();
    /// let expected = Tfz::new([-0.5, -1.], [2., -3., 1.]);
    /// assert_eq!(expected, tfz);
    /// ```
    pub fn normalize_mut(&mut self) {
        self.rf.normalize_mut();
    }
}

impl<T, U: Time> TfGen<T, U>
where
    T: Add<Output = T>
        + AddAssign
        + Clone
        + Inv<Output = T>
        + Mul<Output = T>
        + Neg<Output = T>
        + NumCast
        + One
        + PartialEq
        + Zero,
{
    /// Convert a state-space representation into transfer functions.
    /// Conversion is available for Single Input Single Output system.
    /// If fails if the system is not SISO
    ///
    /// # Arguments
    ///
    /// `ss` - state space linear system
    ///
    /// # Errors
    ///
    /// It returns an error if the linear system is not single input
    /// single output.
    #[allow(clippy::option_if_let_else)]
    pub fn new_from_siso(ss: &SsGen<T, U>) -> Result<Self, Error> {
        let (pc, a_inv) = linear_system::leverrier(&ss.a)?;
        let g = a_inv.left_mul(&ss.c).right_mul(&ss.b);
        let rest = PolyMatrix::multiply(&pc, &ss.d);
        let tf = g + rest;
        if let Some(num) = MatrixOfPoly::from(tf).single() {
            Ok(Self::new_from_poly(num.clone().0, pc.0))
        } else {
            Err(Error::new_internal(ErrorKind::NoSisoSystem))
        }
    }
}

/// Implementation of transfer function negation.
/// Negative sign is transferred to the numerator.
impl<T, U> Neg for &TfGen<T, U>
where
    T: Clone + Neg<Output = T> + PartialEq + Zero,
    U: Time,
{
    type Output = TfGen<T, U>;

    fn neg(self) -> Self::Output {
        Self::Output {
            rf: -&self.rf,
            time: PhantomData,
        }
    }
}

/// Implementation of transfer function negation.
/// Negative sign is transferred to the numerator.
impl<T, U> Neg for TfGen<T, U>
where
    T: Clone + Neg<Output = T> + PartialEq,
    U: Time,
{
    type Output = Self;

    fn neg(mut self) -> Self::Output {
        self.rf = Neg::neg(self.rf);
        self
    }
}

/// Implementation of transfer function addition
#[allow(clippy::suspicious_arithmetic_impl)]
impl<T, U> Add for &TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = TfGen<T, U>;

    fn add(self, rhs: Self) -> Self::Output {
        Self::Output {
            rf: &self.rf + &rhs.rf,
            time: PhantomData,
        }
    }
}

/// Implementation of transfer function addition
impl<T, U> Add for TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = Self;

    fn add(mut self, rhs: Self) -> Self {
        self.rf = self.rf + rhs.rf;
        self
    }
}

/// Implementation of transfer function addition
impl<T, U> Add<T> for TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + PartialEq + Zero,
    U: Time,
{
    type Output = Self;

    fn add(mut self, rhs: T) -> Self {
        self.rf = self.rf + rhs;
        self
    }
}

/// Implementation of transfer function addition
impl<T, U> Add<&T> for TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + PartialEq + Zero,
    U: Time,
{
    type Output = Self;

    fn add(mut self, rhs: &T) -> Self {
        self.rf = self.rf + rhs;
        self
    }
}

/// Implementation of transfer function subtraction
impl<T, U> Sub for &TfGen<T, U>
where
    T: Add<Output = T>
        + Clone
        + Mul<Output = T>
        + Neg<Output = T>
        + One
        + PartialEq
        + Sub<Output = T>
        + Zero,
    U: Time,
{
    type Output = TfGen<T, U>;

    fn sub(self, rhs: Self) -> Self::Output {
        Self::Output {
            rf: Sub::sub(&self.rf, &rhs.rf),
            time: PhantomData,
        }
    }
}

/// Implementation of transfer function subtraction
impl<T, U> Sub for TfGen<T, U>
where
    T: Add<Output = T>
        + Clone
        + Mul<Output = T>
        + Neg<Output = T>
        + One
        + PartialEq
        + Sub<Output = T>
        + Zero,
    U: Time,
{
    type Output = Self;

    fn sub(mut self, rhs: Self) -> Self {
        self.rf = Sub::sub(self.rf, rhs.rf);
        self
    }
}

/// Implementation of transfer function multiplication
impl<T, U> Mul for &TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = TfGen<T, U>;

    fn mul(self, rhs: Self) -> Self::Output {
        Self::Output {
            rf: &self.rf * &rhs.rf,
            time: PhantomData,
        }
    }
}

/// Implementation of transfer function multiplication
impl<T, U> Mul for TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = Self;

    fn mul(mut self, rhs: Self) -> Self {
        self.rf = self.rf * rhs.rf;
        self
    }
}

/// Implementation of transfer function multiplication
impl<T, U> Mul<&TfGen<T, U>> for TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = Self;

    fn mul(mut self, rhs: &TfGen<T, U>) -> Self {
        self.rf = self.rf * &rhs.rf;
        self
    }
}

/// Implementation of transfer function division
impl<T, U> Div for &TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = TfGen<T, U>;

    fn div(self, rhs: Self) -> Self::Output {
        Self::Output {
            rf: Div::div(&self.rf, &rhs.rf),
            time: PhantomData,
        }
    }
}

/// Implementation of transfer function division
impl<T, U> Div for TfGen<T, U>
where
    T: Add<Output = T> + Clone + Mul<Output = T> + One + PartialEq + Zero,
    U: Time,
{
    type Output = Self;

    fn div(mut self, rhs: Self) -> Self {
        self.rf = Div::div(self.rf, rhs.rf);
        self
    }
}

impl<T, U> Zero for TfGen<T, U>
where
    T: Clone + One + PartialEq + Zero,
    U: Time,
{
    fn zero() -> Self {
        Self {
            rf: Rf::zero(),
            time: PhantomData,
        }
    }

    fn is_zero(&self) -> bool {
        self.rf.is_zero()
    }
}

impl<T, U> polynomen::Zero for TfGen<T, U>
where
    T: Clone + One + PartialEq + Zero,
    U: Time,
{
    fn zero() -> Self {
        Self {
            rf: Zero::zero(),
            time: PhantomData,
        }
    }

    fn is_zero(&self) -> bool {
        Zero::is_zero(&self.rf)
    }
}
impl<T, U> TfGen<T, U>
where
    T: Clone + PartialEq,
    U: Time,
{
    /// Evaluate the transfer function.
    ///
    /// # Arguments
    ///
    /// * `s` - Value at which the transfer function is evaluated.
    ///
    /// # Example
    /// ```
    /// use automatica::{Complex as C, Tf};
    /// let tf = Tf::new([1., 2., 3.], [-4., -3., 1.]);
    /// assert_eq!(-8.5, tf.eval_by_val(3.));
    /// assert_eq!(C::new(0.64, -0.98), tf.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,
    {
        self.rf.eval_by_val(s)
    }
}

impl<T, U> TfGen<T, U>
where
    T: Clone + PartialEq,
    U: Time,
{
    /// Evaluate the transfer function.
    ///
    /// # Arguments
    ///
    /// * `s` - Value at which the transfer function is evaluated.
    ///
    /// # Example
    /// ```
    /// use automatica::{Complex as C, Tf};
    /// let tf = Tf::new([1., 2., 3.], [-4., -3., 1.]);
    /// assert_eq!(-8.5, tf.eval(&3.));
    /// assert_eq!(C::new(0.64, -0.98), tf.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,
    {
        self.rf.eval(s)
    }
}

/// Implementation of transfer function printing
impl<T, U> Display for TfGen<T, U>
where
    T: Clone + PartialEq + Display + One + PartialEq + PartialOrd + Zero,
    U: Time,
{
    fn fmt(&self, f: &mut Formatter) -> fmt::Result {
        Display::fmt(&self.rf, f)
    }
}

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

    use crate::{Complex, Continuous, Discrete};

    use super::*;

    #[test]
    #[allow(clippy::float_cmp)]
    fn transfer_function_creation() {
        let num = [1.0_f64, 2., 3.];
        let den = [-4.2, -3.12, 0.0012];
        let tf = TfGen::<f64, Continuous>::new(num, den);
        assert_eq!(num, tf.num().as_slice());
        assert_eq!(den, tf.den().as_slice());
    }

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

    #[test]
    fn evaluation() {
        let tf = TfGen::<_, Continuous>::new([-0.75, 0.25], [0.75, 0.75, 1.]);
        let res = tf.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 tf = TfGen::<_, Continuous>::new([-0.75, 0.25], [0.75, 0.75, 1.]);
        let res1 = tf.eval(&Complex::new(0., 0.9));
        let res2 = tf.eval_by_val(Complex::new(0., 0.9));
        assert_eq!(res1, res2);
    }

    #[test]
    fn tf_inversion() {
        let num1 = [1., 2., 3.];
        let den1 = [-4.2, -3.12, 0.0012];
        let tf1 = TfGen::<_, Discrete>::new(num1, den1);
        let num2 = [-4.2, -3.12, 0.0012];
        let den2 = [1., 2., 3.];
        let mut tf2 = TfGen::new(num2, den2);

        assert_eq!(tf2, (&tf1).inv());
        tf2.inv_mut();
        assert_eq!(tf2, tf1);

        assert_eq!(tf2.inv(), tf1.inv());
    }

    #[test]
    fn poles() {
        let tf = TfGen::<f64, Continuous>::new([1.], [6., -5., 1.]);
        assert_eq!(Some(vec![2., 3.]), tf.real_poles());
    }

    #[test]
    fn complex_poles() {
        let tf = TfGen::<f32, Continuous>::new([1.], [10., -6., 1.]);
        assert_eq!(
            vec![Complex::new(3., -1.), Complex::new(3., 1.)],
            tf.complex_poles()
        );
    }

    #[test]
    fn zeros() {
        let tf = TfGen::<_, Discrete>::new([1.0_f64], [6., -5., 1.]);
        assert_eq!(None, tf.real_zeros());
    }

    #[test]
    fn complex_zeros() {
        let tf = TfGen::<f32, Discrete>::new([3.25, 3., 1.], [10., -3., 1.]);
        assert_eq!(
            vec![Complex::new(-1.5, -1.), Complex::new(-1.5, 1.)],
            tf.complex_zeros()
        );
    }

    proptest! {
        #[test]
        fn qc_tf_negative_feedback(b: f64) {
            let l = TfGen::<_, Continuous>::new([1.], [-b, 1.]);
            let g = TfGen::<_, Continuous>::new([1.], [-b + 1., 1.]);
            assert_eq!(g, l.feedback_n());
        }
    }

    proptest! {
    #[test]
        fn qc_tf_positive_feedback(b: f64) {
            let l = TfGen::<_, Continuous>::new([1.], [-b, 1.]);
            let g = TfGen::<_, Continuous>::new([1.], [-b - 1., 1.]);
            assert_eq!(g, l.feedback_p());
        }
    }

    #[test]
    fn tf_neg() {
        let tf1 = TfGen::<_, Discrete>::new([1., 2.], [1., 5.]);
        let tf2 = TfGen::<_, Discrete>::new([-1., -2.], [1., 5.]);
        assert_eq!(-&tf1, tf2);
        assert_eq!(tf1, -(-(&tf1)));
    }

    #[test]
    fn add_references() {
        let tf1 = TfGen::<_, Continuous>::new([1., 2.], [1., 5.]);
        let tf2 = TfGen::new([3.], [1., 5.]);
        let actual = &tf1 + &tf2;
        let expected = TfGen::new([4., 2.], [1., 5.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, &expected + &TfGen::zero());
        assert_eq!(expected, &TfGen::zero() + &expected);
    }

    #[test]
    fn add_values() {
        let tf1 = TfGen::<_, Discrete>::new([1., 2.], [3., -4.]);
        let tf2 = TfGen::new([3.], [1., 5.]);
        let actual = tf1 + tf2;
        let expected = TfGen::new([10., -5., 10.], [3., 11., -20.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, expected.clone() + TfGen::zero());
        assert_eq!(expected, TfGen::<f64, Discrete>::zero() + expected.clone());
    }

    #[test]
    fn add_multiple_values() {
        let tf1 = TfGen::<_, Discrete>::new([1., 2.], [3., -4.]);
        let tf2 = TfGen::new([3.], [1., 5.]);
        let tf3 = TfGen::new([0., 4.], [3., 11., -20.]);
        let actual = &(&tf1 + &tf2) + &tf3;
        let expected = TfGen::new([10., -1., 10.], [3., 11., -20.]);
        assert_eq!(expected, actual);

        let actual2 = (tf1 + tf2) + tf3;
        assert_eq!(actual, actual2);
    }

    #[test]
    fn add_scalar_value() {
        let tf = TfGen::<_, Discrete>::new([1., 2.], [3., -4.]);
        let actual = tf + 1.;
        let expected = TfGen::new([4., -2.], [3., -4.]);
        assert_eq!(expected, actual);
    }

    #[test]
    #[allow(clippy::op_ref)]
    fn add_scalar_reference() {
        let tf = TfGen::<_, Discrete>::new([1., 2.], [3., -4.]);
        let actual = tf + &2.;
        let expected = TfGen::new([7., -6.], [3., -4.]);
        assert_eq!(expected, actual);
    }

    #[test]
    fn sub_references() {
        let tf1 = TfGen::<_, Continuous>::new([-1., 9.], [4., -1.]);
        let tf2 = TfGen::new([3.], [4., -1.]);
        let actual = &tf1 - &tf2;
        let expected = TfGen::new([-4., 9.], [4., -1.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, &expected - &TfGen::zero());
        assert_eq!(-&expected, &TfGen::zero() - &expected);
    }

    #[test]
    fn sub_values() {
        let tf1 = TfGen::<_, Discrete>::new([1., -2.], [4., -4.]);
        let tf2 = TfGen::new([2.], [2., 3.]);
        let actual = tf1 - tf2;
        let expected = TfGen::new([-6., 7., -6.], [8., 4., -12.]);
        assert_eq!(expected, actual);
        assert_eq!(expected, expected.clone() - TfGen::zero());
        assert_eq!(-&expected, TfGen::zero() - expected);
    }

    #[test]
    fn sub_multiple_values() {
        let tf1 = TfGen::<_, Discrete>::new([1., -2.], [4., -4.]);
        let tf2 = TfGen::new([2.], [2., 3.]);
        let tf3 = TfGen::new([0., 2.], [8., 4., -12.]);
        let actual = &(&tf1 - &tf2) - &tf3;
        let expected = TfGen::new([-6., 5., -6.], [8., 4., -12.]);
        assert_eq!(expected, actual);

        let actual2 = (tf1 - tf2) - tf3;
        assert_eq!(actual, actual2);
    }

    #[test]
    fn mul_references() {
        let tf1 = TfGen::<_, Continuous>::new([1., 2., 3.], [1., 5.]);
        let tf2 = TfGen::new([3.], [1., 6., 5.]);
        let actual = &tf1 * &tf2;
        let expected = TfGen::new([3., 6., 9.], [1., 11., 35., 25.]);
        assert_eq!(expected, actual);
        assert_eq!(TfGen::zero(), &expected * &TfGen::zero());
        assert_eq!(TfGen::zero(), &TfGen::zero() * &expected);
    }

    #[test]
    fn mul_values() {
        let tf1 = TfGen::<_, Continuous>::new([1., 2., 3.], [1., 5.]);
        let tf2 = TfGen::new([-5.], [1., 6., 5.]);
        let actual = tf1 * tf2;
        let expected = TfGen::new([-5., -10., -15.], [1., 11., 35., 25.]);
        assert_eq!(expected, actual);
        assert_eq!(TfGen::zero(), expected.clone() * TfGen::zero());
        assert_eq!(TfGen::zero(), TfGen::zero() * expected);
    }

    #[test]
    fn mul_value_reference() {
        let tf1 = TfGen::<_, Continuous>::new([1., 2., 3.], [1., 5.]);
        let tf2 = TfGen::new([-5.], [1., 6., 5.]);
        let actual = tf1 * &tf2;
        let expected = TfGen::new([-5., -10., -15.], [1., 11., 35., 25.]);
        assert_eq!(expected, actual);
        assert_eq!(TfGen::zero(), expected.clone() * &TfGen::zero());
        assert_eq!(TfGen::zero(), TfGen::zero() * &expected);
    }

    #[test]
    fn div_values() {
        let tf1 = TfGen::<f64, Discrete>::new([1., 2., 3.], [1., 5.]);
        let tf2 = TfGen::new([3.], [1., 6., 5.]);
        let actual = tf2 / tf1;
        let expected = TfGen::new([3., 15.], [1., 8., 20., 28., 15.]);
        assert_eq!(expected, actual);
        assert_eq!(TfGen::zero(), &TfGen::zero() / &expected);
        assert!((&expected / &TfGen::zero()).eval(&1.).is_infinite());
        assert!((&TfGen::<f32, Discrete>::zero() / &TfGen::zero())
            .eval(&1.)
            .is_nan());
    }

    #[test]
    #[allow(clippy::eq_op)]
    fn div_references() {
        let tf1 = TfGen::<f64, Discrete>::new([1., 2., 3.], [1., 5.]);
        let actual = &tf1 / &tf1;
        let expected = TfGen::new([1., 7., 13., 15.], [1., 7., 13., 15.]);
        assert_eq!(expected, actual);
        assert_eq!(TfGen::zero(), TfGen::zero() / expected.clone());
        assert!((expected / TfGen::zero()).eval(&1.).is_infinite());
        assert!((TfGen::<f32, Discrete>::zero() / TfGen::zero())
            .eval(&1.)
            .is_nan());
    }

    #[test]
    fn zero_tf() {
        assert!(TfGen::<f32, Continuous>::zero().is_zero());
        assert!(!TfGen::<_, Discrete>::new([0.], [0.]).is_zero());
    }

    #[test]
    fn zero_polynomen_trait() {
        assert!(<TfGen::<f32, Continuous> as polynomen::Zero>::zero().is_zero());
    }

    #[test]
    fn print() {
        let tf = TfGen::<_, Continuous>::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!("{}", tf)
        );

        let tf2 = TfGen::<_, Continuous>::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}", tf2)
        );
    }

    #[test]
    fn normalization() {
        let tfz = TfGen::<_, Discrete>::new([1., 2.], [-4., 6., -2.]);
        let expected = TfGen::new([-0.5, -1.], [2., -3., 1.]);
        assert_eq!(expected, tfz.normalize());

        let tfz2 = TfGen::<_, Discrete>::new([1.], [0.]);
        assert_eq!(tfz2, tfz2.normalize());
    }

    #[test]
    fn normalization_mutable() {
        let mut tfz = TfGen::<_, Discrete>::new([1., 2.], [-4., 6., -2.]);
        tfz.normalize_mut();
        let expected = TfGen::new([-0.5, -1.], [2., -3., 1.]);
        assert_eq!(expected, tfz);

        let mut tfz2 = TfGen::<_, Discrete>::new([1.], [0.]);
        let tfz3 = tfz2.clone();
        tfz2.normalize_mut();
        assert_eq!(tfz2, tfz3);
    }

    #[test]
    fn failed_conversion_from_ss() {
        let ss = crate::Ss::new_from_slice(
            2,
            2,
            1,
            &[1.0_f32, 1., 1., 1.],
            &[1., 1., 1., 1.],
            &[1., 1.],
            &[1., 1.],
        );
        let res = TfGen::<f32, Continuous>::new_from_siso(&ss);
        assert!(res.is_err());
    }

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