unthreaded/src/linops.rs

/*!
Abstract linear operators.
*/

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

/// Trait for linear operators on `X`.
pub trait Linear<X> {
    /// The range space of the operator.
    type Codomain;
    /// Apply the linear operator to `x`.
    fn apply(&self, x : &X) -> Self::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 Linear<X>>::Codomain> : Linear<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> {
    type FloatType : Float;
    /// 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::FloatType;
}

// 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;
    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;
    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
#[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
pub struct IdOp<X : Clone> (PhantomData<X>);

impl<X> IdOp<X> where X : Clone {
    fn new() -> IdOp<X> { IdOp(PhantomData) }
}

impl<X> Linear<X> for IdOp<X> where X : Clone {
    type Codomain = X;
    fn apply(&self, x : &X) -> X {
        x.clone()
    }
}

#[replace_float_literals(F::cast_from(literal))]
impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> where Y : AXPY<F, X>, X : Clone {
    // 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> BoundedLinear<X> for IdOp<X> where X : Clone {
    type FloatType = float;
    fn opnorm_bound(&self) -> float { 1.0 }
}

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