--- 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