src/linops.rs

/*!
Abstract linear operators.
*/

use crate::direct_product::Pair;
use crate::error::DynResult;
use crate::instance::Instance;
pub use crate::mapping::{ClosedSpace, Composition, DifferentiableImpl, Mapping, Space};
use crate::norms::{HasDual, Linfinity, Norm, NormExponent, PairNorm, L1, L2};
use crate::types::*;
use numeric_literals::replace_float_literals;
use serde::Serialize;
use std::marker::PhantomData;
use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};

/// Trait for linear operators on `X`.
pub trait Linear<X: Space>: Mapping<X> {}

// impl<X: Space, A: Linear<X>> DifferentiableImpl<X> for A {
//     type Derivative = <Self as Mapping<X>>::Codomain;

//     /// Compute the differential of `self` at `x`, consuming the input.
//     fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
//         self.apply(x)
//     }
// }

/// Vector spaces
#[replace_float_literals(Self::Field::cast_from(literal))]
pub trait VectorSpace:
    Space<Principal = Self::PrincipalV>
    + Mul<Self::Field, Output = Self::PrincipalV>
    + Div<Self::Field, Output = Self::PrincipalV>
    + Add<Self, Output = Self::PrincipalV>
    + Add<Self::PrincipalV, Output = Self::PrincipalV>
    + Sub<Self, Output = Self::PrincipalV>
    + Sub<Self::PrincipalV, Output = Self::PrincipalV>
    + Neg
    + for<'b> Add<&'b Self, Output = <Self as VectorSpace>::PrincipalV>
    + for<'b> Sub<&'b Self, Output = <Self as VectorSpace>::PrincipalV>
{
    /// Underlying scalar field
    type Field: Num;

    /// Principal form of the space; always equal to [`Space::Principal`], but with
    /// more traits guaranteed.
    ///
    /// `PrincipalV` is only assumed to be `AXPY` for itself, as [`AXPY`]
    /// uses [`Instance`] to apply all other variants and avoid problems
    /// of choosing multiple implementations of the trait.
    type PrincipalV: ClosedSpace
        + AXPY<
            Self::PrincipalV,
            Field = Self::Field,
            PrincipalV = Self::PrincipalV,
            OwnedVariant = Self::PrincipalV,
            Principal = Self::PrincipalV,
        >;

    /// Return a similar zero as `self`.
    fn similar_origin(&self) -> Self::PrincipalV;
    // {
    //     self.make_origin_generator().make_origin()
    // }

    /// Return a similar zero as `x`.
    fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::PrincipalV {
        x.eval(|xr| xr.similar_origin())
    }
}

/// Efficient in-place summation.
#[replace_float_literals(Self::Field::cast_from(literal))]
pub trait AXPY<X = Self>:
    VectorSpace
    + MulAssign<Self::Field>
    + DivAssign<Self::Field>
    + AddAssign<Self>
    + AddAssign<Self::PrincipalV>
    + SubAssign<Self>
    + SubAssign<Self::PrincipalV>
    + for<'b> AddAssign<&'b Self>
    + for<'b> SubAssign<&'b Self>
where
    X: Space,
{
    /// Computes  `y = βy + αx`, where `y` is `Self`.
    fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field);

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

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

    /// Set self to zero.
    fn set_zero(&mut self);
}

pub trait ClosedVectorSpace: Instance<Self> + VectorSpace<PrincipalV = Self> {}
impl<X: Instance<X> + VectorSpace<PrincipalV = Self>> ClosedVectorSpace for X {}

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

    #[inline]
    /// Computes `y = Ax`, where `A` is `Self`
    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
        self.gemv(y, 1.0, x, 0.0)
    }

    #[inline]
    /// Computes `y += Ax`, where `A` is `Self`
    fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
        self.gemv(y, 1.0, x, 1.0)
    }
}

/// Bounded linear operators
pub trait BoundedLinear<X, XExp, CodExp, F = f64>: Linear<X>
where
    F: Num,
    X: Space + Norm<XExp, F>,
    XExp: NormExponent,
    CodExp: NormExponent,
{
    /// 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. The [`NormExponent`] `xexp`  indicates the norm
    /// in `X`, and `codexp` in the codomain.
    ///
    /// This may fail with an error if the bound is for some reason incalculable.
    fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> DynResult<F>;
}

