rmatrix_ks 0.4.0

//! # Vector
//!
//! some vector manipulations
//!
//! ColumnVector<T, row> = Matrix<T, row, 1>
//!
//! RowVector<T, col> = Matrix<T, 1, col>

use crate::error::IError;
use crate::error::IResult;
use crate::matrix::Matrix;
use crate::num::number::Fractional;
use crate::num::number::Number;
use crate::num::number::One;
use crate::num::number::Zero;
use crate::utils::state::SColumn;
use crate::utils::state::SRow;

/// type alias for row vector
pub type RowVector<T> = Matrix<T, SRow>;

/// type alias for column vector
pub type ColumnVector<T> = Matrix<T, SColumn>;

/// get the row en
///
/// for row vector, e1(3) = {{1, 0, 0}},
/// and e2(4) = {{0, 1, 0, 0}}
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::vector::identity_vector_row;
/// # use rmatrix_ks::vector::RowVector;
/// # fn main() -> IResult<()> {
/// assert_eq!(RowVector::create(1, 3, vec![1.0, 0.0, 0.0])?,
///     identity_vector_row::<f32>(3, 1)?);
/// assert_eq!(RowVector::create(1, 4, vec![0.0, 1.0, 0.0, 0.0])?,
///     identity_vector_row::<f32>(4, 2)?);
/// # Ok(())
/// # }
/// ```
pub fn identity_vector_row<T>(col: usize, index: usize) -> IResult<RowVector<T>>
where
    T: Clone + Zero + One,
{
    let mut vector = RowVector::zeros(1, col)?;
    vector.set_element(1, index, T::one())?;
    Ok(vector)
}

/// get the cloumn en
///
/// for column vector, e1(3) = {{1}, {0}, {0}},
/// and e2(4) = {{0}, {1}, {0}, {0}}
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::identity_vector_column;
/// # fn main() -> IResult<()> {
/// assert_eq!(ColumnVector::create(3, 1, vec![1.0, 0.0, 0.0])?,
///     identity_vector_column::<f32>(3, 1)?);
/// assert_eq!(ColumnVector::create(4, 1, vec![0.0, 1.0, 0.0, 0.0])?,
///     identity_vector_column::<f32>(4, 2)?);
/// # Ok(())
/// # }
/// ```
pub fn identity_vector_column<T>(row: usize, index: usize) -> IResult<ColumnVector<T>>
where
    T: Clone + Zero + One,
{
    let mut vector = ColumnVector::zeros(row, 1)?;
    vector.set_element(index, 1, T::one())?;
    Ok(vector)
}

/// vector cross product
///
/// only 3-d vector has cross product
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::times_c;
/// # fn main() -> IResult<()> {
/// let vec1: ColumnVector<f32> = ColumnVector::create(3, 1, vec![1.0f32, 2.0f32, 3.0f32])?;
/// let vec2: ColumnVector<f32> = ColumnVector::create(3, 1, vec![4.0f32, 5.0f32, 6.0f32])?;
/// assert_eq!(ColumnVector::create(3, 1, vec![-3.0f32, 6.0f32, -3.0f32])?, times_c(vec1, vec2)?);
/// # Ok(())
/// # }
/// ```
pub fn times_c<T>(vec1: ColumnVector<T>, vec2: ColumnVector<T>) -> IResult<ColumnVector<T>>
where
    T: Number,
{
    if vec1.row() != 3 {
        Err(IError::IncompatibleShape((3, 1), (vec1.row(), 1)))
    } else if vec2.row() != 3 {
        Err(IError::IncompatibleShape((3, 1), (vec2.row(), 1)))
    } else {
        ColumnVector::create(
            3,
            1,
            vec![
                vec1.get_element(2, 1)?.clone() * vec2.get_element(3, 1)?.clone()
                    - vec1.get_element(3, 1)?.clone() * vec2.get_element(2, 1)?.clone(),
                vec1.get_element(3, 1)?.clone() * vec2.get_element(1, 1)?.clone()
                    - vec1.get_element(1, 1)?.clone() * vec2.get_element(3, 1)?.clone(),
                vec1.get_element(1, 1)?.clone() * vec2.get_element(2, 1)?.clone()
                    - vec1.get_element(2, 1)?.clone() * vec2.get_element(1, 1)?.clone(),
            ],
        )
    }
}

