src/linops.rs

branch
dev
changeset 78
cebedc4a8331
parent 61
05089fbc0310
child 79
d63e40672dd6
equal deleted inserted replaced
61:05089fbc0310 78:cebedc4a8331
7 use serde::Serialize; 7 use serde::Serialize;
8 use crate::types::*; 8 use crate::types::*;
9 pub use crate::mapping::{Mapping, Space, Composition}; 9 pub use crate::mapping::{Mapping, Space, Composition};
10 use crate::direct_product::Pair; 10 use crate::direct_product::Pair;
11 use crate::instance::Instance; 11 use crate::instance::Instance;
12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm}; 12 use crate::norms::{NormExponent, PairNorm, L1, L2, Linfinity, Norm, HasDual};
13 13
14 /// Trait for linear operators on `X`. 14 /// Trait for linear operators on `X`.
15 pub trait Linear<X : Space> : Mapping<X> 15 pub trait Linear<X : Space> : Mapping<X>
16 { } 16 { }
17 17
78 self.apply(x) 78 self.apply(x)
79 } 79 }
80 }*/ 80 }*/
81 81
82 /// Trait for forming the adjoint operator of `Self`. 82 /// Trait for forming the adjoint operator of `Self`.
83 pub trait Adjointable<X, Yʹ> : Linear<X> 83 pub trait Adjointable<X, F : Float = f64> : Linear<X>
84 where 84 where
85 X : Space, 85 X : HasDual<F>,
86 Yʹ : Space, 86 Self::Codomain : HasDual<F>,
87 { 87 {
88 type AdjointCodomain : Space; 88 type AdjointCodomain : Instance<X::DualSpace>;
89 type Adjoint<'a> : Linear<Yʹ, Codomain=Self::AdjointCodomain> where Self : 'a; 89 type Adjoint<'a> : Linear<
90 <Self::Codomain as HasDual<F>>::DualSpace,
91 Codomain=Self::AdjointCodomain,
92 > where Self : 'a;
90 93
91 /// Form the adjoint operator of `self`. 94 /// Form the adjoint operator of `self`.
92 fn adjoint(&self) -> Self::Adjoint<'_>; 95 fn adjoint(&self) -> Self::Adjoint<'_>;
93 } 96 }
94 97
97 /// For an operator $A$ this is an operator $A\_\*$ 100 /// For an operator $A$ this is an operator $A\_\*$
98 /// such that its adjoint $(A\_\*)^\*=A$. The space `X` is the domain of the `Self` 101 /// such that its adjoint $(A\_\*)^\*=A$. The space `X` is the domain of the `Self`
99 /// operator. The space `Ypre` is the predual of its codomain, and should be the 102 /// operator. The space `Ypre` is the predual of its codomain, and should be the
100 /// domain of the adjointed operator. `Self::Preadjoint` should be 103 /// domain of the adjointed operator. `Self::Preadjoint` should be
101 /// [`Adjointable`]`<'a,Ypre,X>`. 104 /// [`Adjointable`]`<'a,Ypre,X>`.
102 105 ///
103 pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X> { 106 /// We do not set further restrictions on the spacds, to allow preadjointing when `X`
107 /// is on the dual space of `Ypre`, but a subset of it.
108 pub trait Preadjointable<X : Space, Ypre : Space> : Linear<X>
109 //where
110 // Ypre : HasDual<F, DualSpace=Self::Codomain>,
111 {
104 type PreadjointCodomain : Space; 112 type PreadjointCodomain : Space;
105 type Preadjoint<'a> : Adjointable< 113 type Preadjoint<'a> : Linear<Ypre, Codomain=Self::PreadjointCodomain> where Self : 'a;
106 Ypre, X, AdjointCodomain=Self::Codomain, Codomain=Self::PreadjointCodomain
107 > where Self : 'a;
108 114
109 /// Form the adjoint operator of `self`. 115 /// Form the adjoint operator of `self`.
110 fn preadjoint(&self) -> Self::Preadjoint<'_>; 116 fn preadjoint(&self) -> Self::Preadjoint<'_>;
111 } 117 }
112 118
113 /// Adjointable operators $A: X → Y$ between reflexive spaces $X$ and $Y$.
114 pub trait SimplyAdjointable<X : Space> : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
115 impl<'a,X : Space, T> SimplyAdjointable<X> for T
116 where T : Adjointable<X,<Self as Mapping<X>>::Codomain> {}
117 119
118 /// The identity operator 120 /// The identity operator
119 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)] 121 #[derive(Clone,Copy,Debug,Serialize,Eq,PartialEq)]
120 pub struct IdOp<X> (PhantomData<X>); 122 pub struct IdOp<X> (PhantomData<X>);
121 123
157 E : NormExponent 159 E : NormExponent
158 { 160 {
159 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE } 161 fn opnorm_bound(&self, _xexp : E, _codexp : E) -> F { F::ONE }
160 } 162 }
161 163
162 impl<X : Clone + Space> Adjointable<X,X> for IdOp<X> { 164 impl<X : Clone + HasDual<F, DualSpace=X>, F : Float> Adjointable<X, F> for IdOp<X> {
163 type AdjointCodomain=X; 165 type AdjointCodomain = X;
164 type Adjoint<'a> = IdOp<X> where X : 'a; 166 type Adjoint<'a> = IdOp<X::DualSpace> where X : 'a;
165 167
166 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() } 168 fn adjoint(&self) -> Self::Adjoint<'_> { IdOp::new() }
167 } 169 }
168 170
169 impl<X : Clone + Space> Preadjointable<X,X> for IdOp<X> { 171 impl<X : Clone + Space> Preadjointable<X, X> for IdOp<X> {
170 type PreadjointCodomain=X; 172 type PreadjointCodomain = X;
171 type Preadjoint<'a> = IdOp<X> where X : 'a; 173 type Preadjoint<'a> = IdOp<X> where X : 'a;
172 174
173 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() } 175 fn preadjoint(&self) -> Self::Preadjoint<'_> { IdOp::new() }
174 } 176 }
175 177
330 self.1.apply_add(&mut y.1, x); 332 self.1.apply_add(&mut y.1, x);
331 } 333 }
332 } 334 }
333 335
334 336
335 impl<A, B, Yʹ, S, T> Adjointable<Pair<A,B>, Yʹ> for RowOp<S, T> 337 impl<F, A, B, Yʹ, S, T> Adjointable<Pair<A, B>, F> for RowOp<S, T>
336 where 338 where
337 A : Space, 339 F : Float,
338 B : Space, 340 A : HasDual<F>,
341 B : HasDual<F>,
342 S : Adjointable<A, F>,
343 T : Adjointable<B, F>,
339 Yʹ : Space, 344 Yʹ : Space,
340 S : Adjointable<A, Yʹ>, 345 S :: Codomain : HasDual<F, DualSpace=Yʹ>,
341 T : Adjointable<B, Yʹ>, 346 T :: Codomain : HasDual<F, DualSpace=Yʹ>,
342 Self : Linear<Pair<A, B>>, 347 S::Codomain : Add<T::Codomain>,
348 <S::Codomain as Add<T::Codomain>>::Output : HasDual<F, DualSpace=Yʹ>,
349 Self::Codomain : HasDual<F, DualSpace=Yʹ>,
350 //Self : Linear<Pair<A, B>>,
343 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 351 // for<'a> ColOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
344 // Yʹ, 352 // Yʹ,
345 // Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain> 353 // //Codomain=Pair<S::AdjointCodomain, T::AdjointCodomain>
346 // >, 354 // >,
347 { 355 {
348 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>; 356 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
349 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 357 type Adjoint<'a> = ColOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
350 358
358 A : Space, 366 A : Space,
359 B : Space, 367 B : Space,
360 Yʹ : Space, 368 Yʹ : Space,
361 S : Preadjointable<A, Yʹ>, 369 S : Preadjointable<A, Yʹ>,
362 T : Preadjointable<B, Yʹ>, 370 T : Preadjointable<B, Yʹ>,
363 Self : Linear<Pair<A, B>>, 371 S::Codomain : Add<T::Codomain>,
364 for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable< 372 <S::Codomain as Add<T::Codomain>>::Output : Space,
365 Yʹ, Pair<A,B>, 373 //Self : Linear<Pair<A, B>, Codomain=Y>,
366 Codomain=Pair<S::PreadjointCodomain, T::PreadjointCodomain>, 374 //for<'a> ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Yʹ, F>,
367 AdjointCodomain = Self::Codomain,
368 >,
369 { 375 {
370 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>; 376 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
371 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 377 type Preadjoint<'a> = ColOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
372 378
373 fn preadjoint(&self) -> Self::Preadjoint<'_> { 379 fn preadjoint(&self) -> Self::Preadjoint<'_> {
374 ColOp(self.0.preadjoint(), self.1.preadjoint()) 380 ColOp(self.0.preadjoint(), self.1.preadjoint())
375 } 381 }
376 } 382 }
377 383
378 384
379 impl<A, Xʹ, Yʹ, R, S, T> Adjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T> 385 impl<F, A, Aʹ, S, T> Adjointable<A, F> for ColOp<S, T>
380 where 386 where
381 A : Space, 387 F : Float,
388 A : HasDual<F>,
389 S : Adjointable<A, F>,
390 T : Adjointable<A, F>,
391 T::Codomain : HasDual<F>,
392 S::Codomain : HasDual<F>,
393 Aʹ : Space + Instance<A::DualSpace>,
394 <S as Adjointable<A, F>>::AdjointCodomain : Add<<T as Adjointable<A, F>>::AdjointCodomain, Output=Aʹ>,
395 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear<
396 // Pair<<T::Codomain as HasDual<F>>::DualSpace, <S::Codomain as HasDual<F>>::DualSpace>,
397 // Codomain=Aʹ
398 // >,
399 {
400 type AdjointCodomain = Aʹ;
401 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
402
403 fn adjoint(&self) -> Self::Adjoint<'_> {
404 RowOp(self.0.adjoint(), self.1.adjoint())
405 }
406 }
407
408 impl<A, Aʹ, Xʹ, Yʹ, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
409 where
410 A : Space,
411 Aʹ : Space,
382 Xʹ : Space, 412 Xʹ : Space,
383 Yʹ : Space, 413 Yʹ : Space,
384 R : Space + ClosedAdd, 414 S : Preadjointable<A, Xʹ, PreadjointCodomain=Aʹ>,
385 S : Adjointable<A, Xʹ, AdjointCodomain = R>, 415 T : Preadjointable<A, Yʹ, PreadjointCodomain=Aʹ>,
386 T : Adjointable<A, Yʹ, AdjointCodomain = R>, 416 Aʹ : ClosedAdd,
387 Self : Linear<A>, 417 //for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<Pair<Xʹ,Yʹ>, F>,
388 // for<'a> RowOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 418 {
389 // Pair<Xʹ,Yʹ>, 419 type PreadjointCodomain = Aʹ;
390 // Codomain=R,
391 // >,
392 {
393 type AdjointCodomain = R;
394 type Adjoint<'a> = RowOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
395
396 fn adjoint(&self) -> Self::Adjoint<'_> {
397 RowOp(self.0.adjoint(), self.1.adjoint())
398 }
399 }
400
401 impl<A, Xʹ, Yʹ, R, S, T> Preadjointable<A,Pair<Xʹ,Yʹ>> for ColOp<S, T>
402 where
403 A : Space,
404 Xʹ : Space,
405 Yʹ : Space,
406 R : Space + ClosedAdd,
407 S : Preadjointable<A, Xʹ, PreadjointCodomain = R>,
408 T : Preadjointable<A, Yʹ, PreadjointCodomain = R>,
409 Self : Linear<A>,
410 for<'a> RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable<
411 Pair<Xʹ,Yʹ>, A,
412 Codomain = R,
413 AdjointCodomain = Self::Codomain,
414 >,
415 {
416 type PreadjointCodomain = R;
417 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 420 type Preadjoint<'a> = RowOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
418 421
419 fn preadjoint(&self) -> Self::Preadjoint<'_> { 422 fn preadjoint(&self) -> Self::Preadjoint<'_> {
420 RowOp(self.0.preadjoint(), self.1.preadjoint()) 423 RowOp(self.0.preadjoint(), self.1.preadjoint())
421 } 424 }
476 self.0.apply_add(&mut y.0, u); 479 self.0.apply_add(&mut y.0, u);
477 self.1.apply_add(&mut y.1, v); 480 self.1.apply_add(&mut y.1, v);
478 } 481 }
479 } 482 }
480 483
481 impl<A, B, Xʹ, Yʹ, R, S, T> Adjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 484 impl<F, A, B, S, T> Adjointable<Pair<A,B>, F> for DiagOp<S, T>
482 where 485 where
483 A : Space, 486 F: Float,
484 B : Space, 487 A : HasDual<F>,
485 Xʹ: Space, 488 B : HasDual<F>,
486 Yʹ : Space, 489 S : Adjointable<A, F>,
487 R : Space, 490 T : Adjointable<B, F>,
488 S : Adjointable<A, Xʹ>, 491 T::Codomain : HasDual<F>,
489 T : Adjointable<B, Yʹ>, 492 S::Codomain : HasDual<F>,
490 Self : Linear<Pair<A, B>>, 493 {
491 for<'a> DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> : Linear< 494 type AdjointCodomain = Pair<S::AdjointCodomain, T::AdjointCodomain>;
492 Pair<Xʹ,Yʹ>, Codomain=R,
493 >,
494 {
495 type AdjointCodomain = R;
496 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a; 495 type Adjoint<'a> = DiagOp<S::Adjoint<'a>, T::Adjoint<'a>> where Self : 'a;
497 496
498 fn adjoint(&self) -> Self::Adjoint<'_> { 497 fn adjoint(&self) -> Self::Adjoint<'_> {
499 DiagOp(self.0.adjoint(), self.1.adjoint()) 498 DiagOp(self.0.adjoint(), self.1.adjoint())
500 } 499 }
501 } 500 }
502 501
503 impl<A, B, Xʹ, Yʹ, R, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T> 502 impl<A, B, Xʹ, Yʹ, S, T> Preadjointable<Pair<A,B>, Pair<Xʹ,Yʹ>> for DiagOp<S, T>
504 where 503 where
505 A : Space, 504 A : Space,
506 B : Space, 505 B : Space,
507 Xʹ: Space, 506 Xʹ : Space,
508 Yʹ : Space, 507 Yʹ : Space,
509 R : Space,
510 S : Preadjointable<A, Xʹ>, 508 S : Preadjointable<A, Xʹ>,
511 T : Preadjointable<B, Yʹ>, 509 T : Preadjointable<B, Yʹ>,
512 Self : Linear<Pair<A, B>>, 510 {
513 for<'a> DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> : Adjointable< 511 type PreadjointCodomain = Pair<S::PreadjointCodomain, T::PreadjointCodomain>;
514 Pair<Xʹ,Yʹ>, Pair<A, B>,
515 Codomain=R,
516 AdjointCodomain = Self::Codomain,
517 >,
518 {
519 type PreadjointCodomain = R;
520 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a; 512 type Preadjoint<'a> = DiagOp<S::Preadjoint<'a>, T::Preadjoint<'a>> where Self : 'a;
521 513
522 fn preadjoint(&self) -> Self::Preadjoint<'_> { 514 fn preadjoint(&self) -> Self::Preadjoint<'_> {
523 DiagOp(self.0.preadjoint(), self.1.preadjoint()) 515 DiagOp(self.0.preadjoint(), self.1.preadjoint())
524 } 516 }

mercurial