Mercurial > repos > alg_tools / file comparison

comparison: src/linops.rs

src/linops.rs

branch
dev
changeset 140
26df914dda67
parent 139
f78885441218
child 150
c4e394a9c84c
equal deleted inserted replaced
139:f78885441218 140:26df914dda67
45 where 45 where
46 X: Space, 46 X: Space,
47 { 47 {
48 type Field: Num; 48 type Field: Num;
49 type Owned: AXPY<X, Field = Self::Field>; 49 type Owned: AXPY<X, Field = Self::Field>;
50 // type OriginGen<'a>: OriginGenerator<Self::Owned, F>
51 // where
52 // Self: 'a;
53 50
54 /// Computes `y = βy + αx`, where `y` is `Self`. 51 /// Computes `y = βy + αx`, where `y` is `Self`.
55 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field); 52 fn axpy<I: Instance<X>>(&mut self, α: Self::Field, x: I, β: Self::Field);
56 53
57 /// Copies `x` to `self`. 54 /// Copies `x` to `self`.
304 fn adjoint(&self) -> Self::Adjoint<'_> { 301 fn adjoint(&self) -> Self::Adjoint<'_> {
305 SimpleZeroOp 302 SimpleZeroOp
306 } 303 }
307 } 304 }
308 305
306 pub trait OriginGenerator<Y: AXPY> {
307 type Ref<'b>: OriginGenerator<Y>
308 where
309 Self: 'b;
310
311 fn origin(&self) -> Y::Owned;
312 fn as_ref(&self) -> Self::Ref<'_>;
313 }
314
315 impl<Y: AXPY> OriginGenerator<Y> for Y {
316 type Ref<'b>
317 = &'b Y
318 where
319 Self: 'b;
320
321 #[inline]
322 fn origin(&self) -> Y::Owned {
323 return self.similar_origin();
324 }
325
326 #[inline]
327 fn as_ref(&self) -> Self::Ref<'_> {
328 self
329 }
330 }
331
332 impl<'b, Y: AXPY> OriginGenerator<Y> for &'b Y {
333 type Ref<'c>
334 = Self
335 where
336 Self: 'c;
337
338 #[inline]
339 fn origin(&self) -> Y::Owned {
340 return self.similar_origin();
341 }
342
343 #[inline]
344 fn as_ref(&self) -> Self::Ref<'_> {
345 self
346 }
347 }
348
309 /// A zero operator that can be eitherh dualised or predualised (once). 349 /// A zero operator that can be eitherh dualised or predualised (once).
310 /// This is achieved by storing an oppropriate zero. 350 /// This is achieved by storing an oppropriate zero.
311 pub struct ZeroOp<X, Y: AXPY<Field = F>, C, O, F: Float = f64> { 351 pub struct ZeroOp<X, Y: AXPY<Field = F>, OY: OriginGenerator<Y>, O, F: Float = f64> {
312 codomain_sample: C, 352 codomain_origin_generator: OY,
313 other_sample: O, 353 other_origin_generator: O,
314 _phantoms: PhantomData<(X, Y, F)>, 354 _phantoms: PhantomData<(X, Y, F)>,
315 } 355 }
316 356
317 impl<'b, X, Y, F> ZeroOp<X, Y, &'b Y, &'b X::DualSpace, F> 357 impl<X, Y, OY, F> ZeroOp<X, Y, OY, (), F>
318 where 358 where
319 X: HasDual<F>, 359 OY: OriginGenerator<Y>,
320 Y: HasDual<F>, 360 X: AXPY<Field = F>,
321 Y::Owned: Clone, 361 Y: AXPY<Field = F>,
322 F: Float, 362 F: Float,
323 { 363 {
324 pub fn new_dualisable(x_dual_sample: &'b X::DualSpace, y_sample: &'b Y) -> Self { 364 pub fn new(y_og: OY) -> Self {
325 ZeroOp { 365 ZeroOp {
326 codomain_sample: y_sample, 366 codomain_origin_generator: y_og,
327 other_sample: x_dual_sample, 367 other_origin_generator: (),
328 _phantoms: PhantomData, 368 _phantoms: PhantomData,
329 } 369 }
330 } 370 }
331 } 371 }
332 372
333 impl<'b, X, Y, O, F> Mapping<X> for ZeroOp<X, Y, &'b Y, O, F> 373 impl<X, Y, OY, OXprime, Xprime, Yprime, F> ZeroOp<X, Y, OY, OXprime, F>
374 where
375 OY: OriginGenerator<Y>,
376 OXprime: OriginGenerator<Xprime>,
377 X: HasDual<F, DualSpace = Xprime>,
378 Y: HasDual<F, DualSpace = Yprime>,
379 F: Float,
380 Xprime: AXPY<Field = F, Owned = Xprime>,
381 Xprime::Owned: AXPY<Field = F>,
382 Yprime: Space + Instance<Yprime>,
383 {
384 pub fn new_dualisable(y_og: OY, xprime_og: OXprime) -> Self {
385 ZeroOp {
386 codomain_origin_generator: y_og,
387 other_origin_generator: xprime_og,
388 _phantoms: PhantomData,
389 }
390 }
391 }
392
393 impl<X, Y, O, OY, F> Mapping<X> for ZeroOp<X, Y, OY, O, F>
334 where 394 where
335 X: Space + Instance<X>, 395 X: Space + Instance<X>,
336 Y: AXPY<Field = F>, 396 Y: AXPY<Field = F>,
337 F: Float, 397 F: Float,
398 OY: OriginGenerator<Y>,
338 { 399 {
339 type Codomain = Y::Owned; 400 type Codomain = Y::Owned;
340 401
341 fn apply<I: Instance<X>>(&self, _x: I) -> Y::Owned { 402 fn apply<I: Instance<X>>(&self, _x: I) -> Y::Owned {
342 self.codomain_sample.similar_origin() 403 self.codomain_origin_generator.origin()
343 } 404 }
344 } 405 }
345 406
346 impl<'b, X, Y, O, F> Linear<X> for ZeroOp<X, Y, &'b Y, O, F> 407 impl<X, Y, OY, O, F> Linear<X> for ZeroOp<X, Y, OY, O, F>
347 where 408 where
348 X: Space + Instance<X>, 409 X: Space + Instance<X>,
349 Y: AXPY<Field = F>, 410 Y: AXPY<Field = F>,
350 F: Float, 411 F: Float,
412 OY: OriginGenerator<Y>,
351 { 413 {
352 } 414 }
353 415
354 #[replace_float_literals(F::cast_from(literal))] 416 #[replace_float_literals(F::cast_from(literal))]
355 impl<'b, X, Y, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, &'b Y, O, F> 417 impl<X, Y, OY, O, F> GEMV<F, X, Y> for ZeroOp<X, Y, OY, O, F>
356 where 418 where
357 X: Space + Instance<X>, 419 X: Space + Instance<X>,
358 Y: AXPY<Field = F>, 420 Y: AXPY<Field = F, Owned = Y>,
359 F: Float, 421 F: Float,
422 OY: OriginGenerator<Y>,
360 { 423 {
361 // Computes `y = αAx + βy`, where `A` is `Self`. 424 // Computes `y = αAx + βy`, where `A` is `Self`.
362 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) { 425 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
363 *y *= β; 426 *y *= β;
364 } 427 }
366 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) { 429 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
367 y.set_zero(); 430 y.set_zero();
368 } 431 }
369 } 432 }
370 433
371 impl<'b, X, Y, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, &'b Y, O, F> 434 impl<X, Y, OY, O, F, E1, E2> BoundedLinear<X, E1, E2, F> for ZeroOp<X, Y, OY, O, F>
372 where 435 where
373 X: Space + Instance<X> + Norm<E1, F>, 436 X: Space + Instance<X> + Norm<E1, F>,
374 Y: AXPY<Field = F>, 437 Y: AXPY<Field = F>,
375 Y::Owned: Clone, 438 Y::Owned: Clone,
376 F: Float, 439 F: Float,
377 E1: NormExponent, 440 E1: NormExponent,
378 E2: NormExponent, 441 E2: NormExponent,
442 OY: OriginGenerator<Y>,
379 { 443 {
380 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> { 444 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
381 Ok(F::ZERO) 445 Ok(F::ZERO)
382 } 446 }
383 } 447 }
384 448
385 impl<'b, X, Y, Xprime, Yprime, F> Adjointable<X, Yprime> for ZeroOp<X, Y, &'b Y, &'b Xprime, F> 449 impl<'b, X, Y, OY, OXprime, Xprime, Yprime, F> Adjointable<X, Yprime>
450 for ZeroOp<X, Y, OY, OXprime, F>
386 where 451 where
387 X: HasDual<F, DualSpace = Xprime>, 452 X: HasDual<F, DualSpace = Xprime>,
388 Y: HasDual<F, DualSpace = Yprime>, 453 Y: HasDual<F, DualSpace = Yprime>,
389 F: Float, 454 F: Float,
390 Xprime: AXPY<Field = F, Owned = Xprime>, 455 Xprime: AXPY<Field = F, Owned = Xprime>,
391 Xprime::Owned: AXPY<Field = F>, 456 Xprime::Owned: AXPY<Field = F>,
392 Yprime: Space + Instance<Yprime>, 457 Yprime: Space + Instance<Yprime>,
458 OY: OriginGenerator<Y>,
459 OXprime: OriginGenerator<X::DualSpace>,
393 { 460 {
394 type AdjointCodomain = Xprime; 461 type AdjointCodomain = Xprime;
395 type Adjoint<'c> 462 type Adjoint<'c>
396 = ZeroOp<Yprime, Xprime, &'b Xprime, (), F> 463 = ZeroOp<Yprime, Xprime, OXprime::Ref<'c>, (), F>
397 where 464 where
398 Self: 'c; 465 Self: 'c;
399 // () means not (pre)adjointable. 466 // () means not (pre)adjointable.
400 467
401 fn adjoint(&self) -> Self::Adjoint<'_> { 468 fn adjoint(&self) -> Self::Adjoint<'_> {
402 ZeroOp { 469 ZeroOp {
403 codomain_sample: self.other_sample, 470 codomain_origin_generator: self.other_origin_generator.as_ref(),
404 other_sample: (), 471 other_origin_generator: (),
405 _phantoms: PhantomData, 472 _phantoms: PhantomData,
406 } 473 }
407 } 474 }
408 } 475 }
409
410 /*
411 /// Trait for forming origins (zeroes) in vector spaces
412 pub trait OriginGenerator<X, F: Num = f64> {
413 fn make_origin(&self) -> X;
414 }
415
416 /// Origin generator for statically sized Euclidean spaces.
417 pub struct StaticEuclideanOriginGenerator;
418
419 impl<X, F: Float> OriginGenerator<X, F> for StaticEuclideanOriginGenerator
420 where
421 X: StaticEuclidean<F> + Euclidean<F, Output = X>,
422 {
423 #[inline]
424 fn make_origin(&self) -> X {
425 X::origin()
426 }
427 }
428
429 /// Origin generator arbitrary spaces that implement [`AXPY`], based on a sample vector.
430 pub struct SimilarOriginGenerator<'a, X>(&'a X);
431
432 impl<'a, X, F: Float> OriginGenerator<X, F> for SimilarOriginGenerator<'a, X>
433 where
434 X: AXPY<F, Owned = X>,
435 {
436 #[inline]
437 fn make_origin(&self) -> X {
438 self.0.similar_origin()
439 }
440 }
441
442 pub struct NotAnOriginGenerator;
443
444 /// The zero operator
445 #[derive(Clone, Copy, Debug, Serialize, Eq, PartialEq)]
446 pub struct ZeroOp<X, Y, YOrigin, XDOrigin = NotAnOriginGenerator, F: Num = f64> {
447 y_origin: YOrigin,
448 xd_origin: XDOrigin,
449 _phantoms: PhantomData<(X, Y, F)>,
450 }
451
452 // TODO: Need to make Zero in Instance.
453
454 impl<X, Y, F, YOrigin, XDOrigin> ZeroOp<X, Y, YOrigin, XDOrigin, F>
455 where
456 F: Num,
457 YOrigin: OriginGenerator<Y, F>,
458 {
459 pub fn new(y_origin: YOrigin, xd_origin: XDOrigin) -> Self {
460 ZeroOp {
461 y_origin,
462 xd_origin,
463 _phantoms: PhantomData,
464 }
465 }
466 }
467
468 impl<X, Y, F, YOrigin, XDOrigin> Mapping<X> for ZeroOp<X, Y, YOrigin, XDOrigin, F>
469 where
470 F: Num,
471 YOrigin: OriginGenerator<Y, F>,
472 Y: Space,
473 X: Space + Instance<X>,
474 {
475 type Codomain = Y;
476
477 fn apply<I: Instance<X>>(&self, _: I) -> Y {
478 self.y_origin.make_origin()
479 }
480 }
481
482 impl<X, Y, F, YOrigin, XDOrigin> Linear<X> for ZeroOp<X, Y, YOrigin, XDOrigin, F>
483 where
484 F: Num,
485 YOrigin: OriginGenerator<Y, F>,
486 Y: Space,
487 X: Space + Instance<X>,
488 {
489 }
490
491 #[replace_float_literals(F::cast_from(literal))]
492 impl<X, Y, F, YOrigin, XDOrigin> GEMV<F, X, Y> for ZeroOp<X, Y, YOrigin, XDOrigin, F>
493 where
494 F: Num,
495 YOrigin: OriginGenerator<Y, F>,
496 Y: Space + AXPY<F>,
497 X: Space + Instance<X>,
498 {
499 // Computes `y = αAx + βy`, where `A` is `Self`.
500 fn gemv<I: Instance<X>>(&self, y: &mut Y, _α: F, _x: I, β: F) {
501 *y *= β;
502 }
503
504 fn apply_mut<I: Instance<X>>(&self, y: &mut Y, _x: I) {
505 y.set_zero();
506 }
507 }
508
509 impl<X, Y, F, YOrigin, XDOrigin, E1, E2> BoundedLinear<X, E1, E2, F>
510 for ZeroOp<X, Y, YOrigin, XDOrigin, F>
511 where
512 F: Num,
513 YOrigin: OriginGenerator<Y, F>,
514 Y: Space + AXPY<F> + Norm<E2, F>,
515 X: Space + Instance<X> + Norm<E1, F>,
516 E1: NormExponent,
517 E2: NormExponent,
518 {
519 fn opnorm_bound(&self, _xexp: E1, _codexp: E2) -> DynResult<F> {
520 Ok(F::ZERO)
521 }
522 }
523
524 impl<X, Y, Yprime, F, YOrigin, XDOrigin> Adjointable<X, Yprime>
525 for ZeroOp<X, Y, YOrigin, XDOrigin, F>
526 where
527 F: Num,
528 YOrigin: OriginGenerator<Y, F>,
529 Y: Space + AXPY<F>,
530 X: Space + Instance<X> + HasDual<F>,
531 XDOrigin: OriginGenerator<X::DualSpace, F> + Clone,
532 Yprime: Space + Instance<Yprime>,
533 {
534 type AdjointCodomain = X::DualSpace;
535 type Adjoint<'b>
536 = ZeroOp<Yprime, X::DualSpace, XDOrigin, NotAnOriginGenerator, F>
537 where
538 Self: 'b;
539 // () means not (pre)adjointable.
540
541 fn adjoint(&self) -> Self::Adjoint<'_> {
542 ZeroOp::new(self.xd_origin.clone(), NotAnOriginGenerator)
543 }
544 }
545 */
546
547 /*
548 impl<X, Y, Ypre, Xpre, F, YOrigin, XDOrigin> Preadjointable<X, Ypre>
549 for ZeroOp<X, Y, YOrigin, XDOrigin, F>
550 where
551 F: Num,
552 YOrigin: OriginGenerator<Y, F>,
553 Y: Space + AXPY<F>,
554 X: Space + Instance<X> + HasDual<F>,
555 XDOrigin: OriginGenerator<Xpre, F> + Clone,
556 Ypre: Space + Instance<Ypre> + HasDual<DualSpace = X>,
557 {
558 type PreadjointCodomain = Xpre;
559 type Preadjoint<'b>
560 = ZeroOp<Ypre, Xpre, XDOrigin, NotAnOriginGenerator, F>
561 where
562 Self: 'b;
563 // () means not (pre)adjointable.
564
565 fn adjoint(&self) -> Self::Adjoint<'_> {
566 ZeroOp::new(self.xpre_origin.clone(), NotAnOriginGenerator)
567 }
568 }
569 */
570 476
571 impl<S, T, E, X> Linear<X> for Composition<S, T, E> 477 impl<S, T, E, X> Linear<X> for Composition<S, T, E>
572 where 478 where
573 X: Space, 479 X: Space,
574 T: Linear<X>, 480 T: Linear<X>,

Mercurial Repository: alg_tools

mercurial