# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1735653775 18000
# Node ID 848ecc05becf7975deec090896b0ccd5b900dc81
# Parent  9226980e45a704d75ae1b3c8c808e318647792af
More convexity, normed spaces, etc.

diff -r 9226980e45a7 -r 848ecc05becf src/convex.rs
--- a/src/convex.rs	Tue Dec 31 08:30:02 2024 -0500
+++ b/src/convex.rs	Tue Dec 31 09:02:55 2024 -0500
@@ -2,8 +2,10 @@
 Some convex analysis basics
 */
 
+use std::marker::PhantomData;
 use crate::types::*;
 use crate::mapping::{Mapping, Space};
+use crate::linops::IdOp;
 use crate::instance::{Instance, InstanceMut, DecompositionMut};
 use crate::norms::*;
 
@@ -11,16 +13,15 @@
 ///
 /// TODO: should constrain `Mapping::Codomain` to implement a partial order,
 /// but this makes everything complicated with little benefit.
-pub trait ConvexMapping<Domain : Space> : Mapping<Domain>
+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 : Space> : Mapping<Domain> {
-    type DualDomain : Space;
-    type Conjugate<'a> : ConvexMapping<Self::DualDomain> 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<'_>;
 }
@@ -28,10 +29,13 @@
 /// Trait for mappings with a Fenchel preconjugate
 ///
 /// In contrast to [`Conjugable`], the preconjugate need not implement [`ConvexMapping`],
-/// but a `Preconjugable` mapping has be convex.
-pub trait Preconjugable<Domain : Space> : ConvexMapping<Domain> {
-    type PredualDomain : Space;
-    type Preconjugate<'a> : Mapping<Self::PredualDomain> where Self : 'a;
+/// but a `Preconjugable` mapping has to be convex.
+pub trait Preconjugable<Domain, Predual, F : Num = f64> : ConvexMapping<Domain, F>
+where
+    Domain : Space,
+    Predual : HasDual<F>
+{
+    type Preconjugate<'a> : Mapping<Predual> where Self : 'a;
 
     fn preconjugate(&self) -> Self::Preconjugate<'_>;
 }
@@ -65,20 +69,13 @@
 
 pub struct NormConjugate<F : Float, E : NormExponent>(NormMapping<F, E>);
 
-impl<Domain, E, F> ConvexMapping<Domain> for 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> {}
-
-
-impl<Domain, E, F> ConvexMapping<Domain> for NormConjugate<F, E>
-where
-    Domain : Space,
-    E : NormExponent,
-    F : Float,
-    Self : Mapping<Domain, Codomain=F> {}
+    Self : Mapping<Domain, Codomain=F>
+{}
 
 
 impl<F, E, Domain> Mapping<Domain> for NormConjugate<F, E>
@@ -98,20 +95,121 @@
     }
 }
 
+impl<Domain, E, F> ConvexMapping<Domain, F> for NormConjugate<F, E>
+where
+    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>
+{
+    type Conjugate<'a> = NormConjugate<F, E::DualExp> where Self : 'a;
+
+    fn conjugate(&self) -> Self::Conjugate<'_> {
+        NormConjugate(self.exponent.dual_exponent().as_mapping())
+    }
+}
 
 
-impl<E, F, Domain> Conjugable<Domain> for NormMapping<F, E>
-where
-    E : NormExponent + Clone,
-    F : Float,
-    Domain : Norm<F, E> + Space,
-{
+/// The zero mapping
+pub struct Zero<Domain : Space, F : Num>(PhantomData<(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> {
+    type Codomain = F;
 
-    type DualDomain = Domain;
-    type Conjugate<'a> = NormConjugate<F, E> where Self : 'a;
+    /// Compute the value of `self` at `x`.
+    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 : 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<'_> {
-        NormConjugate(self.clone())
+        ZeroIndicator::new()
     }
 }
 
+impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for Zero<Domain, F>
+where
+    Domain : Space,
+    Predual : HasDual<F>
+{
+    type Preconjugate<'a> = ZeroIndicator<Predual, F> where Self : 'a;
+
+    #[inline]
+    fn preconjugate(&self) -> Self::Preconjugate<'_> {
+        ZeroIndicator::new()
+    }
+}
+
+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<'_> {
+        IdOp::new()
+    }
+}
+
+
+/// The zero indicator
+pub struct ZeroIndicator<Domain : Space, F : Num>(PhantomData<(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> {
+    type Codomain = F;
+
+    /// Compute the value of `self` at `x`.
+    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 : 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<'_> {
+        Zero::new()
+    }
+}
+
+impl<Domain, Predual, F : Float> Preconjugable<Domain, Predual, F> for ZeroIndicator<Domain, F>
+where
+    Domain : Normed<F>,
+    Predual : HasDual<F>
+{
+    type Preconjugate<'a> = Zero<Predual, F> where Self : 'a;
+
+    #[inline]
+    fn preconjugate(&self) -> Self::Preconjugate<'_> {
+        Zero::new()
+    }
+}
diff -r 9226980e45a7 -r 848ecc05becf src/direct_product.rs
--- a/src/direct_product.rs	Tue Dec 31 08:30:02 2024 -0500
+++ b/src/direct_product.rs	Tue Dec 31 09:02:55 2024 -0500
@@ -15,7 +15,7 @@
 use crate::mapping::Space;
 use crate::linops::AXPY;
 use crate::loc::Loc;
