# HG changeset patch
# User Tuomo Valkonen <tuomov@iki.fi>
# Date 1747082568 18000
# Node ID d2994e34a5f55ef3dde7ec7f0a72626652b21ae3
# Parent  f75bf34adda0b8df216f055288638f5e00ebae98
Simplify ZeroOp to SimpleZeroOp, only from X to X. Add Prox for ZeroIndicator. Move F parameter to AXPY::Field.

diff -r f75bf34adda0 -r d2994e34a5f5 src/convex.rs
--- a/src/convex.rs	Sun May 11 02:03:45 2025 -0500
+++ b/src/convex.rs	Mon May 12 15:42:48 2025 -0500
@@ -5,7 +5,7 @@
 use crate::error::DynResult;
 use crate::euclidean::Euclidean;
 use crate::instance::{DecompositionMut, Instance, InstanceMut};
-use crate::linops::{IdOp, Scaled};
+use crate::linops::{IdOp, Scaled, SimpleZeroOp, AXPY};
 use crate::mapping::{DifferentiableImpl, LipschitzDifferentiableImpl, Mapping, Space};
 use crate::norms::*;
 use crate::operator_arithmetic::{Constant, Weighted};
@@ -333,6 +333,22 @@
     }
 }
 
+impl<Domain, F> Prox<Domain> for ZeroIndicator<Domain, F>
+where
+    Domain: AXPY<Field = F, Owned = Domain> + Normed<F>,
+    F: Float,
+{
+    type Prox<'a>
+        = SimpleZeroOp
+    where
+        Self: 'a;
+
+    /// Returns a proximal mapping with weight τ
+    fn prox_mapping(&self, _τ: F) -> Self::Prox<'_> {
+        return SimpleZeroOp;
+    }
+}
+
 /// The squared Euclidean norm divided by two
 #[derive(Copy, Clone, Serialize, Deserialize)]
 pub struct Norm222<F: Float>(PhantomData<F>);