// 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>
where
    X: Space,
    Yʹ: Space,
{
    /// Codomain of the adjoint operator.
    type AdjointCodomain: ClosedSpace;
    /// Type of the adjoint operator.
    type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain>
    where
        Self: 'a;

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

/// Variant of [`Adjointable`] where the adjoint does not depend on a lifetime parameter.
/// This exists due to restrictions of Rust's type system: if `A :: Adjointable`, and we make
/// further restrictions on the adjoint operator, through, e.g.
/// ```
/// for<'a> A::Adjoint<'a> : GEMV<F, X, Y>,
/// ```
/// Then `'static` lifetime is forced on `X`. Having `A::SimpleAdjoint` not depend on `'a`
/// avoids this, but makes it impossible for the adjoint to be just a light wrapper around the
/// original operator.
pub trait SimplyAdjointable<X, Yʹ>: Linear<X>
where
    X: Space,
    Yʹ: Space,
{
    /// Codomain of the adjoint operator.
    type AdjointCodomain: ClosedSpace;
    /// Type of the adjoint operator.
    type SimpleAdjoint: Linear<Yʹ, Codomain = Self::AdjointCodomain>;

    /// Form the adjoint operator of `self`.
    fn adjoint(&self) -> Self::SimpleAdjoint;
}

/// 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>`.
/// We do not make additional restrictions on `Self::Preadjoint` (in particular, it
/// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint
/// have the full space as the codomain, etc.
pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>:
    Linear<X>
{
    type PreadjointCodomain: ClosedSpace;
    type Preadjoint<'a>: Linear<Ypre, Codomain = Self::PreadjointCodomain>
    where
        Self: 'a;

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

/// The identity operator
#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct IdOp<X>(PhantomData<X>);

impl<X> IdOp<X> {
    pub fn new() -> IdOp<X> {
        IdOp(PhantomData)
    }
}

impl<X: Space> Mapping<X> for IdOp<X> {
    type Codomain = X::Principal;

    fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain {
        x.own()
    }
}

impl<X: Space> Linear<X> for IdOp<X> {}

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

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

impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
where
    X: Space + Clone + Norm<E, F>,
    F: Num,
    E: NormExponent,
{
    fn opnorm_bound(&self, _xexp: E, _codexp: E) -> DynResult<F> {
        Ok(F::ONE)
    }
}

impl<X: Clone + Space> Adjointable<X, X::Principal> for IdOp<X> {
    type AdjointCodomain = X::Principal;
    type Adjoint<'a>
        = IdOp<X::Principal>
    where
        X: 'a;

    fn adjoint(&self) -> Self::Adjoint<'_> {
        IdOp::new()
    }
}

impl<X: Clone + Space> SimplyAdjointable<X, X::Principal> for IdOp<X> {
    type AdjointCodomain = X::Principal;
    type SimpleAdjoint = IdOp<X::Principal>;

    fn adjoint(&self) -> Self::SimpleAdjoint {
        IdOp::new()
    }
}

impl<X: Clone + Space> Preadjointable<X, X::Principal> for IdOp<X> {
    type PreadjointCodomain = X::Principal;
    type Preadjoint<'a>
        = IdOp<X::Principal>
    where
        X: 'a;

    fn preadjoint(&self) -> Self::Preadjoint<'_> {
        IdOp::new()
    }
}

/// The zero operator from a space to itself
#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct SimpleZeroOp;

impl<X: VectorSpace> Mapping<X> for SimpleZeroOp {
    type Codomain = X::PrincipalV;

    fn apply<I: Instance<X>>(&self, x: I) -> X::PrincipalV {
        X::similar_origin_inst(x)
    }
}

impl<X: VectorSpace> Linear<X> for SimpleZeroOp {}

#[replace_float_literals(F::cast_from(literal))]
impl<X, Y, F> GEMV<F, X, Y> for SimpleZeroOp
where
    F: Num,
    Y: AXPY<Field = F>,
    X: VectorSpace<Field = F> + Instance<X>,
{
    // Computes  `y = αAx + βy`, where `A` is `Self`.
    fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
        *y *= β;
    }

    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
        y.set_zero();
    }
}

impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp
where
    F: Num,
    X: VectorSpace<Field = F> + Norm<E1, F>,
    E1: NormExponent,
    E2: NormExponent,
{
    fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
        Ok(F::ZERO)
    }
}

impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp
where
    F: Num,
    X: VectorSpace<Field = F> + HasDual<F>,
    X::DualSpace: ClosedVectorSpace,
{
    type AdjointCodomain = X::DualSpace;
    type Adjoint<'b>
        = SimpleZeroOp
    where
        Self: 'b;
    // () means not (pre)adjointable.

    fn adjoint(&self) -> Self::Adjoint<'_> {
        SimpleZeroOp
    }
}

pub trait OriginGenerator<Y: VectorSpace> {
    type Ref<'b>: OriginGenerator<Y>
    where
        Self: 'b;

    fn origin(&self) -> Y::PrincipalV;
    fn as_ref(&self) -> Self::Ref<'_>;
}

impl<Y: VectorSpace> OriginGenerator<Y> for Y {
    type Ref<'b>
        = &'b Y
    where
        Self: 'b;

    #[inline]
    fn origin(&self) -> Y::PrincipalV {
        return self.similar_origin();
    }

    #[inline]
    fn as_ref(&self) -> Self::Ref<'_> {
        self
    }
}

impl<'b, Y: VectorSpace> OriginGenerator<Y> for &'b Y {
    type Ref<'c>
        = Self
    where
        Self: 'c;

    #[inline]
    fn origin(&self) -> Y::PrincipalV {
        return self.similar_origin();
    }

    #[inline]
    fn as_ref(&self) -> Self::Ref<'_> {
        self
    }
}

/// A zero operator that can be eitherh dualised or predualised (once).
/// This is achieved by storing an oppropriate zero.
pub struct ZeroOp<X, Y: VectorSpace<Field = F>, OY: OriginGenerator<Y>, O, F: Float = f64> {
    codomain_origin_generator: OY,
    other_origin_generator: O,
    _phantoms: PhantomData<(X, Y, F)>,
}

impl<X, Y, OY, F> ZeroOp<X, Y, OY, (), F>
where
    OY: OriginGenerator<Y>,
    X: VectorSpace<Field = F>,
    Y: VectorSpace<Field = F>,
    F: Float,
{
    pub fn new(y_og: OY) -> Self {
        ZeroOp {
            codomain_origin_generator: y_og,
            other_origin_generator: (),
            _phantoms: PhantomData,
        }
    }
}

impl<X, Y, OY, OXprime, Xprime, Yprime, F> ZeroOp<X, Y, OY, OXprime, F>
where
    OY: OriginGenerator<Y>,
    OXprime: OriginGenerator<Xprime>,
    X: HasDual<F, DualSpace = Xprime>,
    Y: HasDual<F, DualSpace = Yprime>,
    F: Float,
    Xprime: VectorSpace<Field = F, PrincipalV = Xprime>,
    Xprime::PrincipalV: AXPY<Field = F>,
    Yprime: Space + Instance<Yprime>,
{
    pub fn new_dualisable(y_og: OY, xprime_og: OXprime) -> Self {
        ZeroOp {
            codomain_origin_generator: y_og,
            other_origin_generator: xprime_og,
            _phantoms: PhantomData,
        }
    }
}

impl<X, Y, O, OY, F> Mapping<X> for ZeroOp<X, Y, OY, O, F>
where
    X: Space,
    Y: VectorSpace<Field = F>,
    F: Float,
    OY: OriginGenerator<Y>,
{
    type Codomain = Y::PrincipalV;

    fn apply<I: Instance<X>>(&self, _x: I) -> Y::PrincipalV {
        self.codomain_origin_generator.origin()
    }
}

impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F>
where
    X: Space,
    Y: VectorSpace<Field = F>,
    F: Float,
    OY: OriginGenerator<Y>,
{
}

