Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

changeset 198
3868555d135c
parent 197
1f301affeae3
equal deleted inserted replaced
94:1f19c6bbf07b 198:3868555d135c
1 /*! 1 /*!
2 Abstract linear operators. 2 Abstract linear operators.
3 */ 3 */
4 4
5 use crate::direct_product::Pair;
6 use crate::error::DynResult;
7 use crate::euclidean::StaticEuclidean;
8 use crate::instance::Instance;
9 pub use crate::mapping::{ClosedSpace, Composition, DifferentiableImpl, Mapping, Space};
10 use crate::norms::{HasDual, Linfinity, NormExponent, PairNorm, L1, L2};
11 use crate::types::*;
5 use numeric_literals::replace_float_literals; 12 use numeric_literals::replace_float_literals;
13 use serde::Serialize;
6 use std::marker::PhantomData; 14 use std::marker::PhantomData;
7 use serde::Serialize; 15 use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
8 use crate::types::*;
9 pub use crate::mapping::{Mapping, Space, Composition};
10 use crate::direct_product::Pair;
11 use crate::instance::Instance;
12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm};
13 16
14 /// Trait for linear operators on `X`. 17 /// Trait for linear operators on `X`.
15 pub trait Linear<X : Space> : Mapping<X> 18 pub trait Linear<X: Space>: Mapping<X> {}
16 { } 19
20 // impl<X: Space, A: Linear<X>> DifferentiableImpl<X> for A {
21 // type Derivative = <Self as Mapping<X>>::Codomain;
22
23 // /// Compute the differential of `self` at `x`, consuming the input.
24 // fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
25 // self.apply(x)
26 // }
27 // }
28
29 /// Vector spaces
30 #[replace_float_literals(Self::Field::cast_from(literal))]
31 pub trait VectorSpace:
32 Space<Principal = Self::PrincipalV>
33 + Mul<Self::Field, Output = Self::PrincipalV>
34 + Div<Self::Field, Output = Self::PrincipalV>
35 + Add<Self, Output = Self::PrincipalV>
36 + Add<Self::PrincipalV, Output = Self::PrincipalV>
37 + Sub<Self, Output = Self::PrincipalV>
38 + Sub<Self::PrincipalV, Output = Self::PrincipalV>
39 + Neg
40 + for<'b> Add<&'b Self, Output = <Self as VectorSpace>::PrincipalV>
41 + for<'b> Sub<&'b Self, Output = <Self as VectorSpace>::PrincipalV>
42 {
43 /// Underlying scalar field
44 type Field: Num;
45
46 /// Principal form of the space; always equal to [`Space::Principal`], but with
47 /// more traits guaranteed.
48 ///
49 /// `PrincipalV` is only assumed to be `AXPY` for itself, as [`AXPY`]
50 /// uses [`Instance`] to apply all other variants and avoid problems
51 /// of choosing multiple implementations of the trait.
52 type PrincipalV: ClosedSpace
53 + AXPY<
54 Self::PrincipalV,
55 Field = Self::Field,
56 PrincipalV = Self::PrincipalV,
57 OwnedVariant = Self::PrincipalV,
58 Principal = Self::PrincipalV,
59 >;
60
61 /// Return a similar zero as `self`.
62 fn similar_origin(&self) -> Self::PrincipalV;
63 // {
64 // self.make_origin_generator().make_origin()
65 // }
66
67 /// Return a similar zero as `x`.
68 fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::PrincipalV {
69 x.eval(|xr| xr.similar_origin())
70 }
71 }
17 72
18 /// Efficient in-place summation. 73 /// Efficient in-place summation.
19 #[replace_float_literals(F::cast_from(literal))] 74 #[replace_float_literals(Self::Field::cast_from(literal))]
20 pub trait AXPY<F, X = Self> : Space + std::ops::MulAssign<F> 75 pub trait AXPY<X = Self>:
21 where 76 VectorSpace
22 F : Num, 77 + MulAssign<Self::Field>
23 X : Space, 78 + DivAssign<Self::Field>
24 { 79 + AddAssign<Self>
25 type Owned : AXPY<F, X>; 80 + AddAssign<Self::PrincipalV>
26 81 + SubAssign<Self>
82 + SubAssign<Self::PrincipalV>
83 + for<'b> AddAssign<&'b Self>
84 + for<'b> SubAssign<&'b Self>
85 where
86 X: Space,
87 {
27 /// Computes `y = βy + αx`, where `y` is `Self`. 88 /// Computes `y = βy + αx`, where `y` is `Self`.
28 fn axpy<I : Instance<X>>(&mut self, α : F, x : I, β : F); 89 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field);
29 90
30 /// Copies `x` to `self`. 91 /// Copies `x` to `self`.
31 fn copy_from<I : Instance<X>>(&mut self, x : I) { 92 fn copy_from<I: Instance<X>>(&mut self, x: I) {
32 self.axpy(1.0, x, 0.0) 93 self.axpy(1.0, x, 0.0)
33 } 94 }
34 95
35 /// Computes `y = αx`, where `y` is `Self`. 96 /// Computes `y = αx`, where `y` is `Self`.
36 fn scale_from<I : Instance<X>>(&mut self, α : F, x : I) { 97 fn scale_from<I: Instance<X>>(&mut self, α: Self::Field, x: I) {
37 self.axpy(α, x, 0.0) 98 self.axpy(α, x, 0.0)
38 } 99 }
39
40 /// Return a similar zero as `self`.
41 fn similar_origin(&self) -> Self::Owned;
42 100
43 /// Set self to zero. 101 /// Set self to zero.
44 fn set_zero(&mut self); 102 fn set_zero(&mut self);
45 } 103 }
46 104
105 pub trait ClosedVectorSpace: Instance<Self> + VectorSpace<PrincipalV = Self> {}
106 impl<X: Instance<X> + VectorSpace<PrincipalV = Self>> ClosedVectorSpace for X {}
107
47 /// Efficient in-place application for [`Linear`] operators. 108 /// Efficient in-place application for [`Linear`] operators.
48 #[replace_float_literals(F::cast_from(literal))] 109 #[replace_float_literals(F::cast_from(literal))]
49 pub trait GEMV<F : Num, X : Space, Y = <Self as Mapping<X>>::Codomain> : Linear<X> { 110 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> {
50 /// Computes `y = αAx + βy`, where `A` is `Self`. 111 /// Computes `y = αAx + βy`, where `A` is `Self`.
51 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F); 112 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F);
52 113
53 #[inline] 114 #[inline]
54 /// Computes `y = Ax`, where `A` is `Self` 115 /// Computes `y = Ax`, where `A` is `Self`
55 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 116 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
56 self.gemv(y, 1.0, x, 0.0) 117 self.gemv(y, 1.0, x, 0.0)
57 } 118 }
58 119
59 #[inline] 120 #[inline]
60 /// Computes `y += Ax`, where `A` is `Self` 121 /// Computes `y += Ax`, where `A` is `Self`
61 fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ 122 fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
62 self.gemv(y, 1.0, x, 1.0) 123 self.gemv(y, 1.0, x, 1.0)
63 } 124 }
64 } 125 }
65 126
66
67 /// Bounded linear operators 127 /// Bounded linear operators
68 pub trait BoundedLinear<X, XExp, CodExp, F = f64> : Linear<X> 128 pub trait BoundedLinear<X, XExp, CodExp, F = f64>: Linear<X>
69 where 129 where
70 F : Num, 130 F: Num,
71 X : Space + Norm<F, XExp>, 131 X: Space,
72 XExp : NormExponent, 132 XExp: NormExponent,
73 CodExp : NormExponent 133 CodExp: NormExponent,
74 { 134 {
75 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`. 135 /// A bound on the operator norm $\|A\|$ for the linear operator $A$=`self`.
76 /// This is not expected to be the norm, just any bound on it that can be 136 /// This is not expected to be the norm, just any bound on it that can be
77 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm 137 /// reasonably implemented. The [`NormExponent`] `xexp` indicates the norm
78 /// in `X`, and `codexp` in the codomain. 138 /// in `X`, and `codexp` in the codomain.
79 fn opnorm_bound(&self, xexp : XExp, codexp : CodExp) -> F; 139 ///
140 /// This may fail with an error if the bound is for some reason incalculable.
141 fn opnorm_bound(&self, xexp: XExp, codexp: CodExp) -> DynResult<F>;
80 } 142 }
81 143
82 // Linear operator application into mutable target. The [`AsRef`] bound 144 // Linear operator application into mutable target. The [`AsRef`] bound
83 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`; 145 // is used to guarantee compatibility with `Yʹ` and `Self::Codomain`;
84 // the former is assumed to be e.g. a view into the latter. 146 // the former is assumed to be e.g. a view into the latter.
88 self.apply(x) 150 self.apply(x)
89 } 151 }
90 }*/ 152 }*/
91 153
92 /// Trait for forming the adjoint operator of `Self`. 154 /// Trait for forming the adjoint operator of `Self`.
93 pub trait Adjointable<X, Yʹ> : Linear<X> 155 pub trait Adjointable<X, Yʹ>: Linear<X>
94 where 156 where
95 X : Space, 157 X: Space,
96 Yʹ : Space, 158 Yʹ: Space,
97 { 159 {
98 type AdjointCodomain : Space; 160 /// Codomain of the adjoint operator.
99 type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a; 161 type AdjointCodomain: ClosedSpace;
162 /// Type of the adjoint operator.
163 type Adjoint<'a>: Linear<Yʹ, Codomain = Self::AdjointCodomain>
164 where
165 Self: 'a;
100 166
101 /// Form the adjoint operator of `self`. 167 /// Form the adjoint operator of `self`.
102 fn adjoint(&self) -> Self::Adjoint<'_>; 168 fn adjoint(&self) -> Self::Adjoint<'_>;
169 }
170
171 /// Variant of [`Adjointable`] where the adjoint does not depend on a lifetime parameter.
172 /// This exists due to restrictions of Rust's type system: if `A :: Adjointable`, and we make
173 /// further restrictions on the adjoint operator, through, e.g.
174 /// ```
175 /// for<'a> A::Adjoint<'a> : GEMV<F, Z, Y>,
176 /// ```
177 /// Then `'static` lifetime is forced on `X`. Having `A::SimpleAdjoint` not depend on `'a`
178 /// avoids this, but makes it impossible for the adjoint to be just a light wrapper around the
179 /// original operator.
180 pub trait SimplyAdjointable<X, Yʹ>: Linear<X>
181 where
182 X: Space,
183 Yʹ: Space,
184 {
185 /// Codomain of the adjoint operator.
186 type AdjointCodomain: ClosedSpace;
187 /// Type of the adjoint operator.
188 type SimpleAdjoint: Linear<Yʹ, Codomain = Self::AdjointCodomain>;
189
190 /// Form the adjoint operator of `self`.
191 fn adjoint(&self) -> Self::SimpleAdjoint;
103 } 192 }
104 193
105 /// Trait for forming a preadjoint of an operator. 194 /// Trait for forming a preadjoint of an operator.
106 /// 195 ///
107 /// For an operator $A$ this is an operator $A\_\*$ 196 /// For an operator $A$ this is an operator $A\_\*$
110 /// domain of the adjointed operator. `Self::Preadjoint` should be 199 /// domain of the adjointed operator. `Self::Preadjoint` should be
111 /// [`Adjointable`]`<'a,Ypre,X>`. 200 /// [`Adjointable`]`<'a,Ypre,X>`.
112 /// We do not make additional restrictions on `Self::Preadjoint` (in particular, it 201 /// We do not make additional restrictions on `Self::Preadjoint` (in particular, it
113 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint 202 /// does not have to be adjointable) to allow `X` to be a subspace yet the preadjoint
114 /// have the full space as the codomain, etc. 203 /// have the full space as the codomain, etc.
115 pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> { 204 pub trait Preadjointable<X: Space, Ypre: Space = <Self as Mapping<X>>::Codomain>:
116 type PreadjointCodomain : Space; 205 Linear<X>
117 type Preadjoint<'a> : Linear< 206 {
118 Ypre, Codomain=Self::PreadjointCodomain 207 type PreadjointCodomain: ClosedSpace;
119 > where Self : 'a; 208 type Preadjoint<'a>: Linear<Ypre, Codomain = Self::PreadjointCodomain>
209 where
210 Self: 'a;
120 211
121 /// Form the adjoint operator of `self`. 212 /// Form the adjoint operator of `self`.
122 fn preadjoint(&self) -> Self::Preadjoint<'_>; 213 fn preadjoint(&self) -> Self::Preadjoint<'_>;
123 } 214 }
124 215
125 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
126 pub trait SimplyAdjointable<X : Space> : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
127 impl<'a,X : Space, T> SimplyAdjointable<X> for T
128 where T : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
129
130 /// The identity operator 216 /// The identity operator
131 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 217 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
132 pub struct IdOp<X> (PhantomData<X>); 218 pub struct IdOp<X>(PhantomData<X>);
133 219
134 impl<X> IdOp<X> { 220 impl<X> IdOp<X> {
135 pub fn new() -> IdOp<X> { IdOp(PhantomData) } 221 pub fn new() -> IdOp<X> {
136 } 222 IdOp(PhantomData)
137 223 }
138 impl<X : Clone + Space> Mapping<X> for IdOp<X> { 224 }
139 type Codomain = X; 225
140 226 impl<X: Space> Mapping<X> for IdOp<X> {
141 fn apply<I : Instance<X>>(&self, x : I) -> X { 227 type Codomain = X::Principal;
228
229 fn apply<I: Instance<X>>(&self, x: I) -> Self::Codomain {
142 x.own() 230 x.own()
143 } 231 }
144 } 232 }
145 233
146 impl<X : Clone + Space> Linear<X> for IdOp<X> 234 impl<X: Space> Linear<X> for IdOp<X> {}
147 { }
148 235
149 #[replace_float_literals(F::cast_from(literal))] 236 #[replace_float_literals(F::cast_from(literal))]
150 impl<F : Num, X, Y> GEMV<F, X, Y> for IdOp<X> 237 impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X>
151 where 238 where
152 Y : AXPY<F, X>, 239 Y: AXPY<X, Field = F>,
153 X : Clone + Space 240 X: Space,
154 { 241 {
155 // Computes `y = αAx + βy`, where `A` is `Self`. 242 // Computes `y = αAx + βy`, where `A` is `Self`.
156 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { 243 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
157 y.axpy(α, x, β) 244 y.axpy(α, x, β)
158 } 245 }
159 246
160 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 247 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
161 y.copy_from(x); 248 y.copy_from(x);
162 } 249 }
163 } 250 }
164 251
165 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X> 252 impl<F, X, E> BoundedLinear<X, E, E, F> for IdOp<X>
166 where 253 where
167 X : Space + Clone + Norm<F, E>, 254 X: Space + Clone,
168 F : Num, 255 F: Num,
169 E : NormExponent 256 E: NormExponent,
170 { 257 {
171 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE } 258 fn opnorm_bound(&self, _xexp: E, _codexp: E) -> DynResult<F> {
172 } 259 Ok(F::ONE)
173 260 }
174 impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> { 261 }
175 type AdjointCodomain=X; 262
176 type Adjoint<'a> = IdOp<X> where X : 'a; 263 impl<X: Clone + Space> Adjointable<X, X::Principal> for IdOp<X> {
177 264 type AdjointCodomain = X::Principal;
178 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } 265 type Adjoint<'a>
179 } 266 = IdOp<X::Principal>
180 267 where
181 impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> { 268 X: 'a;
182 type PreadjointCodomain=X; 269
183 type Preadjoint<'a> = IdOp<X> where X : 'a; 270 fn adjoint(&self) -> Self::Adjoint<'_> {
184 271 IdOp::new()
185 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() } 272 }
186 } 273 }
187 274
188 275 impl<X: Clone + Space> SimplyAdjointable<X, X::Principal> for IdOp<X> {
189 /// The zero operator 276 type AdjointCodomain = X::Principal;
190 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 277 type SimpleAdjoint = IdOp<X::Principal>;
191 pub struct ZeroOp<'a, X, XD, Y, F> { 278
192 zero : &'a Y, // TODO: don't pass this in `new`; maybe not even store. 279 fn adjoint(&self) -> Self::SimpleAdjoint {
193 dual_or_predual_zero : XD, 280 IdOp::new()
194 _phantoms : PhantomData<(X, Y, F)>, 281 }
195 } 282 }
196 283
197 // TODO: Need to make Zero in Instance. 284 impl<X: Clone + Space> Preadjointable<X, X::Principal> for IdOp<X> {
198 285 type PreadjointCodomain = X::Principal;
199 impl<'a, F : Num, X : Space, XD, Y : Space + Clone> ZeroOp<'a, X, XD, Y, F> { 286 type Preadjoint<'a>
200 pub fn new(zero : &'a Y, dual_or_predual_zero : XD) -> Self { 287 = IdOp<X::Principal>
201 ZeroOp{ zero, dual_or_predual_zero, _phantoms : PhantomData } 288 where
202 } 289 X: 'a;
203 } 290
204 291 fn preadjoint(&self) -> Self::Preadjoint<'_> {
205 impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Mapping<X> for ZeroOp<'a, X, XD, Y, F> { 292 IdOp::new()
206 type Codomain = Y; 293 }
207 294 }
208 fn apply<I : Instance<X>>(&self, _x : I) -> Y { 295
209 self.zero.clone() 296 /// The zero operator from a space to itself
210 } 297 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
211 } 298 pub struct SimpleZeroOp;
212 299
213 impl<'a, F : Num, X : Space, XD, Y : AXPY<F> + Clone> Linear<X> for ZeroOp<'a, X, XD, Y, F> 300 impl<X: VectorSpace> Mapping<X> for SimpleZeroOp {
214 { } 301 type Codomain = X::PrincipalV;
302
303 fn apply<I: Instance<X>>(&self, x: I) -> X::PrincipalV {
304 X::similar_origin_inst(x)
305 }
306 }
307
308 impl<X: VectorSpace> Linear<X> for SimpleZeroOp {}
215 309
216 #[replace_float_literals(F::cast_from(literal))] 310 #[replace_float_literals(F::cast_from(literal))]
217 impl<'a, F, X, XD, Y> GEMV<F, X, Y> for ZeroOp<'a, X, XD, Y, F> 311 impl<X, Y, F> GEMV<F, X, Y> for SimpleZeroOp
218 where 312 where
219 F : Num, 313 F: Num,
220 Y : AXPY<F, Y> + Clone, 314 Y: AXPY<Field = F>,
221 X : Space 315 X: VectorSpace<Field = F> + Instance<X>,
222 { 316 {
223 // Computes `y = αAx + βy`, where `A` is `Self`. 317 // Computes `y = αAx + βy`, where `A` is `Self`.
224 fn gemv<I : Instance<X>>(&self, y : &mut Y, _α : F, _x : I, β : F) { 318 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
225 *y *= β; 319 *y *= β;
226 } 320 }
227 321
228 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, _x : I){ 322 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
229 y.set_zero(); 323 y.set_zero();
230 } 324 }
231 } 325 }
232 326
233 impl<'a, F, X, XD, Y, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<'a, X, XD, Y, F> 327 impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp
234 where 328 where
235 X : Space + Norm<F, E1>, 329 F: Num,
236 Y : AXPY<F> + Clone + Norm<F, E2>, 330 X: VectorSpace<Field = F>,
237 F : Num, 331 E1: NormExponent,
238 E1 : NormExponent, 332 E2: NormExponent,
239 E2 : NormExponent, 333 {
240 { 334 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
241 fn opnorm_bound(&self, _xexp : E1, _codexp : E2) -> F { F::ZERO } 335 Ok(F::ZERO)
242 } 336 }
243 337 }
244 impl<'a, F : Num, X, XD, Y, Yprime : Space> Adjointable<X, Yprime> for ZeroOp<'a, X, XD, Y, F> 338
245 where 339 impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp
246 X : Space, 340 where
247 Y : AXPY<F> + Clone + 'static, 341 F: Num,
248 XD : AXPY<F> + Clone + 'static, 342 X: VectorSpace<Field = F> + HasDual<F>,
249 { 343 X::DualSpace: ClosedVectorSpace,
250 type AdjointCodomain = XD; 344 {
251 type Adjoint<'b> = ZeroOp<'b, Yprime, (), XD, F> where Self : 'b; 345 type AdjointCodomain = X::DualSpace;
252 // () means not (pre)adjointable. 346 type Adjoint<'b>
347 = SimpleZeroOp
348 where
349 Self: 'b;
350 // () means not (pre)adjointable.
253 351
254 fn adjoint(&self) -> Self::Adjoint<'_> { 352 fn adjoint(&self) -> Self::Adjoint<'_> {
255 ZeroOp::new(&self.dual_or_predual_zero, ()) 353 SimpleZeroOp
256 } 354 }
257 } 355 }
258 356
259 impl<'a, F, X, XD, Y, Ypre> Preadjointable<X, Ypre> for ZeroOp<'a, X, XD, Y, F> 357 pub trait OriginGenerator<Y: VectorSpace> {
260 where 358 type Ref<'b>: OriginGenerator<Y>
261 F : Num, 359 where
262 X : Space, 360 Self: 'b;
263 Y : AXPY<F> + Clone, 361
264 Ypre : Space, 362 fn origin(&self) -> Y::PrincipalV;
265 XD : AXPY<F> + Clone + 'static, 363 fn as_ref(&self) -> Self::Ref<'_>;
266 { 364 }
267 type PreadjointCodomain = XD; 365
268 type Preadjoint<'b> = ZeroOp<'b, Ypre, (), XD, F> where Self : 'b; 366 #[derive(Copy, Clone, Debug)]
269 // () means not (pre)adjointable. 367 pub struct StaticEuclideanOriginGenerator;
270 368
271 fn preadjoint(&self) -> Self::Preadjoint<'_> { 369 impl<Y: StaticEuclidean<F, Field = F>, F: Float> OriginGenerator<Y>
272 ZeroOp::new(&self.dual_or_predual_zero, ()) 370 for StaticEuclideanOriginGenerator
371 {
372 type Ref<'b>
373 = Self
374 where
375 Self: 'b;
376
377 #[inline]
378 fn origin(&self) -> Y::PrincipalV {
379 return Y::origin();
380 }
381
382 #[inline]
383 fn as_ref(&self) -> Self::Ref<'_> {
384 *self
385 }
386 }
387
388 impl<Y: VectorSpace> OriginGenerator<Y> for Y {
389 type Ref<'b>
390 = &'b Y
391 where
392 Self: 'b;
393
394 #[inline]
395 fn origin(&self) -> Y::PrincipalV {
396 return self.similar_origin();
397 }
398
399 #[inline]
400 fn as_ref(&self) -> Self::Ref<'_> {
401 self
402 }
403 }
404
405 impl<'b, Y: VectorSpace> OriginGenerator<Y> for &'b Y {
406 type Ref<'c>
407 = Self
408 where
409 Self: 'c;
410
411 #[inline]
412 fn origin(&self) -> Y::PrincipalV {
413 return self.similar_origin();
414 }
415
416 #[inline]
417 fn as_ref(&self) -> Self::Ref<'_> {
418 self
419 }
420 }
421
422 /// A zero operator that can be eitherh dualised or predualised (once).
423 /// This is achieved by storing an oppropriate zero.
424 pub struct ZeroOp<X, Y: VectorSpace<Field = F>, OY: OriginGenerator<Y>, O, F: Float = f64> {
425 codomain_origin_generator: OY,
426 other_origin_generator: O,
427 _phantoms: PhantomData<(X, Y, F)>,
428 }
429
430 impl<X, Y, OY, F> ZeroOp<X, Y, OY, (), F>
431 where
432 OY: OriginGenerator<Y>,
433 X: VectorSpace<Field = F>,
434 Y: VectorSpace<Field = F>,
435 F: Float,
436 {
437 pub fn new(y_og: OY) -> Self {
438 ZeroOp {
439 codomain_origin_generator: y_og,
440 other_origin_generator: (),
441 _phantoms: PhantomData,
442 }
443 }
444 }
445
446 impl<X, Y, OY, OXprime, Xprime, Yprime, F> ZeroOp<X, Y, OY, OXprime, F>
447 where
448 OY: OriginGenerator<Y>,
449 OXprime: OriginGenerator<Xprime>,
450 X: HasDual<F, DualSpace = Xprime>,
451 Y: HasDual<F, DualSpace = Yprime>,
452 F: Float,
453 Xprime: VectorSpace<Field = F, PrincipalV = Xprime>,
454 Xprime::PrincipalV: AXPY<Field = F>,
455 Yprime: Space + Instance<Yprime>,
456 {
457 pub fn new_dualisable(y_og: OY, xprime_og: OXprime) -> Self {
458 ZeroOp {
459 codomain_origin_generator: y_og,
460 other_origin_generator: xprime_og,
461 _phantoms: PhantomData,
462 }
463 }
464 }
465
466 impl<X, Y, O, OY, F> Mapping<X> for ZeroOp<X, Y, OY, O, F>
467 where
468 X: Space,
469 Y: VectorSpace<Field = F>,
470 F: Float,
471 OY: OriginGenerator<Y>,
472 {
473 type Codomain = Y::PrincipalV;
474
475 fn apply<I: Instance<X>>(&self, _x: I) -> Y::PrincipalV {
476 self.codomain_origin_generator.origin()
477 }
478 }
479
480 impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F>
481 where
482 X: Space,
483 Y: VectorSpace<Field = F>,
484 F: Float,
485 OY: OriginGenerator<Y>,
486 {
487 }
488
489 #[replace_float_literals(F::cast_from(literal))]
490 impl<X, Y, OY, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, OY, O, F>
491 where
492 X: Space,
493 Y: AXPY<Field = F, PrincipalV = Y>,
494 F: Float,
495 OY: OriginGenerator<Y>,
496 {
497 // Computes `y = αAx + βy`, where `A` is `Self`.
498 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
499 *y *= β;
500 }
501
502 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
503 y.set_zero();
504 }
505 }
506
507 impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F>
508 where
509 X: Space + Instance<X>,
510 Y: VectorSpace<Field = F>,
511 Y::PrincipalV: Clone,
512 F: Float,
513 E1: NormExponent,
514 E2: NormExponent,
515 OY: OriginGenerator<Y>,
516 {
517 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
518 Ok(F::ZERO)
519 }
520 }
521
522 impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> Adjointable<X, Yprime>
523 for ZeroOp<X, Y, OY, OXprime, F>
524 where
525 X: HasDual<F, DualSpace = Xprime>,
526 Y: HasDual<F, DualSpace = Yprime>,
527 F: Float,
528 Xprime: ClosedVectorSpace<Field = F>,
529 //Xprime::Owned: AXPY<Field = F>,
530 Yprime: ClosedSpace,
531 OY: OriginGenerator<Y>,
532 OXprime: OriginGenerator<X::DualSpace>,
533 {
534 type AdjointCodomain = Xprime;
535 type Adjoint<'c>
536 = ZeroOp<Yprime, Xprime, OXprime::Ref<'c>, (), F>
537 where
538 Self: 'c;
539 // () means not (pre)adjointable.
540
541 fn adjoint(&self) -> Self::Adjoint<'_> {
542 ZeroOp {
543 codomain_origin_generator: self.other_origin_generator.as_ref(),
544 other_origin_generator: (),
545 _phantoms: PhantomData,
546 }
547 }
548 }
549
550 impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> SimplyAdjointable<X, Yprime>
551 for ZeroOp<X, Y, OY, OXprime, F>
552 where
553 X: HasDual<F, DualSpace = Xprime>,
554 Y: HasDual<F, DualSpace = Yprime>,
555 F: Float,
556 Xprime: ClosedVectorSpace<Field = F>,
557 //Xprime::Owned: AXPY<Field = F>,
558 Yprime: ClosedSpace,
559 OY: OriginGenerator<Y>,
560 OXprime: OriginGenerator<X::DualSpace> + Clone,
561 {
562 type AdjointCodomain = Xprime;
563 type SimpleAdjoint = ZeroOp<Yprime, Xprime, OXprime, (), F>;
564 // () means not (pre)adjointable.
565
566 fn adjoint(&self) -> Self::SimpleAdjoint {
567 ZeroOp {
568 codomain_origin_generator: self.other_origin_generator.clone(),
569 other_origin_generator: (),
570 _phantoms: PhantomData,
571 }
273 } 572 }
274 } 573 }
275 574
276 impl<S, T, E, X> Linear<X> for Composition<S, T, E> 575 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
277 where 576 where
278 X : Space, 577 X: Space,
279 T : Linear<X>, 578 T: Linear<X>,
280 S : Linear<T::Codomain> 579 S: Linear<T::Codomain>,
281 { } 580 {
581 }
282 582
283 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E> 583 impl<F, S, T, E, X, Y> GEMV<F, X, Y> for Composition<S, T, E>
284 where 584 where
285 F : Num, 585 F: Num,
286 X : Space, 586 X: Space,
287 T : Linear<X>, 587 T: Linear<X>,
288 S : GEMV<F, T::Codomain, Y>, 588 S: GEMV<F, T::Codomain, Y>,
289 { 589 {
290 fn gemv<I : Instance<X>>(&self, y : &mut Y, α : F, x : I, β : F) { 590 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
291 self.outer.gemv(y, α, self.inner.apply(x), β) 591 self.outer.gemv(y, α, self.inner.apply(x), β)
292 } 592 }
293 593
294 /// Computes `y = Ax`, where `A` is `Self` 594 /// Computes `y = Ax`, where `A` is `Self`
295 fn apply_mut<I : Instance<X>>(&self, y : &mut Y, x : I){ 595 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, x: I) {
296 self.outer.apply_mut(y, self.inner.apply(x)) 596 self.outer.apply_mut(y, self.inner.apply(x))
297 } 597 }
298 598
299 /// Computes `y += Ax`, where `A` is `Self` 599 /// Computes `y += Ax`, where `A` is `Self`
300 fn apply_add<I : Instance<X>>(&self, y : &mut Y, x : I){ 600 fn apply_add<I: Instance<X>>(&self, y: &mut Y, x: I) {
301 self.outer.apply_add(y, self.inner.apply(x)) 601 self.outer.apply_add(y, self.inner.apply(x))
302 } 602 }
303 } 603 }
304 604
305 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp> 605 impl<F, S, T, X, Z, Xexp, Yexp, Zexp> BoundedLinear<X, Xexp, Yexp, F> for Composition<S, T, Zexp>
306 where 606 where
307 F : Num, 607 F: Num,
308 X : Space + Norm<F, Xexp>, 608 X: Space,
309 Z : Space + Norm<F, Zexp>, 609 Z: Space,
310 Xexp : NormExponent, 610 Xexp: NormExponent,
311 Yexp : NormExponent, 611 Yexp: NormExponent,
312 Zexp : NormExponent, 612 Zexp: NormExponent,
313 T : BoundedLinear<X, Xexp, Zexp, F, Codomain=Z>, 613 T: BoundedLinear<X, Xexp, Zexp, F, Codomain = Z>,
314 S : BoundedLinear<Z, Zexp, Yexp, F>, 614 S: BoundedLinear<Z, Zexp, Yexp, F>,
315 { 615 {
316 fn opnorm_bound(&self, xexp : Xexp, yexp : Yexp) -> F { 616 fn opnorm_bound(&self, xexp: Xexp, yexp: Yexp) -> DynResult<F> {
317 let zexp = self.intermediate_norm_exponent; 617 let zexp = self.intermediate_norm_exponent;
318 self.outer.opnorm_bound(zexp, yexp) * self.inner.opnorm_bound(xexp, zexp) 618 Ok(self.outer.opnorm_bound(zexp, yexp)? * self.inner.opnorm_bound(xexp, zexp)?)
319 } 619 }
320 } 620 }
321 621
322 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$. 622 /// “Row operator” $(S, T)$; $(S, T)(x, y)=Sx + Ty$.
623 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
323 pub struct RowOp<S, T>(pub S, pub T); 624 pub struct RowOp<S, T>(pub S, pub T);
324 625
325 use std::ops::Add;
326
327 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T> 626 impl<A, B, S, T> Mapping<Pair<A, B>> for RowOp<S, T>
328 where 627 where
329 A : Space, 628 A: Space,
330 B : Space, 629 B: Space,
331 S : Mapping<A>, 630 S: Mapping<A>,
332 T : Mapping<B>, 631 T: Mapping<B>,
333 S::Codomain : Add<T::Codomain>, 632 S::Codomain: Add<T::Codomain>,
334 <S::Codomain as Add<T::Codomain>>::Output : Space, 633 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
335
336 { 634 {
337 type Codomain = <S::Codomain as Add<T::Codomain>>::Output; 635 type Codomain = <S::Codomain as Add<T::Codomain>>::Output;
338 636
339 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { 637 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
340 let Pair(a, b) = x.decompose(); 638 x.eval_decompose(|Pair(a, b)| self.0.apply(a) + self.1.apply(b))
341 self.0.apply(a) + self.1.apply(b)
342 } 639 }
343 } 640 }
344 641
345 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T> 642 impl<A, B, S, T> Linear<Pair<A, B>> for RowOp<S, T>
346 where 643 where
347 A : Space, 644 A: Space,
348 B : Space, 645 B: Space,
349 S : Linear<A>, 646 S: Linear<A>,
350 T : Linear<B>, 647 T: Linear<B>,
351 S::Codomain : Add<T::Codomain>, 648 S::Codomain: Add<T::Codomain>,
352 <S::Codomain as Add<T::Codomain>>::Output : Space, 649 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
353 { } 650 {
354 651 }
355 652
356 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T> 653 impl<'b, F, S, T, Y, U, V> GEMV<F, Pair<U, V>, Y> for RowOp<S, T>
357 where 654 where
358 U : Space, 655 U: Space,
359 V : Space, 656 V: Space,
360 S : GEMV<F, U, Y>, 657 S: GEMV<F, U, Y>,
361 T : GEMV<F, V, Y>, 658 T: GEMV<F, V, Y>,
362 F : Num, 659 F: Num,
363 Self : Linear<Pair<U, V>, Codomain=Y> 660 Self: Linear<Pair<U, V>, Codomain = Y>,
364 { 661 {
365 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Y, α : F, x : I, β : F) { 662 fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Y, α: F, x: I, β: F) {
366 let Pair(u, v) = x.decompose(); 663 x.eval_decompose(|Pair(u, v)| {
367 self.0.gemv(y, α, u, β); 664 self.0.gemv(y, α, u, β);
368 self.1.gemv(y, α, v, F::ONE); 665 self.1.gemv(y, α, v, F::ONE);
369 } 666 })
370 667 }
371 fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) { 668
372 let Pair(u, v) = x.decompose(); 669 fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
373 self.0.apply_mut(y, u); 670 x.eval_decompose(|Pair(u, v)| {
374 self.1.apply_add(y, v); 671 self.0.apply_mut(y, u);
672 self.1.apply_add(y, v);
673 })
375 } 674 }
376 675
377 /// Computes `y += Ax`, where `A` is `Self` 676 /// Computes `y += Ax`, where `A` is `Self`
378 fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Y, x : I) { 677 fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Y, x: I) {
379 let Pair(u, v) = x.decompose(); 678 x.eval_decompose(|Pair(u, v)| {
380 self.0.apply_add(y, u); 679 self.0.apply_add(y, u);
381 self.1.apply_add(y, v); 680 self.1.apply_add(y, v);
681 })
382 } 682 }
383 } 683 }
384 684
385 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$. 685 /// “Column operator” $(S; T)$; $(S; T)x=(Sx, Tx)$.
686 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
386 pub struct ColOp<S, T>(pub S, pub T); 687 pub struct ColOp<S, T>(pub S, pub T);
387 688
388 impl<A, S, T> Mapping<A> for ColOp<S, T> 689 impl<A, S, T> Mapping<A> for ColOp<S, T>
389 where 690 where
390 A : Space, 691 A: Space,
391 S : Mapping<A>, 692 S: Mapping<A>,
392 T : Mapping<A>, 693 T: Mapping<A>,
393 { 694 {
394 type Codomain = Pair<S::Codomain, T::Codomain>; 695 type Codomain = Pair<S::Codomain, T::Codomain>;
395 696
396 fn apply<I : Instance<A>>(&self, a : I) -> Self::Codomain { 697 fn apply<I: Instance<A>>(&self, a: I) -> Self::Codomain {
397 Pair(self.0.apply(a.ref_instance()), self.1.apply(a)) 698 Pair(a.eval_ref(|r| self.0.apply(r)), self.1.apply(a))
398 } 699 }
399 } 700 }
400 701
401 impl<A, S, T> Linear<A> for ColOp<S, T> 702 impl<A, S, T> Linear<A> for ColOp<S, T>
402 where 703 where
403 A : Space, 704 A: Space,
404 S : Mapping<A>, 705 S: Mapping<A>,
405 T : Mapping<A>, 706 T: Mapping<A>,
406 { } 707 {
708 }
407 709
408 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T> 710 impl<F, S, T, A, B, X> GEMV<F, X, Pair<A, B>> for ColOp<S, T>
409 where 711 where
410 X : Space, 712 X: Space,
411 S : GEMV<F, X, A>, 713 S: GEMV<F, X, A>,
412 T : GEMV<F, X, B>, 714 T: GEMV<F, X, B>,
413 F : Num, 715 F: Num,
414 Self : Linear<X, Codomain=Pair<A, B>> 716 Self: Linear<X, Codomain = Pair<A, B>>,
415 { 717 {
416 fn gemv<I : Instance<X>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { 718 fn gemv<I: Instance<X>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
417 self.0.gemv(&mut y.0, α, x.ref_instance(), β); 719 x.eval_ref(|r| self.0.gemv(&mut y.0, α, r, β));
418 self.1.gemv(&mut y.1, α, x, β); 720 self.1.gemv(&mut y.1, α, x, β);
419 } 721 }
420 722
421 fn apply_mut<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){ 723 fn apply_mut<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
422 self.0.apply_mut(&mut y.0, x.ref_instance()); 724 x.eval_ref(|r| self.0.apply_mut(&mut y.0, r));
423 self.1.apply_mut(&mut y.1, x); 725 self.1.apply_mut(&mut y.1, x);
424 } 726 }
425 727
426 /// Computes `y += Ax`, where `A` is `Self` 728 /// Computes `y += Ax`, where `A` is `Self`
427 fn apply_add<I : Instance<X>>(&self, y : &mut Pair<A, B>, x : I){ 729 fn apply_add<I: Instance<X>>(&self, y: &mut Pair<A, B>, x: I) {
428 self.0.apply_add(&mut y.0, x.ref_instance()); 730 x.eval_ref(|r| self.0.apply_add(&mut y.0, r));
429 self.1.apply_add(&mut y.1, x); 731 self.1.apply_add(&mut y.1, x);
430 } 732 }
431 } 733 }
432 734
433 735 impl<A, B, Yʹ, S, T> Adjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
434 impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T> 736 where
435 where 737 A: Space,
436 A : Space, 738 B: Space,
437 B : Space, 739 Yʹ: Space,
438 Yʹ : Space, 740 S: Adjointable<A, Yʹ>,
439 S : Adjointable<A, Yʹ>, 741 T: Adjointable<B, Yʹ>,
440 T : Adjointable<B, Yʹ>, 742 Self: Linear<Pair<A, B>>,
441 Self : Linear<Pair<A, B>>,
442 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 743 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
443 // Yʹ, 744 // Yʹ,
444 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain> 745 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
445 // >, 746 // >,
446 { 747 {
447 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>; 748 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
448 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 749 type Adjoint<'a>
750 = ColOp<S::Adjoint<'a>, T::Adjoint<'a>>
751 where
752 Self: 'a;
449 753
450 fn adjoint(&self) -> Self::Adjoint<'_> { 754 fn adjoint(&self) -> Self::Adjoint<'_> {
451 ColOp(self.0.adjoint(), self.1.adjoint()) 755 ColOp(self.0.adjoint(), self.1.adjoint())
452 } 756 }
453 } 757 }
454 758
455 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A,B>, Yʹ> for RowOp<S, T> 759 impl<A, B, Yʹ, S, T> SimplyAdjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
456 where 760 where
457 A : Space, 761 A: Space,
458 B : Space, 762 B: Space,
459 Yʹ : Space, 763 Yʹ: Space,
460 S : Preadjointable<A, Yʹ>, 764 S: SimplyAdjointable<A, Yʹ>,
461 T : Preadjointable<B, Yʹ>, 765 T: SimplyAdjointable<B, Yʹ>,
462 Self : Linear<Pair<A, B>>, 766 Self: Linear<Pair<A, B>>,
463 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< 767 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
464 Yʹ, Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>, 768 // Yʹ,
465 >, 769 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
770 // >,
771 {
772 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
773 type SimpleAdjoint = ColOp<S::SimpleAdjoint, T::SimpleAdjoint>;
774
775 fn adjoint(&self) -> Self::SimpleAdjoint {
776 ColOp(self.0.adjoint(), self.1.adjoint())
777 }
778 }
779
780 impl<A, B, Yʹ, S, T> Preadjointable<Pair<A, B>, Yʹ> for RowOp<S, T>
781 where
782 A: Space,
783 B: Space,
784 Yʹ: Space,
785 S: Preadjointable<A, Yʹ>,
786 T: Preadjointable<B, Yʹ>,
787 Self: Linear<Pair<A, B>>,
788 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>:
789 Linear<Yʹ, Codomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>>,
466 { 790 {
467 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>; 791 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
468 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 792 type Preadjoint<'a>
793 = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
794 where
795 Self: 'a;
469 796
470 fn preadjoint(&self) -> Self::Preadjoint<'_> { 797 fn preadjoint(&self) -> Self::Preadjoint<'_> {
471 ColOp(self.0.preadjoint(), self.1.preadjoint()) 798 ColOp(self.0.preadjoint(), self.1.preadjoint())
472 } 799 }
473 } 800 }
474 801
475 802 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
476 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 803 where
477 where 804 A: Space,
478 A : Space, 805 Xʹ: Space,
479 Xʹ : Space, 806 Yʹ: Space,
480 Yʹ : Space, 807 R: ClosedSpace + ClosedAdd,
481 R : Space + ClosedAdd, 808 S: Adjointable<A, Xʹ, AdjointCodomain = R>,
482 S : Adjointable<A, Xʹ, AdjointCodomain = R>, 809 T: Adjointable<A, Yʹ, AdjointCodomain = R>,
483 T : Adjointable<A, Yʹ, AdjointCodomain = R>, 810 Self: Linear<A>,
484 Self : Linear<A>,
485 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 811 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
486 // Pair<Xʹ,Yʹ>, 812 // Pair<Xʹ,Yʹ>,
487 // Codomain=R, 813 // Codomain=R,
488 // >, 814 // >,
489 { 815 {
490 type AdjointCodomain = R; 816 type AdjointCodomain = R;
491 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 817 type Adjoint<'a>
818 = RowOp<S::Adjoint<'a>, T::Adjoint<'a>>
819 where
820 Self: 'a;
492 821
493 fn adjoint(&self) -> Self::Adjoint<'_> { 822 fn adjoint(&self) -> Self::Adjoint<'_> {
494 RowOp(self.0.adjoint(), self.1.adjoint()) 823 RowOp(self.0.adjoint(), self.1.adjoint())
495 } 824 }
496 } 825 }
497 826
498 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 827 impl<A, Xʹ, Yʹ, R, S, T> SimplyAdjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
499 where 828 where
500 A : Space, 829 A: Space,
501 Xʹ : Space, 830 Xʹ: Space,
502 Yʹ : Space, 831 Yʹ: Space,
503 R : Space + ClosedAdd, 832 R: ClosedSpace + ClosedAdd,
504 S : Preadjointable<A, Xʹ, PreadjointCodomain = R>, 833 S: SimplyAdjointable<A, Xʹ, AdjointCodomain = R>,
505 T : Preadjointable<A, Yʹ, PreadjointCodomain = R>, 834 T: SimplyAdjointable<A, Yʹ, AdjointCodomain = R>,
506 Self : Linear<A>, 835 Self: Linear<A>,
507 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< 836 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
508 Pair<Xʹ,Yʹ>, Codomain = R, 837 // Pair<Xʹ,Yʹ>,
509 >, 838 // Codomain=R,
839 // >,
840 {
841 type AdjointCodomain = R;
842 type SimpleAdjoint = RowOp<S::SimpleAdjoint, T::SimpleAdjoint>;
843
844 fn adjoint(&self) -> Self::SimpleAdjoint {
845 RowOp(self.0.adjoint(), self.1.adjoint())
846 }
847 }
848
849 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A, Pair<Xʹ, Yʹ>> for ColOp<S, T>
850 where
851 A: Space,
852 Xʹ: Space,
853 Yʹ: Space,
854 R: ClosedSpace + ClosedAdd,
855 S: Preadjointable<A, Xʹ, PreadjointCodomain = R>,
856 T: Preadjointable<A, Yʹ, PreadjointCodomain = R>,
857 Self: Linear<A>,
858 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
510 { 859 {
511 type PreadjointCodomain = R; 860 type PreadjointCodomain = R;
512 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 861 type Preadjoint<'a>
862 = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
863 where
864 Self: 'a;
513 865
514 fn preadjoint(&self) -> Self::Preadjoint<'_> { 866 fn preadjoint(&self) -> Self::Preadjoint<'_> {
515 RowOp(self.0.preadjoint(), self.1.preadjoint()) 867 RowOp(self.0.preadjoint(), self.1.preadjoint())
516 } 868 }
517 } 869 }
518 870
519 /// Diagonal operator 871 /// Diagonal operator
872 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
520 pub struct DiagOp<S, T>(pub S, pub T); 873 pub struct DiagOp<S, T>(pub S, pub T);
521 874
522 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T> 875 impl<A, B, S, T> Mapping<Pair<A, B>> for DiagOp<S, T>
523 where 876 where
524 A : Space, 877 A: Space,
525 B : Space, 878 B: Space,
526 S : Mapping<A>, 879 S: Mapping<A>,
527 T : Mapping<B>, 880 T: Mapping<B>,
528 { 881 {
529 type Codomain = Pair<S::Codomain, T::Codomain>; 882 type Codomain = Pair<S::Codomain, T::Codomain>;
530 883
531 fn apply<I : Instance<Pair<A, B>>>(&self, x : I) -> Self::Codomain { 884 fn apply<I: Instance<Pair<A, B>>>(&self, x: I) -> Self::Codomain {
532 let Pair(a, b) = x.decompose(); 885 x.eval_decompose(|Pair(a, b)| Pair(self.0.apply(a), self.1.apply(b)))
533 Pair(self.0.apply(a), self.1.apply(b))
534 } 886 }
535 } 887 }
536 888
537 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T> 889 impl<A, B, S, T> Linear<Pair<A, B>> for DiagOp<S, T>
538 where 890 where
539 A : Space, 891 A: Space,
540 B : Space, 892 B: Space,
541 S : Linear<A>, 893 S: Linear<A>,
542 T : Linear<B>, 894 T: Linear<B>,
543 { } 895 {
896 }
544 897
545 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T> 898 impl<F, S, T, A, B, U, V> GEMV<F, Pair<U, V>, Pair<A, B>> for DiagOp<S, T>
546 where 899 where
547 A : Space, 900 A: Space,
548 B : Space, 901 B: Space,
549 U : Space, 902 U: Space,
550 V : Space, 903 V: Space,
551 S : GEMV<F, U, A>, 904 S: GEMV<F, U, A>,
552 T : GEMV<F, V, B>, 905 T: GEMV<F, V, B>,
553 F : Num, 906 F: Num,
554 Self : Linear<Pair<U, V>, Codomain=Pair<A, B>>, 907 Self: Linear<Pair<U, V>, Codomain = Pair<A, B>>,
555 { 908 {
556 fn gemv<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, α : F, x : I, β : F) { 909 fn gemv<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, α: F, x: I, β: F) {
557 let Pair(u, v) = x.decompose(); 910 x.eval_decompose(|Pair(u, v)| {
558 self.0.gemv(&mut y.0, α, u, β); 911 self.0.gemv(&mut y.0, α, u, β);
559 self.1.gemv(&mut y.1, α, v, β); 912 self.1.gemv(&mut y.1, α, v, β);
560 } 913 })
561 914 }
562 fn apply_mut<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){ 915
563 let Pair(u, v) = x.decompose(); 916 fn apply_mut<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
564 self.0.apply_mut(&mut y.0, u); 917 x.eval_decompose(|Pair(u, v)| {
565 self.1.apply_mut(&mut y.1, v); 918 self.0.apply_mut(&mut y.0, u);
919 self.1.apply_mut(&mut y.1, v);
920 })
566 } 921 }
567 922
568 /// Computes `y += Ax`, where `A` is `Self` 923 /// Computes `y += Ax`, where `A` is `Self`
569 fn apply_add<I : Instance<Pair<U, V>>>(&self, y : &mut Pair<A, B>, x : I){ 924 fn apply_add<I: Instance<Pair<U, V>>>(&self, y: &mut Pair<A, B>, x: I) {
570 let Pair(u, v) = x.decompose(); 925 x.eval_decompose(|Pair(u, v)| {
571 self.0.apply_add(&mut y.0, u); 926 self.0.apply_add(&mut y.0, u);
572 self.1.apply_add(&mut y.1, v); 927 self.1.apply_add(&mut y.1, v);
573 } 928 })
574 } 929 }
575 930 }
576 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 931
577 where 932 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
578 A : Space, 933 where
579 B : Space, 934 A: Space,
935 B: Space,
580 Xʹ: Space, 936 Xʹ: Space,
581 Yʹ : Space, 937 Yʹ: Space,
582 R : Space, 938 R: ClosedSpace,
583 S : Adjointable<A, Xʹ>, 939 S: Adjointable<A, Xʹ>,
584 T : Adjointable<B, Yʹ>, 940 T: Adjointable<B, Yʹ>,
585 Self : Linear<Pair<A, B>>, 941 Self: Linear<Pair<A, B>>,
586 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 942 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
587 Pair<Xʹ,Yʹ>, Codomain=R,
588 >,
589 { 943 {
590 type AdjointCodomain = R; 944 type AdjointCodomain = R;
591 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 945 type Adjoint<'a>
946 = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>>
947 where
948 Self: 'a;
592 949
593 fn adjoint(&self) -> Self::Adjoint<'_> { 950 fn adjoint(&self) -> Self::Adjoint<'_> {
594 DiagOp(self.0.adjoint(), self.1.adjoint()) 951 DiagOp(self.0.adjoint(), self.1.adjoint())
595 } 952 }
596 } 953 }
597 954
598 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 955 impl<A, B, Xʹ, Yʹ, R, S, T> SimplyAdjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
599 where 956 where
600 A : Space, 957 A: Space,
601 B : Space, 958 B: Space,
602 Xʹ: Space, 959 Xʹ: Space,
603 Yʹ : Space, 960 Yʹ: Space,
604 R : Space, 961 R: ClosedSpace,
605 S : Preadjointable<A, Xʹ>, 962 S: SimplyAdjointable<A, Xʹ>,
606 T : Preadjointable<B, Yʹ>, 963 T: SimplyAdjointable<B, Yʹ>,
607 Self : Linear<Pair<A, B>>, 964 Self: Linear<Pair<A, B>>,
608 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Linear< 965 for<'a> DiagOp<S::SimpleAdjoint, T::SimpleAdjoint>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
609 Pair<Xʹ,Yʹ>, Codomain=R, 966 {
610 >, 967 type AdjointCodomain = R;
968 type SimpleAdjoint = DiagOp<S::SimpleAdjoint, T::SimpleAdjoint>;
969
970 fn adjoint(&self) -> Self::SimpleAdjoint {
971 DiagOp(self.0.adjoint(), self.1.adjoint())
972 }
973 }
974
975 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A, B>, Pair<Xʹ, Yʹ>> for DiagOp<S, T>
976 where
977 A: Space,
978 B: Space,
979 Xʹ: Space,
980 Yʹ: Space,
981 R: ClosedSpace,
982 S: Preadjointable<A, Xʹ>,
983 T: Preadjointable<B, Yʹ>,
984 Self: Linear<Pair<A, B>>,
985 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>: Linear<Pair<Xʹ, Yʹ>, Codomain = R>,
611 { 986 {
612 type PreadjointCodomain = R; 987 type PreadjointCodomain = R;
613 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 988 type Preadjoint<'a>
989 = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>>
990 where
991 Self: 'a;
614 992
615 fn preadjoint(&self) -> Self::Preadjoint<'_> { 993 fn preadjoint(&self) -> Self::Preadjoint<'_> {
616 DiagOp(self.0.preadjoint(), self.1.preadjoint()) 994 DiagOp(self.0.preadjoint(), self.1.preadjoint())
617 } 995 }
618 } 996 }
619 997
620 /// Block operator 998 /// Block operator
621 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>; 999 pub type BlockOp<S11, S12, S21, S22> = ColOp<RowOp<S11, S12>, RowOp<S21, S22>>;
622
623 1000
624 macro_rules! pairnorm { 1001 macro_rules! pairnorm {
625 ($expj:ty) => { 1002 ($expj:ty) => {
626 impl<F, A, B, S, T, ExpA, ExpB, ExpR> 1003 impl<F, A, B, S, T, ExpA, ExpB, ExpR>
627 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> 1004 BoundedLinear<Pair<A, B>, PairNorm<ExpA, ExpB, $expj>, ExpR, F> for RowOp<S, T>
628 for RowOp<S, T>
629 where 1005 where
630 F : Float, 1006 F: Float,
631 A : Space + Norm<F, ExpA>, 1007 A: Space,
632 B : Space + Norm<F, ExpB>, 1008 B: Space,
633 S : BoundedLinear<A, ExpA, ExpR, F>, 1009 S: BoundedLinear<A, ExpA, ExpR, F>,
634 T : BoundedLinear<B, ExpB, ExpR, F>, 1010 T: BoundedLinear<B, ExpB, ExpR, F>,
635 S::Codomain : Add<T::Codomain>, 1011 S::Codomain: Add<T::Codomain>,
636 <S::Codomain as Add<T::Codomain>>::Output : Space, 1012 <S::Codomain as Add<T::Codomain>>::Output: ClosedSpace,
637 ExpA : NormExponent, 1013 ExpA: NormExponent,
638 ExpB : NormExponent, 1014 ExpB: NormExponent,
639 ExpR : NormExponent, 1015 ExpR: NormExponent,
640 { 1016 {
641 fn opnorm_bound( 1017 fn opnorm_bound(
642 &self, 1018 &self,
643 PairNorm(expa, expb, _) : PairNorm<ExpA, ExpB, $expj>, 1019 PairNorm(expa, expb, _): PairNorm<ExpA, ExpB, $expj>,
644 expr : ExpR 1020 expr: ExpR,
645 ) -> F { 1021 ) -> DynResult<F> {
646 // An application of the triangle inequality bounds the norm by the maximum 1022 // An application of the triangle inequality bounds the norm by the maximum
647 // of the individual norms. A simple observation shows this to be exact. 1023 // of the individual norms. A simple observation shows this to be exact.
648 let na = self.0.opnorm_bound(expa, expr); 1024 let na = self.0.opnorm_bound(expa, expr)?;
649 let nb = self.1.opnorm_bound(expb, expr); 1025 let nb = self.1.opnorm_bound(expb, expr)?;
650 na.max(nb) 1026 Ok(na.max(nb))
651 } 1027 }
652 } 1028 }
653 1029
654 impl<F, A, S, T, ExpA, ExpS, ExpT> 1030 impl<F, A, S, T, ExpA, ExpS, ExpT> BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F>
655 BoundedLinear<A, ExpA, PairNorm<ExpS, ExpT, $expj>, F> 1031 for ColOp<S, T>
656 for ColOp<S, T>
657 where 1032 where
658 F : Float, 1033 F: Float,
659 A : Space + Norm<F, ExpA>, 1034 A: Space,
660 S : BoundedLinear<A, ExpA, ExpS, F>, 1035 S: BoundedLinear<A, ExpA, ExpS, F>,
661 T : BoundedLinear<A, ExpA, ExpT, F>, 1036 T: BoundedLinear<A, ExpA, ExpT, F>,
662 ExpA : NormExponent, 1037 ExpA: NormExponent,
663 ExpS : NormExponent, 1038 ExpS: NormExponent,
664 ExpT : NormExponent, 1039 ExpT: NormExponent,
665 { 1040 {
666 fn opnorm_bound( 1041 fn opnorm_bound(
667 &self, 1042 &self,
668 expa : ExpA, 1043 expa: ExpA,
669 PairNorm(exps, expt, _) : PairNorm<ExpS, ExpT, $expj> 1044 PairNorm(exps, expt, _): PairNorm<ExpS, ExpT, $expj>,
670 ) -> F { 1045 ) -> DynResult<F> {
671 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence, 1046 // This is based on the rule for RowOp and ‖A^*‖ = ‖A‖, hence,
672 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖} 1047 // for A=[S; T], ‖A‖=‖[S^*, T^*]‖ ≤ max{‖S^*‖, ‖T^*‖} = max{‖S‖, ‖T‖}
673 let ns = self.0.opnorm_bound(expa, exps); 1048 let ns = self.0.opnorm_bound(expa, exps)?;
674 let nt = self.1.opnorm_bound(expa, expt); 1049 let nt = self.1.opnorm_bound(expa, expt)?;
675 ns.max(nt) 1050 Ok(ns.max(nt))
676 } 1051 }
677 } 1052 }
678 } 1053 };
679 } 1054 }
680 1055
681 pairnorm!(L1); 1056 pairnorm!(L1);
682 pairnorm!(L2); 1057 pairnorm!(L2);
683 pairnorm!(Linfinity); 1058 pairnorm!(Linfinity);
684 1059
1060 /// The simplest linear mapping, scaling by a scalar.
1061 ///
1062 /// TODO: redefined/replace `Weighted` by composition with [`Scaled`].
1063 pub struct Scaled<F: Float>(pub F);
1064
1065 impl<Domain, F> Mapping<Domain> for Scaled<F>
1066 where
1067 F: Float,
1068 Domain: Space,
1069 Domain::Principal: Mul<F>,
1070 <Domain::Principal as Mul<F>>::Output: ClosedSpace,
1071 {
1072 type Codomain = <Domain::Principal as Mul<F>>::Output;
1073
1074 /// Compute the value of `self` at `x`.
1075 fn apply<I: Instance<Domain>>(&self, x: I) -> Self::Codomain {
1076 x.own() * self.0
1077 }
1078 }
1079
1080 impl<Domain, F> Linear<Domain> for Scaled<F>
1081 where
1082 F: Float,
1083 Domain: Space,
1084 Domain::Principal: Mul<F>,
1085 <Domain::Principal as Mul<F>>::Output: ClosedSpace,
1086 {
1087 }

Mercurial Repository: alg_tools

mercurial