diff -r f75bf34adda0 -r d2994e34a5f5 src/direct_product.rs
--- a/src/direct_product.rs	Sun May 11 02:03:45 2025 -0500
+++ b/src/direct_product.rs	Mon May 12 15:42:48 2025 -0500
@@ -244,13 +244,13 @@
 impl_binary_mut!(SubAssign, sub_assign);
 
 macro_rules! impl_scalar_mut {
-    ($trait:ident, $fn:ident, $F:ty) => {
-        impl<'a, A, B> $trait<$F> for Pair<A, B>
+    ($trait:ident, $fn:ident) => {
+        impl<'a, A, B, F: Num> $trait<F> for Pair<A, B>
         where
-            A: $trait<$F>,
-            B: $trait<$F>,
+            A: $trait<F>,
+            B: $trait<F>,
         {
-            fn $fn(&mut self, t: $F) {
+            fn $fn(&mut self, t: F) {
                 let Pair(ref mut a, ref mut b) = self;
                 a.$fn(t);
                 b.$fn(t);
@@ -259,10 +259,8 @@
     };
 }
 
-impl_scalar_mut!(MulAssign, mul_assign, f32);
-impl_scalar_mut!(MulAssign, mul_assign, f64);
-impl_scalar_mut!(DivAssign, div_assign, f32);
-impl_scalar_mut!(DivAssign, div_assign, f64);
+impl_scalar_mut!(MulAssign, mul_assign);
+impl_scalar_mut!(DivAssign, div_assign);
 
 /// We only support 'closed' `Euclidean` `Pair`s, as more general ones cause
 /// compiler overflows.
@@ -303,17 +301,20 @@
     }
 }
 
-impl<F, A, B, U, V> AXPY<F, Pair<U, V>> for Pair<A, B>
+impl<F, A, B, U, V> AXPY<Pair<U, V>> for Pair<A, B>
 where
     U: Space,
     V: Space,
-    A: AXPY<F, U>,
-    B: AXPY<F, V>,
+    A: AXPY<U, Field = F>,
+    B: AXPY<V, Field = F>,
     F: Num,
     Self: MulAssign<F>,
     Pair<A, B>: MulAssign<F>,
-    Pair<A::Owned, B::Owned>: AXPY<F, Pair<U, V>>,
+    //A::Owned: MulAssign<F>,
+    //B::Owned: MulAssign<F>,
+    //Pair<A::Owned, B::Owned>: AXPY<Pair<U, V>, Field = F>,
 {
+    type Field = F;
     type Owned = Pair<A::Owned, B::Owned>;
 
     fn axpy<I: Instance<Pair<U, V>>>(&mut self, α: F, x: I, β: F) {
diff -r f75bf34adda0 -r d2994e34a5f5 src/linops.rs
--- a/src/linops.rs	Sun May 11 02:03:45 2025 -0500
+++ b/src/linops.rs	Mon May 12 15:42:48 2025 -0500
@@ -6,7 +6,7 @@
 use crate::error::DynResult;
 use crate::instance::Instance;
 pub use crate::mapping::{Composition, DifferentiableImpl, Mapping, Space};
-use crate::norms::{Linfinity, Norm, NormExponent, PairNorm, L1, L2};
+use crate::norms::{HasDual, Linfinity, Norm, NormExponent, PairNorm, L1, L2};
 use crate::types::*;
 use numeric_literals::replace_float_literals;
 use serde::Serialize;
@@ -26,16 +26,19 @@
 // }
 
 /// Efficient in-place summation.
-#[replace_float_literals(F::cast_from(literal))]
-pub trait AXPY<F, X = Self>: Space + std::ops::MulAssign<F>
+#[replace_float_literals(Self::Field::cast_from(literal))]
+pub trait AXPY<X = Self>: Space + std::ops::MulAssign<Self::Field>
 where
-    F: Num,
     X: Space,
 {
-    type Owned: AXPY<F, X>;
+    type Field: Num;
+    type Owned: AXPY<X, Field = Self::Field>;
+    // type OriginGen<'a>: OriginGenerator<Self::Owned, F>
+    // where
+    //     Self: 'a;
 
     /// Computes  `y = βy + αx`, where `y` is `Self`.
-    fn axpy<I: Instance<X>>(&mut self, α: F, x: I, β: F);
+    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) {
@@ -43,15 +46,25 @@
     }
 
     /// Computes  `y = αx`, where `y` is `Self`.
-    fn scale_from<I: Instance<X>>(&mut self, α: F, x: I) {
+    fn scale_from<I: Instance<X>>(&mut self, α: Self::Field, x: I) {
         self.axpy(α, x, 0.0)
     }
 
     /// Return a similar zero as `self`.
     fn similar_origin(&self) -> Self::Owned;
+    // {
+    //     self.make_origin_generator().make_origin()
+    // }
+
+    /// Return a similar zero as `x`.
+    fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::Owned {
+        x.eval(|xr| xr.similar_origin())
+    }
 
     /// Set self to zero.
     fn set_zero(&mut self);
+
+    //fn make_origin_generator(&self) -> Self::OriginGen<'_>;
 }
 
 /// Efficient in-place application for [`Linear`] operators.
@@ -167,7 +180,7 @@
 #[replace_float_literals(F::cast_from(literal))]
 impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X>
 where
-    Y: AXPY<F, X>,
+    Y: AXPY<X, Field = F>,
     X: Clone + Space,
 {
     // Computes  `y = αAx + βy`, where `A` is `Self`.
@@ -215,42 +228,29 @@
     }
 }
 
-/// The zero operator
+/// The zero operator from a space to itself
 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
-pub struct ZeroOp<'a, X, XD, Y, F> {
-    zero: &'a Y, // TODO: don't pass this in `new`; maybe not even store.
-    dual_or_predual_zero: XD,
-    _phantoms: PhantomData<(X, Y, F)>,
-}
+pub struct SimpleZeroOp;
 
-// TODO: Need to make Zero in Instance.
+impl<X> Mapping<X> for SimpleZeroOp
+where
+    X: AXPY + Instance<X>,
+{
+    type Codomain = X::Owned;
 
-impl<'a, F: Num, X: Space, XD, Y: Space + Clone> ZeroOp<'a, X, XD, Y, F> {
-    pub fn new(zero: &'a Y, dual_or_predual_zero: XD) -> Self {
-        ZeroOp {
-            zero,
-            dual_or_predual_zero,
-            _phantoms: PhantomData,
-        }
+    fn apply<I: Instance<X>>(&self, x: I) -> X::Owned {
+        X::similar_origin_inst(x)
     }
 }
 
-impl<'a, F: Num, X: Space, XD, Y: AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> {
-    type Codomain = Y;
-
-    fn apply<I: Instance<X>>(&self, _x: I) -> Y {
-        self.zero.clone()
-    }
-}
-
-impl<'a, F: Num, X: Space, XD, Y: AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F> {}
+impl<X> Linear<X> for SimpleZeroOp where X: AXPY + Instance<X> {}
 
 #[replace_float_literals(F::cast_from(literal))]
-impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F>
+impl<X, Y, F> GEMV<F, X, Y> for SimpleZeroOp
 where
     F: Num,
-    Y: AXPY<F, Y> + Clone,
-    X: Space,
+    Y: AXPY<Field = F>,
+    X: AXPY<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) {
@@ -262,11 +262,10 @@
     }
 }
 
-impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F>
+impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp
 where
-    X: Space + Norm<E1, F>,
-    Y: AXPY<F> + Clone + Norm<E2, F>,
     F: Num,
+    X: AXPY<Field = F> + Instance<X> + Norm<E1, F>,
     E1: NormExponent,
     E2: NormExponent,
 {
@@ -275,43 +274,184 @@
     }
 }
 