#[replace_float_literals(F::cast_from(literal))]
impl<X, Y, OY, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, OY, O, F>
where
    X: Space,
    Y: AXPY<Field = F, PrincipalV = Y>,
    F: Float,
    OY: OriginGenerator<Y>,
{
    // Computes  `y = αAx + βy`, where `A` is `Self`.
    fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
        *y *= β;
    }

    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
        y.set_zero();
    }
}

impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F>
where
    X: Space + Instance<X> + Norm<E1, F>,
    Y: VectorSpace<Field = F>,
    Y::PrincipalV: Clone,
    F: Float,
    E1: NormExponent,
    E2: NormExponent,
    OY: OriginGenerator<Y>,
{
    fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
        Ok(F::ZERO)
    }
}

impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> Adjointable<X, Yprime>
    for ZeroOp<X, Y, OY, OXprime, F>
where
    X: HasDual<F, DualSpace = Xprime>,
    Y: HasDual<F, DualSpace = Yprime>,
    F: Float,
    Xprime: ClosedVectorSpace<Field = F>,
    //Xprime::Owned: AXPY<Field = F>,
    Yprime: ClosedSpace,
    OY: OriginGenerator<Y>,
    OXprime: OriginGenerator<X::DualSpace>,
{
    type AdjointCodomain = Xprime;
    type Adjoint<'c>
        = ZeroOp<Yprime, Xprime, OXprime::Ref<'c>, (), F>
    where
        Self: 'c;
    // () means not (pre)adjointable.

    fn adjoint(&self) -> Self::Adjoint<'_> {
        ZeroOp {
            codomain_origin_generator: self.other_origin_generator.as_ref(),
            other_origin_generator: (),
            _phantoms: PhantomData,
        }
    }
}

impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> SimplyAdjointable<X, Yprime>
    for ZeroOp<X, Y, OY, OXprime, F>
where
    X: HasDual<F, DualSpace = Xprime>,
    Y: HasDual<F, DualSpace = Yprime>,
    F: Float,
    Xprime: ClosedVectorSpace<Field = F>,
    //Xprime::Owned: AXPY<Field = F>,
    Yprime: ClosedSpace,
    OY: OriginGenerator<Y>,
    OXprime: OriginGenerator<X::DualSpace> + Clone,
{
    type AdjointCodomain = Xprime;
    type SimpleAdjoint = ZeroOp<Yprime, Xprime, OXprime, (), F>;
    // () means not (pre)adjointable.

    fn adjoint(&self) -> Self::SimpleAdjoint {
        ZeroOp {
            codomain_origin_generator: self.other_origin_generator.clone(),
            other_origin_generator: (),
            _phantoms: PhantomData,
        }
    }
}

impl<S, T, E, X> Linear<X> for Composition<S, T, E>
where
    X: Space,
    T: Linear<X>,
    S: Linear<T::Codomain>,
{
}

impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
where
    F: Num,
    X: Space,
    T: Linear<X>,
    S: GEMV<F, T::Codomain, Y>,
{
    fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
        self.outer.gemv(y, α, self.inner.apply(x), β)
    }

    /// Computes `y = Ax`, where `A` is `Self`
    fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
        self.outer.apply_mut(y, self.inner.apply(x))
    }

    /// Computes `y += Ax`, where `A` is `Self`
    fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
        self.outer.apply_add(y, self.inner.apply(x))
    }
}

impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
where
    F: Num,
    X: Space + Norm<Xexp, F>,
    Z: Space + Norm<Zexp, F>,
    Xexp: NormExponent,
    Yexp: NormExponent,
    Zexp: NormExponent,
    T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>,
    S: BoundedLinear<Z, Zexp, Yexp, F>,
{
    fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> DynResult<F> {
        let zexp = self.intermediate_norm_exponent;
        Ok(self.outer.opnorm_bound(zexp, yexp)? * self.inner.opnorm_bound(xexp, zexp)?)
    }
}

/// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct RowOp<S, T>(pub S, pub T);

impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
where
    A: Space,
    B: Space,
    S: Mapping<A>,
    T: Mapping<B>,
    S::Codomain: Add<T::Codomain>,
    <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
{
    type Codomain = <S::Codomain as Add<T::Codomain>>::Output;

    fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
        x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b))
    }
}

impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
where
    A: Space,
    B: Space,
    S: Linear<A>,
    T: Linear<B>,
    S::Codomain: Add<T::Codomain>,
    <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
{
}

impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
where
    U: Space,
    V: Space,
    S: GEMV<F, U, Y>,
    T: GEMV<F, V, Y>,
    F: Num,
    Self: Linear<Pair<U, V>, Codomain = Y>,
{
    fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Y, α: F, x: I, β: F) {
        x.eval_decompose(|Pair(u, v)| {
            self.0.gemv(y, α, u, β);
            self.1.gemv(y, α, v, F::ONE);
        })
    }

    fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
        x.eval_decompose(|Pair(u, v)| {
            self.0.apply_mut(y, u);
            self.1.apply_add(y, v);
        })
    }

    /// Computes `y += Ax`, where `A` is `Self`
    fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
        x.eval_decompose(|Pair(u, v)| {
            self.0.apply_add(y, u);
            self.1.apply_add(y, v);
        })
    }
}

/// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct ColOp<S, T>(pub S, pub T);

impl<A, S, T> Mapping<A> for ColOp<S, T>
where
    A: Space,
    S: Mapping<A>,
    T: Mapping<A>,
{
    type Codomain = Pair<S::Codomain, T::Codomain>;

    fn apply<I: Instance<A>>(&self, a: I) -> Self::Codomain {
        Pair(a.eval_ref(|r| self.0.apply(r)), self.1.apply(a))
    }
}

impl<A, S, T> Linear<A> for ColOp<S, T>
where
    A: Space,
    S: Mapping<A>,
    T: Mapping<A>,
{
}

impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
where
    X: Space,
    S: GEMV<F, X, A>,
    T: GEMV<F, X, B>,
    F: Num,
    Self: Linear<X, Codomain = Pair<A, B>>,
{
    fn gemv<I: Instance<X>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
        x.eval_ref(|r| self.0.gemv(&mut y.0, α, r, β));
        self.1.gemv(&mut y.1, α, x, β);
    }

    fn apply_mut<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
        x.eval_ref(|r| self.0.apply_mut(&mut y.0, r));
        self.1.apply_mut(&mut y.1, x);
    }

    /// Computes `y += Ax`, where `A` is `Self`
    fn apply_add<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
        x.eval_ref(|r| self.0.apply_add(&mut y.0, r));
        self.1.apply_add(&mut y.1, x);
    }
}

impl<A, B, Yʹ, S, T> Adjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
where
    A: Space,
    B: Space,
    Yʹ: Space,
    S: Adjointable<A, Yʹ>,
    T: Adjointable<B, Yʹ>,
    Self: Linear<Pair<A, B>>,
    // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
    //     Yʹ,
    //     Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
    // >,
{
    type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
    type Adjoint<'a>
        = ColOp<S::Adjoint<'a>, T::Adjoint<'a>>
    where
        Self: 'a;

    fn adjoint(&self) -> Self::Adjoint<'_> {
        ColOp(self.0.adjoint(), self.1.adjoint())
    }
}

impl<A, B, Yʹ, S, T> SimplyAdjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
where
    A: Space,
    B: Space,
    Yʹ: Space,
    S: SimplyAdjointable<A, Yʹ>,
    T: SimplyAdjointable<B, Yʹ>,
    Self: Linear<Pair<A, B>>,
    // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
    //     Yʹ,
    //     Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
    // >,
{
    type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
    type SimpleAdjoint = ColOp<S::SimpleAdjoint, T::SimpleAdjoint>;

    fn adjoint(&self) -> Self::SimpleAdjoint {
        ColOp(self.0.adjoint(), self.1.adjoint())
    }
}

impl<A, B, Yʹ, S, T> Preadjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
where
    A: Space,
    B: Space,
    Yʹ: Space,
    S: Preadjointable<A, Yʹ>,
    T: Preadjointable<B, Yʹ>,
    Self: Linear<Pair<A, B>>,
    for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>:
        Linear<Yʹ, Codomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>>,
{
    type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
    type Preadjoint<'a>
        = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
    where
        Self: 'a;

    fn preadjoint(&self) -> Self::Preadjoint<'_> {
        ColOp(self.0.preadjoint(), self.1.preadjoint())
    }
}

impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
where
    A: Space,
    Xʹ: Space,
    Yʹ: Space,
    R: ClosedSpace + ClosedAdd,
    S: Adjointable<A, Xʹ, AdjointCodomain = R>,
    T: Adjointable<A, Yʹ, AdjointCodomain = R>,
    Self: Linear<A>,
    // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
    //     Pair<Xʹ,Yʹ>,
    //     Codomain=R,
    // >,
{
    type AdjointCodomain = R;
    type Adjoint<'a>
        = RowOp<S::Adjoint<'a>, T::Adjoint<'a>>
    where
        Self: 'a;

    fn adjoint(&self) -> Self::Adjoint<'_> {
        RowOp(self.0.adjoint(), self.1.adjoint())
    }
}

impl<A, Xʹ, Yʹ, R, S, T> SimplyAdjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
where
    A: Space,
    Xʹ: Space,
    Yʹ: Space,
    R: ClosedSpace + ClosedAdd,
    S: SimplyAdjointable<A, Xʹ, AdjointCodomain = R>,
    T: SimplyAdjointable<A, Yʹ, AdjointCodomain = R>,
    Self: Linear<A>,
    // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
    //     Pair<Xʹ,Yʹ>,
    //     Codomain=R,
    // >,
{
    type AdjointCodomain = R;
    type SimpleAdjoint = RowOp<S::SimpleAdjoint, T::SimpleAdjoint>;

    fn adjoint(&self) -> Self::SimpleAdjoint {
        RowOp(self.0.adjoint(), self.1.adjoint())
    }
}

impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
where
    A: Space,
    Xʹ: Space,
    Yʹ: Space,
    R: ClosedSpace + ClosedAdd,
    S: Preadjointable<A, Xʹ, PreadjointCodomain = R>,
    T: Preadjointable<A, Yʹ, PreadjointCodomain = R>,
    Self: Linear<A>,
    for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
{
    type PreadjointCodomain = R;
    type Preadjoint<'a>
        = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
    where
        Self: 'a;

    fn preadjoint(&self) -> Self::Preadjoint<'_> {
        RowOp(self.0.preadjoint(), self.1.preadjoint())
    }
}

/// Diagonal operator
#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
pub struct DiagOp<S, T>(pub S, pub T);

impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
where
    A: Space,
    B: Space,
    S: Mapping<A>,
    T: Mapping<B>,
{
    type Codomain = Pair<S::Codomain, T::Codomain>;

    fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
        x.eval_decompose(|Pair(a, b)| Pair(self.0.apply(a), self.1.apply(b)))
    }
}

impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
where
    A: Space,
    B: Space,
    S: Linear<A>,
    T: Linear<B>,
{
}

impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
where
    A: Space,
    B: Space,
    U: Space,
    V: Space,
    S: GEMV<F, U, A>,
    T: GEMV<F, V, B>,
    F: Num,
    Self: Linear<Pair<U, V>, Codomain = Pair<A, B>>,
{
    fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
        x.eval_decompose(|Pair(u, v)| {
            self.0.gemv(&mut y.0, α, u, β);
            self.1.gemv(&mut y.1, α, v, β);
        })
    }

    fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
        x.eval_decompose(|Pair(u, v)| {
            self.0.apply_mut(&mut y.0, u);
            self.1.apply_mut(&mut y.1, v);
        })
    }

    /// Computes `y += Ax`, where `A` is `Self`
    fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
        x.eval_decompose(|Pair(u, v)| {
            self.0.apply_add(&mut y.0, u);
            self.1.apply_add(&mut y.1, v);
        })
    }
}

impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
where
    A: Space,
    B: Space,
    Xʹ: Space,
    Yʹ: Space,
    R: ClosedSpace,
    S: Adjointable<A, Xʹ>,
    T: Adjointable<B, Yʹ>,
    Self: Linear<Pair<A, B>>,
    for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
{
    type AdjointCodomain = R;
    type Adjoint<'a>
        = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>
    where
        Self: 'a;

    fn adjoint(&self) -> Self::Adjoint<'_> {
        DiagOp(self.0.adjoint(), self.1.adjoint())
    }
}

