Mercurial > repos > alg_tools / changeset

--- a/src/convex.rs	Mon Apr 28 23:16:56 2025 -0500
+++ b/src/convex.rs	Tue Apr 29 00:03:12 2025 -0500
@@ -2,27 +2,29 @@
 Some convex analysis basics
 */
 
-use std::marker::PhantomData;
-use crate::types::*;
+use crate::euclidean::Euclidean;
+use crate::instance::{DecompositionMut, Instance, InstanceMut};
+use crate::linops::{IdOp, Scaled};
 use crate::mapping::{Mapping, Space};
-use crate::linops::IdOp;
-use crate::instance::{Instance, InstanceMut, DecompositionMut};
+use crate::norms::*;
 use crate::operator_arithmetic::{Constant, Weighted};
-use crate::norms::*;
+use crate::types::*;
+use std::marker::PhantomData;
 
 /// Trait for convex mappings. Has no features, just serves as a constraint
 ///
 /// TODO: should constrain `Mapping::Codomain` to implement a partial order,
 /// but this makes everything complicated with little benefit.
-pub trait ConvexMapping<Domain : Space, F : Num = f64> : Mapping<Domain, Codomain = F>
-{}
+pub trait ConvexMapping<Domain: Space, F: Num = f64>: Mapping<Domain, Codomain = F> {}
 
 /// Trait for mappings with a Fenchel conjugate
 ///
 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
 /// not be convex.