/// vector dot product
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::RowVector;
/// # use rmatrix_ks::vector::times_d;
/// # fn main() -> IResult<()> {
/// let vec1: RowVector<f32> = RowVector::create(1, 3, vec![1.0f32, 2.0f32, 3.0f32])?;
/// let vec2: ColumnVector<f32> = ColumnVector::create(3, 1, vec![4.0f32, 5.0f32, 6.0f32])?;
/// assert_eq!(32.0f32, times_d(vec1, vec2)?);
/// # Ok(())
/// # }
/// ```
pub fn times_d<T>(vec1: RowVector<T>, vec2: ColumnVector<T>) -> IResult<T>
where
    T: Number,
{
    if vec1.column() != vec2.row() {
        Err(IError::IncompatibleShape(
            (vec1.column(), 1),
            (vec2.row(), 1),
        ))
    } else {
        Ok(vec1
            .get_inner()
            .zip(vec2.get_inner())
            .map(|(e1, e2)| e1.clone() * e2.clone())
            .sum())
    }
}

/// vector convolution
///
/// ```rust
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::convolution;
/// # use rmatrix_ks::error::IResult;
/// # fn main() -> IResult<()> {
/// let vec1: ColumnVector<f32> = ColumnVector::create(3, 1, vec![1.0f32, 2.0f32, 3.0f32])?;
/// let vec2: ColumnVector<f32> = ColumnVector::create(3, 1, vec![4.0f32, 5.0f32, 6.0f32])?;
/// assert_eq!(ColumnVector::create(5, 1, vec![4.0f32, 13.0f32, 28.0f32, 27.0f32, 18.0f32])?,
///     convolution(vec1, vec2)?);
/// # Ok(())
/// # }
/// ```
pub fn convolution<T>(vec1: ColumnVector<T>, vec2: ColumnVector<T>) -> IResult<ColumnVector<T>>
where
    T: Number,
{
    let edge: usize = vec1.row() + vec2.row() - 1;
    let mut conv: ColumnVector<T> = ColumnVector::zeros(edge, 1)?;

    for k in 1..=edge {
        for i in 1..=(k.min(vec1.row())) {
            // k + 1 = i + j
            let j = k + 1 - i;
            if (1..=vec2.row()).contains(&j) {
                conv.set_element(
                    k,
                    1,
                    conv.get_element(k, 1)?.clone()
                        + vec1.get_element(i, 1)?.clone() * vec2.get_element(j, 1)?.clone(),
                )?;
            }
        }
    }

    Ok(conv)
}

/// vector product
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::matrix::Matrix;
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::RowVector;
/// # use rmatrix_ks::vector::times_v;
/// # fn main() -> IResult<()> {
/// let vec1: ColumnVector<f32> = ColumnVector::create(3, 1, vec![1.0f32, 2.0f32, 3.0f32])?;
/// let vec2: RowVector<f32> = RowVector::create(1, 3, vec![4.0f32, 5.0f32, 6.0f32])?;
/// assert_eq!(Matrix::create(3, 3,
///     vec![
///         4.0f32, 5.0f32, 6.0f32,
///         8.0f32, 10.0f32, 12.0f32,
///         12.0f32, 15.0f32, 18.0f32])?,
///     times_v(vec1, vec2)?);
/// # Ok(())
/// # }
/// ```
pub fn times_v<T>(vector_c: ColumnVector<T>, vector_r: RowVector<T>) -> IResult<Matrix<T>>
where
    T: Number,
{
    let mut mat = Matrix::<T>::create(
        vector_c.row(),
        vector_r.column(),
        vec![T::default(); vector_c.row() * vector_r.column()],
    )?;
    for r in 1..=vector_c.row() {
        for c in 1..=vector_r.column() {
            mat.set_element(
                r,
                c,
                vector_c.get_element(r, 1)?.clone() * vector_r.get_element(1, c)?.clone(),
            )?;
        }
    }
    Ok(mat)
}

/// euclid norm for vector
///
/// aka l2-norm
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::euclid_norm;
/// # fn main() -> IResult<()> {
/// let v: ColumnVector<f32> = ColumnVector::create(2, 1, vec![3.0f32, 4.0f32])?;
/// assert_eq!(5.0f32, euclid_norm(v)?);
/// # Ok(())
/// # }
/// ```
pub fn euclid_norm<T>(vc: ColumnVector<T>) -> IResult<T>
where
    T: Fractional,
{
    times_d(vc.conjugate_transpose()?, vc)?.nsqrt()
}

/// root mean square for a vector
///
/// ```rust
/// # use rmatrix_ks::error::IResult;
/// # use rmatrix_ks::vector::ColumnVector;
/// # use rmatrix_ks::vector::root_mean_square;
/// # fn main() -> IResult<()> {
/// let v: ColumnVector<f32> = ColumnVector::create(2, 1, vec![3.0f32, 4.0f32])?;
/// assert_eq!(5.0f32 / 2.0f32.sqrt(), root_mean_square(v)?);
/// # Ok(())
/// # }
/// ```
pub fn root_mean_square<T>(vc: ColumnVector<T>) -> IResult<T>
where
    T: Fractional,
{
    let row = vc.row();
    euclid_norm(vc)?.ndiv(T::from_usize(row).nsqrt()?)
}