Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 150
c4e394a9c84c
parent 140
26df914dda67
child 151
402d717bb5c0
equal deleted inserted replaced
149:2f1798c65fd6 150:c4e394a9c84c
3 */ 3 */
4 4
5 use crate::direct_product::Pair; 5 use crate::direct_product::Pair;
6 use crate::error::DynResult; 6 use crate::error::DynResult;
7 use crate::instance::Instance; 7 use crate::instance::Instance;
8 pub use crate::mapping::{Composition, DifferentiableImpl, Mapping, Space}; 8 pub use crate::mapping::{ClosedSpace, Composition, DifferentiableImpl, Mapping, Space};
9 use crate::norms::{HasDual, Linfinity, Norm, NormExponent, PairNorm, L1, L2}; 9 use crate::norms::{HasDual, Linfinity, Norm, NormExponent, PairNorm, L1, L2};
10 use crate::types::*; 10 use crate::types::*;
11 use numeric_literals::replace_float_literals; 11 use numeric_literals::replace_float_literals;
12 use serde::Serialize; 12 use serde::Serialize;
13 use std::marker::PhantomData; 13 use std::marker::PhantomData;
23 // fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative { 23 // fn differential_impl<I: Instance<X>>(&self, x: I) -> Self::Derivative {
24 // self.apply(x) 24 // self.apply(x)
25 // } 25 // }
26 // } 26 // }
27 27
28 /// Efficient in-place summation. 28 /// Vector spaces
29 #[replace_float_literals(Self::Field::cast_from(literal))] 29 #[replace_float_literals(Self::Field::cast_from(literal))]
30 pub trait AXPY<X = Self>: 30 pub trait VectorSpace:
31 Space 31 Space<OwnedSpace = Self::Owned>
32 + MulAssign<Self::Field>
33 + DivAssign<Self::Field>
34 + AddAssign<Self>
35 + AddAssign<Self::Owned>
36 + SubAssign<Self>
37 + SubAssign<Self::Owned>
38 + Mul<Self::Field, Output = Self::Owned> 32 + Mul<Self::Field, Output = Self::Owned>
39 + Div<Self::Field, Output = Self::Owned> 33 + Div<Self::Field, Output = Self::Owned>
40 + Add<Self, Output = Self::Owned> 34 + Add<Self, Output = Self::Owned>
41 + Add<Self::Owned, Output = Self::Owned> 35 + Add<Self::Owned, Output = Self::Owned>
42 + Sub<Self, Output = Self::Owned> 36 + Sub<Self, Output = Self::Owned>
43 + Sub<Self::Owned, Output = Self::Owned> 37 + Sub<Self::Owned, Output = Self::Owned>
44 + Neg 38 + Neg
45 where
46 X: Space,
47 { 39 {
48 type Field: Num; 40 type Field: Num;
49 type Owned: AXPY<X, Field = Self::Field>; 41 type Owned: ClosedSpace
50 42 + AXPY<
51 /// Computes `y = βy + αx`, where `y` is `Self`. 43 Self,
52 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field); 44 Field = Self::Field,
53 45 Owned = Self::Owned,
54 /// Copies `x` to `self`. 46 OwnedVariant = Self::Owned,
55 fn copy_from<I: Instance<X>>(&mut self, x: I) { 47 OwnedSpace = Self::Owned,
56 self.axpy(1.0, x, 0.0) 48 >;
57 }
58
59 /// Computes `y = αx`, where `y` is `Self`.
60 fn scale_from<I: Instance<X>>(&mut self, α: Self::Field, x: I) {
61 self.axpy(α, x, 0.0)
62 }
63 49
64 /// Return a similar zero as `self`. 50 /// Return a similar zero as `self`.
65 fn similar_origin(&self) -> Self::Owned; 51 fn similar_origin(&self) -> Self::Owned;
66 // { 52 // {
67 // self.make_origin_generator().make_origin() 53 // self.make_origin_generator().make_origin()
69 55
70 /// Return a similar zero as `x`. 56 /// Return a similar zero as `x`.
71 fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::Owned { 57 fn similar_origin_inst<I: Instance<Self>>(x: I) -> Self::Owned {
72 x.eval(|xr| xr.similar_origin()) 58 x.eval(|xr| xr.similar_origin())
73 } 59 }
60 }
61
62 /// Efficient in-place summation.
63 #[replace_float_literals(Self::Field::cast_from(literal))]
64 pub trait AXPY<X = Self>:
65 VectorSpace
66 + MulAssign<Self::Field>
67 + DivAssign<Self::Field>
68 + AddAssign<Self>
69 + AddAssign<Self::Owned>
70 + SubAssign<Self>
71 + SubAssign<Self::Owned>
72 where
73 X: Space,
74 {
75 /// Computes `y = βy + αx`, where `y` is `Self`.
76 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field);
77
78 /// Copies `x` to `self`.
79 fn copy_from<I: Instance<X>>(&mut self, x: I) {
80 self.axpy(1.0, x, 0.0)
81 }
82
83 /// Computes `y = αx`, where `y` is `Self`.
84 fn scale_from<I: Instance<X>>(&mut self, α: Self::Field, x: I) {
85 self.axpy(α, x, 0.0)
86 }
74 87
75 /// Set self to zero. 88 /// Set self to zero.
76 fn set_zero(&mut self); 89 fn set_zero(&mut self);
77
78 //fn make_origin_generator(&self) -> Self::OriginGen<'_>;
79 } 90 }
80 91
81 /// Efficient in-place application for [`Linear`] operators. 92 /// Efficient in-place application for [`Linear`] operators.
82 #[replace_float_literals(F::cast_from(literal))] 93 #[replace_float_literals(F::cast_from(literal))]
83 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> { 94 pub trait GEMV<F: Num, X: Space, Y = <Self as Mapping<X>>::Codomain>: Linear<X> {
176 pub fn new() -> IdOp<X> { 187 pub fn new() -> IdOp<X> {
177 IdOp(PhantomData) 188 IdOp(PhantomData)
178 } 189 }
179 } 190 }
180 191
181 impl<X: Clone + Space> Mapping<X> for IdOp<X> { 192 impl<X: Space> Mapping<X> for IdOp<X> {
182 type Codomain = X; 193 type Codomain = X::OwnedVariant;
183 194
184 fn apply<I: Instance<X>>(&self, x: I) -> X { 195 fn apply<I: Instance<X>>(&self, x: I) -> X {
185 x.own() 196 x.own()
186 } 197 }
187 } 198 }
188 199
189 impl<X: Clone + Space> Linear<X> for IdOp<X> {} 200 impl<X: Space> Linear<X> for IdOp<X> {}
190 201
191 #[replace_float_literals(F::cast_from(literal))] 202 #[replace_float_literals(F::cast_from(literal))]
192 impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X> 203 impl<F: Num, X, Y> GEMV<F, X, Y> for IdOp<X>
193 where 204 where
194 Y: AXPY<X, Field = F>, 205 Y: AXPY<X, Field = F>,
195 X: Clone + Space, 206 X: Space,
196 { 207 {
197 // Computes `y = αAx + βy`, where `A` is `Self`. 208 // Computes `y = αAx + βy`, where `A` is `Self`.
198 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) { 209 fn gemv<I: Instance<X>>(&self, y: &mut Y, α: F, x: I, β: F) {
199 y.axpy(α, x, β) 210 y.axpy(α, x, β)
200 } 211 }
241 252
242 /// The zero operator from a space to itself 253 /// The zero operator from a space to itself
243 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)] 254 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
244 pub struct SimpleZeroOp; 255 pub struct SimpleZeroOp;
245 256
246 impl<X> Mapping<X> for SimpleZeroOp 257 impl<X: VectorSpace> Mapping<X> for SimpleZeroOp {
247 where
248 X: AXPY + Instance<X>,
249 {
250 type Codomain = X::Owned; 258 type Codomain = X::Owned;
251 259
252 fn apply<I: Instance<X>>(&self, x: I) -> X::Owned { 260 fn apply<I: Instance<X>>(&self, x: I) -> X::Owned {
253 X::similar_origin_inst(x) 261 X::similar_origin_inst(x)
254 } 262 }
255 } 263 }
256 264
257 impl<X> Linear<X> for SimpleZeroOp where X: AXPY + Instance<X> {} 265 impl<X: VectorSpace> Linear<X> for SimpleZeroOp {}
258 266
259 #[replace_float_literals(F::cast_from(literal))] 267 #[replace_float_literals(F::cast_from(literal))]
260 impl<X, Y, F> GEMV<F, X, Y> for SimpleZeroOp 268 impl<X, Y, F> GEMV<F, X, Y> for SimpleZeroOp
261 where 269 where
262 F: Num, 270 F: Num,
263 Y: AXPY<Field = F>, 271 Y: AXPY<Field = F>,
264 X: AXPY<Field = F> + Instance<X>, 272 X: VectorSpace<Field = F> + Instance<X>,
265 { 273 {
266 // Computes `y = αAx + βy`, where `A` is `Self`. 274 // Computes `y = αAx + βy`, where `A` is `Self`.
267 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) { 275 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
268 *y *= β; 276 *y *= β;
269 } 277 }
274 } 282 }
275 283
276 impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp 284 impl<X, F, E1, E2> BoundedLinear<X, E1, E2, F> for SimpleZeroOp
277 where 285 where
278 F: Num, 286 F: Num,
279 X: AXPY<Field = F> + Instance<X> + Norm<E1, F>, 287 X: VectorSpace<Field = F> + Norm<E1, F>,
280 E1: NormExponent, 288 E1: NormExponent,
281 E2: NormExponent, 289 E2: NormExponent,
282 { 290 {
283 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> { 291 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
284 Ok(F::ZERO) 292 Ok(F::ZERO)
286 } 294 }
287 295
288 impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp 296 impl<X, F> Adjointable<X, X::DualSpace> for SimpleZeroOp
289 where 297 where
290 F: Num, 298 F: Num,
291 X: AXPY<Field = F> + Instance<X> + HasDual<F>, 299 X: VectorSpace<Field = F> + HasDual<F>,
292 X::DualSpace: AXPY<Owned = X::DualSpace>, 300 X::DualSpace: VectorSpace<Owned = X::DualSpace>,
293 { 301 {
294 type AdjointCodomain = X::DualSpace; 302 type AdjointCodomain = X::DualSpace;
295 type Adjoint<'b> 303 type Adjoint<'b>
296 = SimpleZeroOp 304 = SimpleZeroOp
297 where 305 where
301 fn adjoint(&self) -> Self::Adjoint<'_> { 309 fn adjoint(&self) -> Self::Adjoint<'_> {
302 SimpleZeroOp 310 SimpleZeroOp
303 } 311 }
304 } 312 }
305 313
306 pub trait OriginGenerator<Y: AXPY> { 314 pub trait OriginGenerator<Y: VectorSpace> {
307 type Ref<'b>: OriginGenerator<Y> 315 type Ref<'b>: OriginGenerator<Y>
308 where 316 where
309 Self: 'b; 317 Self: 'b;
310 318
311 fn origin(&self) -> Y::Owned; 319 fn origin(&self) -> Y::Owned;
312 fn as_ref(&self) -> Self::Ref<'_>; 320 fn as_ref(&self) -> Self::Ref<'_>;
313 } 321 }
314 322
315 impl<Y: AXPY> OriginGenerator<Y> for Y { 323 impl<Y: VectorSpace> OriginGenerator<Y> for Y {
316 type Ref<'b> 324 type Ref<'b>
317 = &'b Y 325 = &'b Y
318 where 326 where
319 Self: 'b; 327 Self: 'b;
320 328
327 fn as_ref(&self) -> Self::Ref<'_> { 335 fn as_ref(&self) -> Self::Ref<'_> {
328 self 336 self
329 } 337 }
330 } 338 }
331 339
332 impl<'b, Y: AXPY> OriginGenerator<Y> for &'b Y { 340 impl<'b, Y: VectorSpace> OriginGenerator<Y> for &'b Y {
333 type Ref<'c> 341 type Ref<'c>
334 = Self 342 = Self
335 where 343 where
336 Self: 'c; 344 Self: 'c;
337 345
346 } 354 }
347 } 355 }
348 356
349 /// A zero operator that can be eitherh dualised or predualised (once). 357 /// A zero operator that can be eitherh dualised or predualised (once).
350 /// This is achieved by storing an oppropriate zero. 358 /// This is achieved by storing an oppropriate zero.
351 pub struct ZeroOp<X, Y: AXPY<Field = F>, OY: OriginGenerator<Y>, O, F: Float = f64> { 359 pub struct ZeroOp<X, Y: VectorSpace<Field = F>, OY: OriginGenerator<Y>, O, F: Float = f64> {
352 codomain_origin_generator: OY, 360 codomain_origin_generator: OY,
353 other_origin_generator: O, 361 other_origin_generator: O,
354 _phantoms: PhantomData<(X, Y, F)>, 362 _phantoms: PhantomData<(X, Y, F)>,
355 } 363 }
356 364
357 impl<X, Y, OY, F> ZeroOp<X, Y, OY, (), F> 365 impl<X, Y, OY, F> ZeroOp<X, Y, OY, (), F>
358 where 366 where
359 OY: OriginGenerator<Y>, 367 OY: OriginGenerator<Y>,
360 X: AXPY<Field = F>, 368 X: VectorSpace<Field = F>,
361 Y: AXPY<Field = F>, 369 Y: VectorSpace<Field = F>,
362 F: Float, 370 F: Float,
363 { 371 {
364 pub fn new(y_og: OY) -> Self { 372 pub fn new(y_og: OY) -> Self {
365 ZeroOp { 373 ZeroOp {
366 codomain_origin_generator: y_og, 374 codomain_origin_generator: y_og,
375 OY: OriginGenerator<Y>, 383 OY: OriginGenerator<Y>,
376 OXprime: OriginGenerator<Xprime>, 384 OXprime: OriginGenerator<Xprime>,
377 X: HasDual<F, DualSpace = Xprime>, 385 X: HasDual<F, DualSpace = Xprime>,
378 Y: HasDual<F, DualSpace = Yprime>, 386 Y: HasDual<F, DualSpace = Yprime>,
379 F: Float, 387 F: Float,
380 Xprime: AXPY<Field = F, Owned = Xprime>, 388 Xprime: VectorSpace<Field = F, Owned = Xprime>,
381 Xprime::Owned: AXPY<Field = F>, 389 Xprime::Owned: AXPY<Field = F>,
382 Yprime: Space + Instance<Yprime>, 390 Yprime: Space + Instance<Yprime>,
383 { 391 {
384 pub fn new_dualisable(y_og: OY, xprime_og: OXprime) -> Self { 392 pub fn new_dualisable(y_og: OY, xprime_og: OXprime) -> Self {
385 ZeroOp { 393 ZeroOp {
390 } 398 }
391 } 399 }
392 400
393 impl<X, Y, O, OY, F> Mapping<X> for ZeroOp<X, Y, OY, O, F> 401 impl<X, Y, O, OY, F> Mapping<X> for ZeroOp<X, Y, OY, O, F>
394 where 402 where
395 X: Space + Instance<X>, 403 X: Space,
396 Y: AXPY<Field = F>, 404 Y: VectorSpace<Field = F>,
397 F: Float, 405 F: Float,
398 OY: OriginGenerator<Y>, 406 OY: OriginGenerator<Y>,
399 { 407 {
400 type Codomain = Y::Owned; 408 type Codomain = Y::Owned;
401 409
404 } 412 }
405 } 413 }
406 414
407 impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F> 415 impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F>
408 where 416 where
409 X: Space + Instance<X>, 417 X: Space,
410 Y: AXPY<Field = F>, 418 Y: VectorSpace<Field = F>,
411 F: Float, 419 F: Float,
412 OY: OriginGenerator<Y>, 420 OY: OriginGenerator<Y>,
413 { 421 {
414 } 422 }
415 423
416 #[replace_float_literals(F::cast_from(literal))] 424 #[replace_float_literals(F::cast_from(literal))]
417 impl<X, Y, OY, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, OY, O, F> 425 impl<X, Y, OY, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, OY, O, F>
418 where 426 where
419 X: Space + Instance<X>, 427 X: Space,
420 Y: AXPY<Field = F, Owned = Y>, 428 Y: AXPY<Field = F, Owned = Y>,
421 F: Float, 429 F: Float,
422 OY: OriginGenerator<Y>, 430 OY: OriginGenerator<Y>,
423 { 431 {
424 // Computes `y = αAx + βy`, where `A` is `Self`. 432 // Computes `y = αAx + βy`, where `A` is `Self`.
432 } 440 }
433 441
434 impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F> 442 impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F>
435 where 443 where
436 X: Space + Instance<X> + Norm<E1, F>, 444 X: Space + Instance<X> + Norm<E1, F>,
437 Y: AXPY<Field = F>, 445 Y: VectorSpace<Field = F>,
438 Y::Owned: Clone, 446 Y::Owned: Clone,
439 F: Float, 447 F: Float,
440 E1: NormExponent, 448 E1: NormExponent,
441 E2: NormExponent, 449 E2: NormExponent,
442 OY: OriginGenerator<Y>, 450 OY: OriginGenerator<Y>,
450 for ZeroOp<X, Y, OY, OXprime, F> 458 for ZeroOp<X, Y, OY, OXprime, F>
451 where 459 where
452 X: HasDual<F, DualSpace = Xprime>, 460 X: HasDual<F, DualSpace = Xprime>,
453 Y: HasDual<F, DualSpace = Yprime>, 461 Y: HasDual<F, DualSpace = Yprime>,
454 F: Float, 462 F: Float,
455 Xprime: AXPY<Field = F, Owned = Xprime>, 463 Xprime: VectorSpace<Field = F, Owned = Xprime>,
456 Xprime::Owned: AXPY<Field = F>, 464 Xprime::Owned: AXPY<Field = F>,
457 Yprime: Space + Instance<Yprime>, 465 Yprime: Space + Instance<Yprime>,
458 OY: OriginGenerator<Y>, 466 OY: OriginGenerator<Y>,
459 OXprime: OriginGenerator<X::DualSpace>, 467 OXprime: OriginGenerator<X::DualSpace>,
460 { 468 {

Mercurial Repository: alg_tools

mercurial