-use crate::norms::{Norm, PairNorm, NormExponent};
+use crate::norms::{Norm, PairNorm, NormExponent, Normed, HasDual, L2};
 
 #[derive(Debug,Clone,Copy,PartialEq,Eq,Serialize,Deserialize)]
 pub struct Pair<A, B> (pub A, pub B);
@@ -468,3 +468,29 @@
 }
 
 
+impl<F : Float, A, B> Normed<F> for Pair<A,B>
+where
+    A : Normed<F>,
+    B : Normed<F>,
+{
+    type NormExp = PairNorm<A::NormExp, B::NormExp, L2>;
+
+    #[inline]
+    fn norm_exponent(&self) -> Self::NormExp {
+        PairNorm(self.0.norm_exponent(), self.1.norm_exponent(), L2)
+    }
+
+    #[inline]
+    fn is_zero(&self) -> bool {
+        self.0.is_zero() && self.1.is_zero()
+    }
+}
+
+impl<F : Float, A, B> HasDual<F> for Pair<A,B>
+where
+    A : HasDual<F>,
+    B : HasDual<F>,
+
+{
+    type DualSpace = Pair<A::DualSpace, B::DualSpace>;
+}
diff -r 9226980e45a7 -r 848ecc05becf src/loc.rs
--- a/src/loc.rs	Tue Dec 31 08:30:02 2024 -0500
+++ b/src/loc.rs	Tue Dec 31 09:02:55 2024 -0500
@@ -578,6 +578,31 @@
     }
 }
 
+
+/// The default norm for `Loc` is [`L2`].
+impl<F : Float, const N : usize> Normed<F> for Loc<F, N> {
+    type NormExp = L2;
+
+    #[inline]
+    fn norm_exponent(&self) -> Self::NormExp {
+        L2
+    }
+
+    // #[inline]
+    // fn similar_origin(&self) -> Self {
+    //     [F::ZERO; N].into()
+    // }
+
+    #[inline]
+    fn is_zero(&self) -> bool {
+        self.norm2_squared() == F::ZERO
+    }
+}
+
+impl<F : Float, const N : usize> HasDual<F> for Loc<F, N> {
+    type DualSpace = Self;
+}
+
 impl<F : Float, const N : usize> Norm<F, L2> for Loc<F, N> {
     #[inline]
     fn norm(&self, _ : L2) -> F { self.norm2() }
@@ -588,6 +613,17 @@
     fn dist(&self, other : &Self, _ : L2) -> F { self.dist2(other) }
 }
 
+/* Implemented automatically as Euclidean.
+impl<F : Float, const N : usize> Projection<F, L2> for Loc<F, N> {
+    #[inline]
+    fn proj_ball_mut(&mut self, ρ : F, nrm : L2) {
+        let n = self.norm(nrm);
+        if n > ρ {
+            *self *= ρ/n;
+        }
+    }
+}*/
+
 impl<F : Float, const N : usize> Norm<F, L1> for Loc<F, N> {
     /// This implementation is not stabilised as it's meant to be used for very small vectors.
     /// Use [`nalgebra`] for larger vectors.
diff -r 9226980e45a7 -r 848ecc05becf src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Tue Dec 31 08:30:02 2024 -0500
+++ b/src/nalgebra_support.rs	Tue Dec 31 09:02:55 2024 -0500
@@ -101,6 +101,20 @@
     }
 }
 
+/* Implemented automatically as Euclidean.
+impl<SM,M,E> Projection<E, L2> for Vector<E,M,SM>
+where SM: StorageMut<E,M> + Clone,
+      M : Dim, E : Scalar + Zero + One + Float + RealField,
+      DefaultAllocator : Allocator<M> {
+    #[inline]
+    fn proj_ball_mut(&mut self, ρ : E, _ : L2) {
+        let n = self.norm(L2);
+        if n > ρ {
+            self.iter_mut().for_each(|v| *v *= ρ/n)
+        }
+    }
+}*/
+
 impl<SM,M,E> Projection<E, Linfinity> for Vector<E,M,SM>
 where SM: StorageMut<E,M> + Clone,
       M : Dim, E : Scalar + Zero + One + Float + RealField,