-impl<'a, F: Num, X, XD, Y, Yprime: Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F>
+impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp
 where
-    X: Space,
-    Y: AXPY<F> + Clone + 'static,
-    XD: AXPY<F> + Clone + 'static,
+    F: Num,
+    X: AXPY<Field = F> + Instance<X> + HasDual<F>,
+    X::DualSpace: AXPY<Owned = X::DualSpace>,
 {
-    type AdjointCodomain = XD;
+    type AdjointCodomain = X::DualSpace;
     type Adjoint<'b>
-        = ZeroOp<'b, Yprime, (), XD, F>
+        = SimpleZeroOp
     where
         Self: 'b;
     // () means not (pre)adjointable.
 
     fn adjoint(&self) -> Self::Adjoint<'_> {
-        ZeroOp::new(&self.dual_or_predual_zero, ())
+        SimpleZeroOp
+    }
+}
+
+/*
+/// Trait for forming origins (zeroes) in vector spaces
+pub trait OriginGenerator<X, F: Num = f64> {
+    fn make_origin(&self) -> X;
+}
+
+/// Origin generator for statically sized Euclidean spaces.
+pub struct StaticEuclideanOriginGenerator;
+
+impl<X, F: Float> OriginGenerator<X, F> for StaticEuclideanOriginGenerator
+where
+    X: StaticEuclidean<F> + Euclidean<F, Output = X>,
+{
+    #[inline]
+    fn make_origin(&self) -> X {
+        X::origin()
+    }
+}
+
+/// Origin generator arbitrary spaces that implement [`AXPY`], based on a sample vector.
+pub struct SimilarOriginGenerator<'a, X>(&'a X);
+
+impl<'a, X, F: Float> OriginGenerator<X, F> for SimilarOriginGenerator<'a, X>
+where
+    X: AXPY<F, Owned = X>,
+{
+    #[inline]
+    fn make_origin(&self) -> X {
+        self.0.similar_origin()
+    }
+}
+
+pub struct NotAnOriginGenerator;
+
+/// The zero operator
+#[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
+pub struct ZeroOp<X, Y, YOrigin, XDOrigin = NotAnOriginGenerator, F: Num = f64> {
+    y_origin: YOrigin,
+    xd_origin: XDOrigin,
+    _phantoms: PhantomData<(X, Y, F)>,
+}
+
+// TODO: Need to make Zero in Instance.
+
+impl<X, Y, F, YOrigin, XDOrigin> ZeroOp<X, Y, YOrigin, XDOrigin, F>
+where
+    F: Num,
+    YOrigin: OriginGenerator<Y, F>,
+{
+    pub fn new(y_origin: YOrigin, xd_origin: XDOrigin) -> Self {
+        ZeroOp {
+            y_origin,
+            xd_origin,
+            _phantoms: PhantomData,
+        }
     }
 }
 
-impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F>
+impl<X, Y, F, YOrigin, XDOrigin> Mapping<X> for ZeroOp<X, Y, YOrigin, XDOrigin, F>
+where
+    F: Num,
+    YOrigin: OriginGenerator<Y, F>,
+    Y: Space,
+    X: Space + Instance<X>,
+{
+    type Codomain = Y;
+
+    fn apply<I: Instance<X>>(&self, _: I) -> Y {
+        self.y_origin.make_origin()
+    }
+}
+
+impl<X, Y, F, YOrigin, XDOrigin> Linear<X> for ZeroOp<X, Y, YOrigin, XDOrigin, F>
+where
+    F: Num,
+    YOrigin: OriginGenerator<Y, F>,
+    Y: Space,
+    X: Space + Instance<X>,
+{
+}
+
+#[replace_float_literals(F::cast_from(literal))]
+impl<X, Y, F, YOrigin, XDOrigin> GEMV<F, X, Y> for ZeroOp<X, Y, YOrigin, XDOrigin, F>
 where
     F: Num,
-    X: Space,
-    Y: AXPY<F> + Clone,
-    Ypre: Space,
-    XD: AXPY<F> + Clone + 'static,
+    YOrigin: OriginGenerator<Y, F>,
+    Y: Space + AXPY<F>,
+    X: Space + Instance<X>,
 {
-    type PreadjointCodomain = XD;
-    type Preadjoint<'b>
-        = ZeroOp<'b, Ypre, (), XD, F>
+    // 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, F, YOrigin, XDOrigin, E1, E2> BoundedLinear<X, E1, E2, F>
+    for ZeroOp<X, Y, YOrigin, XDOrigin, F>
+where
+    F: Num,
+    YOrigin: OriginGenerator<Y, F>,
+    Y: Space + AXPY<F> + Norm<E2, F>,
+    X: Space + Instance<X> + Norm<E1, F>,
+    E1: NormExponent,
+    E2: NormExponent,
+{
+    fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
+        Ok(F::ZERO)
+    }
+}
+
+impl<X, Y, Yprime, F, YOrigin, XDOrigin> Adjointable<X, Yprime>
+    for ZeroOp<X, Y, YOrigin, XDOrigin, F>
+where
+    F: Num,
+    YOrigin: OriginGenerator<Y, F>,
+    Y: Space + AXPY<F>,
+    X: Space + Instance<X> + HasDual<F>,
+    XDOrigin: OriginGenerator<X::DualSpace, F> + Clone,
+    Yprime: Space + Instance<Yprime>,
+{
+    type AdjointCodomain = X::DualSpace;
+    type Adjoint<'b>
+        = ZeroOp<Yprime, X::DualSpace, XDOrigin, NotAnOriginGenerator, F>
     where
         Self: 'b;
     // () means not (pre)adjointable.
 
-    fn preadjoint(&self) -> Self::Preadjoint<'_> {
-        ZeroOp::new(&self.dual_or_predual_zero, ())
+    fn adjoint(&self) -> Self::Adjoint<'_> {
+        ZeroOp::new(self.xd_origin.clone(), NotAnOriginGenerator)
     }
 }
+*/
+
+/*
+impl<X, Y, Ypre, Xpre, F, YOrigin, XDOrigin> Preadjointable<X, Ypre>
+    for ZeroOp<X, Y, YOrigin, XDOrigin, F>
+where
+    F: Num,
+    YOrigin: OriginGenerator<Y, F>,
+    Y: Space + AXPY<F>,
+    X: Space + Instance<X> + HasDual<F>,
+    XDOrigin: OriginGenerator<Xpre, F> + Clone,
+    Ypre: Space + Instance<Ypre> + HasDual<DualSpace = X>,
+{
+    type PreadjointCodomain = Xpre;
+    type Preadjoint<'b>
+        = ZeroOp<Ypre, Xpre, XDOrigin, NotAnOriginGenerator, F>
+    where
+        Self: 'b;
+    // () means not (pre)adjointable.
+
+    fn adjoint(&self) -> Self::Adjoint<'_> {
+        ZeroOp::new(self.xpre_origin.clone(), NotAnOriginGenerator)
+    }
+}
+*/
 
 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
 where
