src/linops.rs

/*!
Abstract linear operators.
*/

use numeric_literals::replace_float_literals;
use std::marker::PhantomData;
use crate::types::*;
use serde::Serialize;
pub use crate::mapping::Apply;

/// Trait for linear operators on `X`.
pub trait Linear<X> : Apply<X, Output=Self::Codomain>
                      + for<'a> Apply<&'a X, Output=Self::Codomain>
                      + HasScalarField {
    type Codomain;
}

/// Efficient in-place summation.
#[replace_float_literals(F::cast_from(literal))]
pub trait AXPY<F : Num, X = Self> {
    /// Computes  `y = βy + αx`, where `y` is `Self`.
    fn axpy(&mut self, α : F, x : X, β : F);

    /// Copies `x` to `self`.
    fn copy_from(&mut self, x : X) {
        self.axpy(1.0, x, 0.0)
    }

    /// Computes  `y = αx`, where `y` is `Self`.
    fn scale_from(&mut self, α : F, x : X) {
        self.axpy(α, x, 0.0)
    }
}

/// Efficient in-place application for [`Linear`] operators.
#[replace_float_literals(F::cast_from(literal))]
pub trait GEMV<F : Num, X, Y = <Self as Apply<X>>::Output> : Apply<X> {
    /// Computes  `y = αAx + βy`, where `A` is `Self`.
    fn gemv(&self, y : &mut Y, α : F, x : X, β : F);

    /// Computes `y = Ax`, where `A` is `Self`
    fn apply_mut(&self, y : &mut Y, x : X){
        self.gemv(y, 1.0, x, 0.0)
    }

    /// Computes `y += Ax`, where `A` is `Self`
    fn apply_add(&self, y : &mut Y, x : X){
        self.gemv(y, 1.0, x, 1.0)
    }
}


/// Bounded linear operators
pub trait BoundedLinear<X> : Linear<X> {
    /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
    /// This is not expected to be the norm, just any bound on it that can be
    /// reasonably implemented.
    fn opnorm_bound(&self) -> Self::Field;
}

// Linear operator application into mutable target. The [`AsRef`] bound
// is used to guarantee compatibility with `Yʹ` and `Self::Codomain`;
// the former is assumed to be e.g. a view into the latter.

/*impl<X,Y,T> Fn(&X) -> Y for T where T : Linear<X,Codomain=Y> {
    fn call(&self, x : &X) -> Y {
        self.apply(x)
    }
}*/

/// Trait for forming the adjoint operator of `Self`.
pub trait Adjointable<X,Yʹ> : Linear<X> {
    type AdjointCodomain : HasScalarField<Field=Self::Field>;
    type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a;

    /// Form the adjoint operator of `self`.
    fn adjoint(&self) -> Self::Adjoint<'_>;

    /*fn adjoint_apply(&self, y : &Yʹ) -> Self::AdjointCodomain {
        self.adjoint().apply(y)
    }*/
}

/// Trait for forming a preadjoint of an operator.
///
/// For an operator $A$ this is an operator $A\_\*$
/// such that its adjoint $(A\_\*)^\*=A$. The space `X` is the domain of the `Self`
/// operator. The space `Ypre` is the predual of its codomain, and should be the
/// domain of the adjointed operator. `Self::Preadjoint` should be
/// [`Adjointable`]`<'a,Ypre,X>`.
pub trait Preadjointable<X,Ypre> : Linear<X> {
    type PreadjointCodomain : HasScalarField<Field=Self::Field>;
    type Preadjoint<'a> : Adjointable<Ypre, X, Codomain=Self::PreadjointCodomain> where Self : 'a;

    /// Form the preadjoint operator of `self`.
    fn preadjoint(&self) -> Self::Preadjoint<'_>;
}

/// Adjointable operators $A: X → Y$ on between reflexive spaces $X$ and $Y$.
pub trait SimplyAdjointable<X> : Adjointable<X,<Self as Linear<X>>::Codomain> {}
impl<'a,X,T> SimplyAdjointable<X> for T where T : Adjointable<X,<Self as Linear<X>>::Codomain> {}

/// The identity operator on `X` with scalar field `F`.
#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
pub struct IdOp<X, F : Num> (PhantomData<(X, F)>);

impl<X, F : Num> HasScalarField for IdOp<X, F> {
    type Field = F;
}

impl<X, F : Num> IdOp<X, F> {
    fn new() -> IdOp<X, F> { IdOp(PhantomData) }
}

impl<X, F : Num, T : CloneIfNeeded<X>> Apply<T> for IdOp<X, F> {
    type Output = X;

    fn apply(&self, x : T) -> X {
        x.clone_if_needed()
    }
}

impl<X : Clone, F : Num> Linear<X> for IdOp<X, F> {
    type Codomain = X;
}

#[replace_float_literals(F::cast_from(literal))]
impl<F : Num, X : Clone, Y> GEMV<F, X, Y> for IdOp<X, F> where Y : AXPY<F, X> {
    // Computes  `y = αAx + βy`, where `A` is `Self`.
    fn gemv(&self, y : &mut Y, α : F, x : X, β : F) {
        y.axpy(α, x, β)
    }

    fn apply_mut(&self, y : &mut Y, x : X){
        y.copy_from(x);
    }
}

impl<X, F : Num> BoundedLinear<X> for IdOp<X, F> where X : Clone {
    fn opnorm_bound(&self) -> F { F::ONE }
}

impl<X, F : Num> Adjointable<X,X> for IdOp<X, F> where X : Clone + HasScalarField<Field=F> {
    type AdjointCodomain=X;
    type Adjoint<'a> = IdOp<X, F> where X : 'a;
    fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
}