--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/mapping/dataterm.rs Fri May 15 14:46:30 2026 -0500 @@ -0,0 +1,155 @@ +/*! +General deata terms of the form $g(Ax-b)$ for an operator $A$ +to a [`Euclidean`] space, and a function g on that space. +*/ + +#![allow(non_snake_case)] + +use super::{DifferentiableImpl, DifferentiableMapping, LipschitzDifferentiableImpl, Mapping}; +use crate::convex::ConvexMapping; +use crate::error::DynResult; +use crate::instance::{ClosedSpace, Instance, Space}; +use crate::linops::{BoundedLinear, Linear, Preadjointable}; +use crate::norms::{Normed, L2}; +use crate::types::Float; +use std::ops::Sub; + +//use serde::{Deserialize, Serialize}; + +/// Functions of the form $g(Ax-b)$ for an operator $A$, data $b$, and fidelity $g$. +pub struct DataTerm< + F: Float, + Domain: Space, + A: Mapping<Domain>, + G: Mapping<A::Codomain, Codomain = F>, +> { + // The operator A + opA: A, + // The data b + b: <A as Mapping<Domain>>::Codomain, + // The outer fidelity + g: G, +} + +// Derive has troubles with `b`. +impl<F, Domain, A, G> Clone for DataTerm<F, Domain, A, G> +where + F: Float, + Domain: Space, + A: Mapping<Domain> + Clone, + <A as Mapping<Domain>>::Codomain: Clone, + G: Mapping<A::Codomain, Codomain = F> + Clone, +{ + fn clone(&self) -> Self { + DataTerm { opA: self.opA.clone(), b: self.b.clone(), g: self.g.clone() } + } +} + +#[allow(non_snake_case)] +impl<F: Float, Domain: Space, A: Mapping<Domain>, G: Mapping<A::Codomain, Codomain = F>> + DataTerm<F, Domain, A, G> +{ + pub fn new(opA: A, b: A::Codomain, g: G) -> Self { + DataTerm { opA, b, g } + } + + pub fn operator(&self) -> &'_ A { + &self.opA + } + + pub fn data(&self) -> &'_ <A as Mapping<Domain>>::Codomain { + &self.b + } + + pub fn fidelity(&self) -> &'_ G { + &self.g + } + + /// Returns the residual $Ax-b$. + pub fn residual<'a, 'b>(&'b self, x: &'a Domain) -> <A as Mapping<Domain>>::Codomain + where + &'a Domain: Instance<Domain>, + <A as Mapping<Domain>>::Codomain: + Sub<&'b <A as Mapping<Domain>>::Codomain, Output = <A as Mapping<Domain>>::Codomain>, + { + self.opA.apply(x) - &self.b + } +} + +//+ AdjointProductBoundedBy<RNDM<N, F>, P, FloatType = F>, + +impl<F, X, A, G> Mapping<X> for DataTerm<F, X, A, G> +where + F: Float, + X: Space, + A: Mapping<X>, + G: Mapping<A::Codomain, Codomain = F>, + A::Codomain: ClosedSpace + for<'a> Sub<&'a A::Codomain, Output = A::Codomain>, +{ + type Codomain = F; + + fn apply<I: Instance<X>>(&self, x: I) -> F { + // TODO: possibly (if at all more effcient) use GEMV once generalised + // to not require preallocation. However, Rust should be pretty efficient + // at not doing preallocations or anything here, as the result of self.opA.apply() + // can be consumed, so maybe GEMV is no use. + self.g.apply(self.opA.apply(x) - &self.b) + } +} + +impl<F, X, A, G> ConvexMapping<X, F> for DataTerm<F, X, A, G> +where + F: Float, + X: Normed<F>, + A: Linear<X>, + G: ConvexMapping<A::Codomain, F>, + A::Codomain: ClosedSpace + Normed<F> + for<'a> Sub<&'a A::Codomain, Output = A::Codomain>, +{ +} + +impl<F, X, Y, A, G> DifferentiableImpl<X> for DataTerm<F, X, A, G> +where + F: Float, + X: Space, + Y: Space + Instance<Y> + for<'a> Sub<&'a Y, Output = Y>, + A: Linear<X, Codomain = Y> + Preadjointable<X, G::DerivativeDomain>, + G::DerivativeDomain: Instance<G::DerivativeDomain>, + A::PreadjointCodomain: ClosedSpace, + G: DifferentiableMapping<Y, Codomain = F>, + Self: Mapping<X, Codomain = F>, +{ + type Derivative = A::PreadjointCodomain; + + fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative { + // TODO: possibly (if at all more effcient) use GEMV once generalised + // to not require preallocation. However, Rust should be pretty efficient + // at not doing preallocations or anything here, as the result of self.opA.apply() + // can be consumed, so maybe GEMV is no use. + //self.opA.preadjoint().apply(self.opA.apply(x) - self.b) + self.opA + .preadjoint() + .apply(self.g.differential(self.opA.apply(x) - &self.b)) + } + + fn apply_and_differential_impl<I: Instance<X>>(&self, x: I) -> (F, Self::Derivative) { + let j = self.opA.apply(x) - &self.b; + let (v, d) = self.g.apply_and_differential(j); + (v, self.opA.preadjoint().apply(d)) + } +} + +impl<'a, F, X, Y, A, G> LipschitzDifferentiableImpl<X, X::NormExp> for DataTerm<F, X, A, G> +where + F: Float, + X: Normed<F>, + Y: Normed<F>, + A: BoundedLinear<X, X::NormExp, L2, F, Codomain = Y>, + G: Mapping<Y, Codomain = F> + LipschitzDifferentiableImpl<Y, Y::NormExp>, + Self: DifferentiableImpl<X>, +{ + type FloatType = F; + + fn diff_lipschitz_factor(&self, seminorm: X::NormExp) -> DynResult<F> { + Ok(self.opA.opnorm_bound(seminorm, L2)?.powi(2)) + } +}