impl<A, B, Xʹ, Yʹ, R, S, T> SimplyAdjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
where
    A: Space,
    B: Space,
    Xʹ: Space,
    Yʹ: Space,
    R: ClosedSpace,
    S: SimplyAdjointable<A, Xʹ>,
    T: SimplyAdjointable<B, Yʹ>,
    Self: Linear<Pair<A, B>>,
    for<'a> DiagOp<S::SimpleAdjoint, T::SimpleAdjoint>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
{
    type AdjointCodomain = R;
    type SimpleAdjoint = DiagOp<S::SimpleAdjoint, T::SimpleAdjoint>;

    fn adjoint(&self) -> Self::SimpleAdjoint {
        DiagOp(self.0.adjoint(), self.1.adjoint())
    }
}

impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
where
    A: Space,
    B: Space,
    Xʹ: Space,
    Yʹ: Space,
    R: ClosedSpace,
    S: Preadjointable<A, Xʹ>,
    T: Preadjointable<B, Yʹ>,
    Self: Linear<Pair<A, B>>,
    for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
{
    type PreadjointCodomain = R;
    type Preadjoint<'a>
        = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
    where
        Self: 'a;

    fn preadjoint(&self) -> Self::Preadjoint<'_> {
        DiagOp(self.0.preadjoint(), self.1.preadjoint())
    }
}

/// Block operator
pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;

macro_rules! pairnorm {
    ($expj:ty) => {
        impl<F, A, B, S, T, ExpA, ExpB, ExpR>
            BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> for RowOp<S, T>
        where
            F: Float,
            A: Space + Norm<ExpA, F>,
            B: Space + Norm<ExpB, F>,
            S: BoundedLinear<A, ExpA, ExpR, F>,
            T: BoundedLinear<B, ExpB, ExpR, F>,
            S::Codomain: Add<T::Codomain>,
            <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
            ExpA: NormExponent,
            ExpB: NormExponent,
            ExpR: NormExponent,
        {
            fn opnorm_bound(
                &self,
                PairNorm(expa, expb, _): PairNorm<ExpA, ExpB, $expj>,
                expr: ExpR,
            ) -> DynResult<F> {
                // An application of the triangle inequality bounds the norm by the maximum
                // of the individual norms. A simple observation shows this to be exact.
                let na = self.0.opnorm_bound(expa, expr)?;
                let nb = self.1.opnorm_bound(expb, expr)?;
                Ok(na.max(nb))
            }
        }

        impl<F, A, S, T, ExpA, ExpS, ExpT> BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
            for ColOp<S, T>
        where
            F: Float,
            A: Space + Norm<ExpA, F>,
            S: BoundedLinear<A, ExpA, ExpS, F>,
            T: BoundedLinear<A, ExpA, ExpT, F>,
            ExpA: NormExponent,
            ExpS: NormExponent,
            ExpT: NormExponent,
        {
            fn opnorm_bound(
                &self,
                expa: ExpA,
                PairNorm(exps, expt, _): PairNorm<ExpS, ExpT, $expj>,
            ) -> DynResult<F> {
                // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
                // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
                let ns = self.0.opnorm_bound(expa, exps)?;
                let nt = self.1.opnorm_bound(expa, expt)?;
                Ok(ns.max(nt))
            }
        }
    };
}

pairnorm!(L1);
pairnorm!(L2);
pairnorm!(Linfinity);

/// The simplest linear mapping, scaling by a scalar.
///
/// TODO: redefined/replace [`Weighted`] by composition with [`Scaled`].
pub struct Scaled<F: Float>(pub F);

impl<Domain, F> Mapping<Domain> for Scaled<F>
where
    F: Float,
    Domain: Space,
    Domain::Principal: Mul<F>,
    <Domain::Principal as Mul<F>>::Output: ClosedSpace,
{
    type Codomain = <Domain::Principal as Mul<F>>::Output;

    /// Compute the value of `self` at `x`.
    fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
        x.own() * self.0
    }
}

impl<Domain, F> Linear<Domain> for Scaled<F>
where
    F: Float,
    Domain: Space,
    Domain::Principal: Mul<F>,
    <Domain::Principal as Mul<F>>::Output: ClosedSpace,
{
}