@@ -205,10 +219,41 @@
     }
 }
 
+/// The default norm for `Vector` is [`L2`].
+impl<E,M,S> Normed<E>
+for Vector<E,M,S>
+where M : Dim,
+      S : Storage<E,M> + Clone,
+      E : Float + Scalar + Zero + One + RealField,
+      DefaultAllocator : Allocator<M> {
+
+    type NormExp = L2;
+
+    #[inline]
+    fn norm_exponent(&self) -> Self::NormExp {
+        L2
+    }
+
+    #[inline]
+    fn is_zero(&self) -> bool {
+        Vector::<E,M,S>::norm_squared(self) == E::ZERO
+    }
+}
+
+impl<E,M,S> HasDual<E>
+for Vector<E,M,S>
+where M : Dim,
+      S : Storage<E,M> + Clone,
+      E : Float + Scalar + Zero + One + RealField,
+      DefaultAllocator : Allocator<M> {
+    // TODO: Doesn't work with different storage formats.
+    type DualSpace = Self;
+}
+
 impl<E,M,S> Norm<E, L1>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : Storage<E,M>,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
 
@@ -221,7 +266,7 @@
 impl<E,M,S> Dist<E, L1>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : Storage<E,M>,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
     #[inline]
@@ -233,7 +278,7 @@
 impl<E,M,S> Norm<E, L2>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : Storage<E,M>,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
 
@@ -246,7 +291,7 @@
 impl<E,M,S> Dist<E, L2>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : Storage<E,M>,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
     #[inline]
@@ -258,7 +303,7 @@
 impl<E,M,S> Norm<E, Linfinity>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : Storage<E,M>,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
 
@@ -271,7 +316,7 @@
 impl<E,M,S> Dist<E, Linfinity>
 for Vector<E,M,S>
 where M : Dim,
-      S : StorageMut<E,M>,
+      S : Storage<E,M>,
       E : Float + Scalar + Zero + One + RealField,
       DefaultAllocator : Allocator<M> {
     #[inline]
diff -r 9226980e45a7 -r 848ecc05becf src/norms.rs
--- a/src/norms.rs	Tue Dec 31 08:30:02 2024 -0500
+++ b/src/norms.rs	Tue Dec 31 09:02:55 2024 -0500
@@ -119,7 +119,7 @@
 ///
 /// println!("{:?}, {:?}", x.proj_ball(1.0, L2), x.proj_ball(0.5, Linfinity));
 /// ```
-pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Euclidean<F>
+pub trait Projection<F : Num, Exponent : NormExponent> : Norm<F, Exponent> + Sized
 where F : Float {
     /// Projection of `self` to the `q`-norm-ball of radius ρ.
     fn proj_ball(mut self, ρ : F, q : Exponent) -> Self {
@@ -128,7 +128,7 @@
     }
 
     /// In-place projection of `self` to the `q`-norm-ball of radius ρ.
-    fn proj_ball_mut(&mut self, ρ : F, _q : Exponent);
+    fn proj_ball_mut(&mut self, ρ : F, q : Exponent);
 }
 
 /*impl<F : Float, E : Euclidean<F>> Norm<F, L2> for E {
@@ -202,3 +202,65 @@
     }
 }
 
+pub trait Normed<F : Num = f64> : Space + Norm<F, Self::NormExp> {
+    type NormExp : NormExponent;
+
+    fn norm_exponent(&self) -> Self::NormExp;
+
+    #[inline]
+    fn norm_(&self) -> F {
+        self.norm(self.norm_exponent())
+    }
+
+    // fn similar_origin(&self) -> Self;
+
+    fn is_zero(&self) -> bool;
+}
+
+pub trait HasDual<F : Num = f64> : Normed<F> {
+    type DualSpace : Normed<F>;
+}
+
+pub trait Reflexive<F : Num = f64> : HasDual<F>
+where
+    Self::DualSpace : HasDual<F, DualSpace = Self>
+{ }
+
+
+pub trait HasDualExponent : NormExponent {
+    type DualExp : NormExponent;
+
+    fn dual_exponent(&self) -> Self::DualExp;
+}
+
+impl HasDualExponent for L2 {
+    type DualExp = L2;
+    
+    #[inline]
+    fn dual_exponent(&self) -> Self::DualExp {
+        L2
+    }
+}
+
+impl HasDualExponent for L1 {
+    type DualExp = Linfinity;
+    
+    #[inline]
+    fn dual_exponent(&self) -> Self::DualExp {
+        Linfinity
+    }
+}
+
+
+impl HasDualExponent for Linfinity {
+    type DualExp = L1;
+    
+    #[inline]
+    fn dual_exponent(&self) -> Self::DualExp {
+        L1
+    }
+}
+
+
+
+