diff -r f75bf34adda0 -r d2994e34a5f5 src/loc.rs
--- a/src/loc.rs	Sun May 11 02:03:45 2025 -0500
+++ b/src/loc.rs	Mon May 12 15:42:48 2025 -0500
@@ -710,7 +710,8 @@
 
 impl<F: Float, const N: usize> Linear<Loc<N, F>> for Loc<N, F> {}
 
-impl<F: Float, const N: usize> AXPY<F, Loc<N, F>> for Loc<N, F> {
+impl<F: Float, const N: usize> AXPY<Loc<N, F>> for Loc<N, F> {
+    type Field = F;
     type Owned = Self;
 
     #[inline]
@@ -729,12 +730,22 @@
         x.eval(|x̃| map2_mut(self, x̃, |yi, xi| *yi = *xi))
     }
 
+    // #[inline]
+    // fn make_origin_generator(&self) -> StaticEuclideanOriginGenerator {
+    //     StaticEuclideanOriginGenerator
+    // }
+
     #[inline]
     fn similar_origin(&self) -> Self::Owned {
         Self::ORIGIN
     }
 
     #[inline]
+    fn similar_origin_inst<I: Instance<Self>>(_: I) -> Self::Owned {
+        Self::ORIGIN
+    }
+
+    #[inline]
     fn set_zero(&mut self) {
         *self = Self::ORIGIN;
     }
diff -r f75bf34adda0 -r d2994e34a5f5 src/nalgebra_support.rs
--- a/src/nalgebra_support.rs	Sun May 11 02:03:45 2025 -0500
+++ b/src/nalgebra_support.rs	Mon May 12 15:42:48 2025 -0500
@@ -103,7 +103,7 @@
     }
 }
 
-impl<SM, SV1, M, E> AXPY<E, Vector<E, M, SV1>> for Vector<E, M, SM>
+impl<SM, SV1, M, E> AXPY<Vector<E, M, SV1>> for Vector<E, M, SM>
 where
     SM: StorageMut<E, M> + Clone,
     SV1: Storage<E, M> + Clone,
@@ -111,6 +111,7 @@
     E: Scalar + Zero + One + Float,
     DefaultAllocator: Allocator<M>,
 {
+    type Field = E;
     type Owned = OVector<E, M>;
 
     #[inline]