-pub trait Conjugable<Domain : HasDual<F>, F : Num = f64> : Mapping<Domain> {
-    type Conjugate<'a> : ConvexMapping<Domain::DualSpace, F> where Self : 'a;
+pub trait Conjugable<Domain: HasDual<F>, F: Num = f64>: Mapping<Domain> {
+    type Conjugate<'a>: ConvexMapping<Domain::DualSpace, F>
+    where
+        Self: 'a;
 
     fn conjugate(&self) -> Self::Conjugate<'_>;
 }
@@ -31,12 +33,14 @@
 ///
 /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`],
 /// but a `Preconjugable` mapping has to be convex.
-pub trait Preconjugable<Domain, Predual, F : Num = f64> : ConvexMapping<Domain, F>
+pub trait Preconjugable<Domain, Predual, F: Num = f64>: ConvexMapping<Domain, F>
 where
-    Domain : Space,
-    Predual : HasDual<F>
+    Domain: Space,
+    Predual: HasDual<F>,
 {
-    type Preconjugate<'a> : Mapping<Predual> where Self : 'a;
+    type Preconjugate<'a>: Mapping<Predual>
+    where
+        Self: 'a;
 
     fn preconjugate(&self) -> Self::Preconjugate<'_>;
 }
@@ -45,53 +49,54 @@
 ///
 /// The conjugate type has to implement [`ConvexMapping`], but a `Conjugable` mapping need
 /// not be convex.
-pub trait Prox<Domain : Space> : Mapping<Domain> {
-    type Prox<'a> : Mapping<Domain, Codomain=Domain> where Self : 'a;
+pub trait Prox<Domain: Space>: Mapping<Domain> {
+    type Prox<'a>: Mapping<Domain, Codomain = Domain>
+    where
+        Self: 'a;
 
     /// Returns a proximal mapping with weight τ
-    fn prox_mapping(&self, τ : Self::Codomain) -> Self::Prox<'_>;
+    fn prox_mapping(&self, τ: Self::Codomain) -> Self::Prox<'_>;
 
     /// Calculate the proximal mapping with weight τ
-    fn prox<I : Instance<Domain>>(&self, τ : Self::Codomain, z : I) -> Domain {
+    fn prox<I: Instance<Domain>>(&self, τ: Self::Codomain, z: I) -> Domain {
         self.prox_mapping(τ).apply(z)
     }
 
     /// Calculate the proximal mapping with weight τ in-place
-    fn prox_mut<'b>(&self, τ : Self::Codomain, y : &'b mut Domain)
+    fn prox_mut<'b>(&self, τ: Self::Codomain, y: &'b mut Domain)
     where
-        &'b mut Domain : InstanceMut<Domain>,
-        Domain:: Decomp : DecompositionMut<Domain>,
-        for<'a> &'a Domain : Instance<Domain>,
+        &'b mut Domain: InstanceMut<Domain>,
+        Domain::Decomp: DecompositionMut<Domain>,
+        for<'a> &'a Domain: Instance<Domain>,
     {
         *y = self.prox(τ, &*y);
     }
 }
 
-
 /// Constraint to the unit ball of the norm described by `E`.
-pub struct NormConstraint<F : Float, E : NormExponent> {
-    radius : F,
-    norm : NormMapping<F, E>,
+pub struct NormConstraint<F: Float, E: NormExponent> {
+    radius: F,
+    norm: NormMapping<F, E>,
 }
 
 impl<Domain, E, F> ConvexMapping<Domain, F> for NormMapping<F, E>
 where
-    Domain : Space,
-    E : NormExponent,
-    F : Float,
-    Self : Mapping<Domain, Codomain=F>
-{}
-
+    Domain: Space,
+    E: NormExponent,
+    F: Float,
+    Self: Mapping<Domain, Codomain = F>,
+{
+}
 
 impl<F, E, Domain> Mapping<Domain> for NormConstraint<F, E>
 where
-    Domain : Space + Norm<F, E>,
-    F : Float,
-    E : NormExponent,
+    Domain: Space + Norm<F, E>,
+    F: Float,
+    E: NormExponent,
 {
     type Codomain = F;
 
-    fn apply<I : Instance<Domain>>(&self, d : I) -> F {
+    fn apply<I: Instance<Domain>>(&self, d: I) -> F {
         if d.eval(|x| x.norm(self.norm.exponent)) <= self.radius {
             F::ZERO
         } else {
@@ -102,68 +107,80 @@
 
 impl<Domain, E, F> ConvexMapping<Domain, F> for NormConstraint<F, E>
 where
-    Domain : Space,
-    E : NormExponent,
-    F : Float,
-    Self : Mapping<Domain, Codomain=F>
-{}
-
+    Domain: Space,
+    E: NormExponent,
+    F: Float,
+    Self: Mapping<Domain, Codomain = F>,
+{
+}
 
 impl<E, F, Domain> Conjugable<Domain, F> for NormMapping<F, E>
 where
-    E : HasDualExponent,
-    F : Float,
-    Domain : HasDual<F> + Norm<F, E> + Space,
-    <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp>
+    E: HasDualExponent,
+    F: Float,
+    Domain: HasDual<F> + Norm<F, E> + Space,
+    <Domain as HasDual<F>>::DualSpace: Norm<F, E::DualExp>,
 {
-    type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a;
+    type Conjugate<'a>
+        = NormConstraint<F, E::DualExp>
+    where
+        Self: 'a;
 
     fn conjugate(&self) -> Self::Conjugate<'_> {
         NormConstraint {
-            radius : F::ONE,
-            norm : self.exponent.dual_exponent().as_mapping()
+            radius: F::ONE,
+            norm: self.exponent.dual_exponent().as_mapping(),
         }
     }
 }
 
 impl<C, E, F, Domain> Conjugable<Domain, F> for Weighted<NormMapping<F, E>, C>
 where
-    C : Constant<Type = F>,
-    E : HasDualExponent,
-    F : Float,
-    Domain : HasDual<F> + Norm<F, E> + Space,
-    <Domain as HasDual<F>>::DualSpace : Norm<F, E::DualExp>
+    C: Constant<Type = F>,
+    E: HasDualExponent,
+    F: Float,
+    Domain: HasDual<F> + Norm<F, E> + Space,
+    <Domain as HasDual<F>>::DualSpace: Norm<F, E::DualExp>,
 {
-    type Conjugate<'a> = NormConstraint<F, E::DualExp> where Self : 'a;
+    type Conjugate<'a>
+        = NormConstraint<F, E::DualExp>
+    where
+        Self: 'a;
 
     fn conjugate(&self) -> Self::Conjugate<'_> {
         NormConstraint {
-            radius : self.weight.value(),
-            norm : self.base_fn.exponent.dual_exponent().as_mapping()
+            radius: self.weight.value(),
+            norm: self.base_fn.exponent.dual_exponent().as_mapping(),
         }
     }
 }
 
 impl<Domain, E, F> Prox<Domain> for NormConstraint<F, E>
 where
-    Domain : Space + Norm<F, E>,
-    E : NormExponent,
-    F : Float,
-    NormProjection<F, E> : Mapping<Domain, Codomain=Domain>,
+    Domain: Space + Norm<F, E>,
+    E: NormExponent,
+    F: Float,
+    NormProjection<F, E>: Mapping<Domain, Codomain = Domain>,
 {
-    type Prox<'a> = NormProjection<F, E> where Self : 'a;
+    type Prox<'a>
+        = NormProjection<F, E>
+    where
+        Self: 'a;
 
     #[inline]
-    fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> {
+    fn prox_mapping(&self, _τ: Self::Codomain) -> Self::Prox<'_> {
         assert!(self.radius >= F::ZERO);
-        NormProjection{ radius : self.radius, exponent : self.norm.exponent }
+        NormProjection {
+            radius: self.radius,
+            exponent: self.norm.exponent,
+        }
     }
 }
 
 /// Projection to the unit ball of the norm described by `E`.
-pub struct NormProjection<F : Float, E : NormExponent> {
-    radius : F,
-    exponent : E,
+pub struct NormProjection<F: Float, E: NormExponent> {
+    radius: F,
+    exponent: E,
 }
 
 /*
@@ -182,41 +199,42 @@
 
 impl<F, E, Domain> Mapping<Domain> for NormProjection<F, E>
 where
-    Domain : Space + Projection<F, E>,
-    F : Float,
-    E : NormExponent,
+    Domain: Space + Projection<F, E>,
+    F: Float,
+    E: NormExponent,
 {
     type Codomain = Domain;
 
-    fn apply<I : Instance<Domain>>(&self, d : I) -> Domain {
+    fn apply<I: Instance<Domain>>(&self, d: I) -> Domain {
         d.own().proj_ball(self.radius, self.exponent)
     }
 }
 
+/// The zero mapping
+pub struct Zero<Domain: Space, F: Num>(PhantomData<(Domain, F)>);
 
-/// The zero mapping
-pub struct Zero<Domain : Space, F : Num>(PhantomData<(Domain, F)>);
-
-impl<Domain : Space, F : Num> Zero<Domain, F> {
+impl<Domain: Space, F: Num> Zero<Domain, F> {
     pub fn new() -> Self {
         Zero(PhantomData)
     }
 }
 
-impl<Domain : Space, F : Num> Mapping<Domain> for Zero<Domain, F> {
+impl<Domain: Space, F: Num> Mapping<Domain> for Zero<Domain, F> {
     type Codomain = F;
 
     /// Compute the value of `self` at `x`.
-    fn apply<I : Instance<Domain>>(&self, _x : I) -> Self::Codomain {
+    fn apply<I: Instance<Domain>>(&self, _x: I) -> Self::Codomain {
         F::ZERO
     }
 }
 
-impl<Domain : Space, F : Num> ConvexMapping<Domain, F> for Zero<Domain, F> { }
-
+impl<Domain: Space, F: Num> ConvexMapping<Domain, F> for Zero<Domain, F> {}
 
-impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for Zero<Domain, F> {
-    type Conjugate<'a> = ZeroIndicator<Domain::DualSpace, F> where Self : 'a;
+impl<Domain: HasDual<F>, F: Float> Conjugable<Domain, F> for Zero<Domain, F> {
+    type Conjugate<'a>
+        = ZeroIndicator<Domain::DualSpace, F>
+    where
+        Self: 'a;
 
     #[inline]
     fn conjugate(&self) -> Self::Conjugate<'_> {
@@ -224,12 +242,15 @@
     }
 }
 
-impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F>
+impl<Domain, Predual, F: Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F>
 where
-    Domain : Space,
-    Predual : HasDual<F>
+    Domain: Space,
+    Predual: HasDual<F>,
 {
-    type Preconjugate<'a> = ZeroIndicator<Predual, F> where Self : 'a;
+    type Preconjugate<'a>
+        = ZeroIndicator<Predual, F>
+    where
+        Self: 'a;
 
     #[inline]
     fn preconjugate(&self) -> Self::Preconjugate<'_> {
@@ -237,38 +258,43 @@
     }
 }
 
-impl<Domain : Space + Clone, F : Num> Prox<Domain> for Zero<Domain, F> {
-    type Prox<'a> = IdOp<Domain> where Self : 'a;
+impl<Domain: Space + Clone, F: Num> Prox<Domain> for Zero<Domain, F> {
+    type Prox<'a>
+        = IdOp<Domain>
+    where
+        Self: 'a;
 
     #[inline]
-    fn prox_mapping(&self, _τ : Self::Codomain) -> Self::Prox<'_> {
+    fn prox_mapping(&self, _τ: Self::Codomain) -> Self::Prox<'_> {
         IdOp::new()
     }
 }
 
+/// The zero indicator
+pub struct ZeroIndicator<Domain: Space, F: Num>(PhantomData<(Domain, F)>);
 
-/// The zero indicator
-pub struct ZeroIndicator<Domain : Space, F : Num>(PhantomData<(Domain, F)>);
-
-impl<Domain : Space, F : Num> ZeroIndicator<Domain, F> {
+impl<Domain: Space, F: Num> ZeroIndicator<Domain, F> {
     pub fn new() -> Self {
         ZeroIndicator(PhantomData)
     }
 }
 
-impl<Domain : Normed<F>, F : Float> Mapping<Domain> for ZeroIndicator<Domain, F> {
+impl<Domain: Normed<F>, F: Float> Mapping<Domain> for ZeroIndicator<Domain, F> {
     type Codomain = F;
 
     /// Compute the value of `self` at `x`.
-    fn apply<I : Instance<Domain>>(&self, x : I) -> Self::Codomain {
+    fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
         x.eval(|x̃| if x̃.is_zero() { F::ZERO } else { F::INFINITY })
     }
 }
 
-impl<Domain : Normed<F>, F : Float> ConvexMapping<Domain, F> for ZeroIndicator<Domain, F> { }
+impl<Domain: Normed<F>, F: Float> ConvexMapping<Domain, F> for ZeroIndicator<Domain, F> {}
 
-impl<Domain : HasDual<F>, F : Float> Conjugable<Domain, F> for ZeroIndicator<Domain, F> {
-    type Conjugate<'a> = Zero<Domain::DualSpace, F> where Self : 'a;
+impl<Domain: HasDual<F>, F: Float> Conjugable<Domain, F> for ZeroIndicator<Domain, F> {
+    type Conjugate<'a>
+        = Zero<Domain::DualSpace, F>
+    where
+        Self: 'a;
 
     #[inline]
     fn conjugate(&self) -> Self::Conjugate<'_> {
@@ -276,15 +302,77 @@
     }
 }
 
-impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F>
+impl<Domain, Predual, F: Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F>
 where
-    Domain : Normed<F>,
-    Predual : HasDual<F>
+    Domain: Normed<F>,
+    Predual: HasDual<F>,
 {
-    type Preconjugate<'a> = Zero<Predual, F> where Self : 'a;
+    type Preconjugate<'a>
+        = Zero<Predual, F>
+    where
+        Self: 'a;
 
     #[inline]
     fn preconjugate(&self) -> Self::Preconjugate<'_> {
         Zero::new()
     }
 }
+
+/// The squared Euclidean norm divided by two
+pub struct Norm222<Domain: Space, F: Float>(PhantomData<(Domain, F)>);
+
+impl<Domain: Euclidean<F>, F: Float> Norm222<Domain, F> {
+    pub fn new() -> Self {
+        Norm222(PhantomData)
+    }
+}
+
+impl<Domain: Euclidean<F>, F: Float> Mapping<Domain> for Norm222<Domain, F> {
+    type Codomain = F;
+
+    /// Compute the value of `self` at `x`.
+    fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
+        x.eval(|z| z.norm2_squared() / F::TWO)
+    }
+}
+
+impl<Domain: Euclidean<F>, F: Float> ConvexMapping<Domain, F> for Norm222<Domain, F> {}
+
+impl<Domain: Euclidean<F>, F: Float> Conjugable<Domain, F> for Norm222<Domain, F> {
+    type Conjugate<'a>
+        = Self
+    where
+        Self: 'a;
+
+    #[inline]
+    fn conjugate(&self) -> Self::Conjugate<'_> {
+        Self::new()
+    }
+}
+
+impl<Domain: Euclidean<F>, F: Float> Preconjugable<Domain, Domain, F> for Norm222<Domain, F> {
+    type Preconjugate<'a>
+        = Self
+    where
+        Self: 'a;
+
+    #[inline]
+    fn preconjugate(&self) -> Self::Preconjugate<'_> {
+        Self::new()
+    }
+}
+
+impl<Domain, F> Prox<Domain> for Norm222<Domain, F>
+where
+    F: Float,
+    Domain: Euclidean<F, Output = Domain>,
+{
+    type Prox<'a>
+        = Scaled<F>
+    where
+        Self: 'a;
+
+    fn prox_mapping(&self, τ: F) -> Self::Prox<'_> {
+        Scaled(F::ONE / (F::ONE + τ))
+    }
+}

--- a/src/linops.rs	Mon Apr 28 23:16:56 2025 -0500
+++ b/src/linops.rs	Tue Apr 29 00:03:12 2025 -0500
@@ -10,6 +10,7 @@
 use numeric_literals::replace_float_literals;
 use serde::Serialize;
 use std::marker::PhantomData;
+use std::ops::Mul;
 
 /// Trait for linear operators on `X`.
 pub trait Linear<X: Space>: Mapping<X> {}
@@ -728,3 +729,28 @@
 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 + ClosedMul<F>,
+{
+    type Codomain = <Domain 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 + ClosedMul<F>,
+{
+}