1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

//! Linear operator algebra

use crate::generate::hstack;
use crate::types::*;
use ndarray::*;

/// Abstracted linear operator as an action to vector (`ArrayBase<S, Ix1>`) and matrix
/// (`ArrayBase<S, Ix2`)
pub trait LinearOperator {
    type Elem: Scalar;

    /// Apply operator out-place
    fn apply<S>(&self, a: &ArrayBase<S, Ix1>) -> Array1<S::Elem>
    where
        S: Data<Elem = Self::Elem>,
    {
        let mut a = a.to_owned();
        self.apply_mut(&mut a);
        a
    }

    /// Apply operator in-place
    fn apply_mut<S>(&self, a: &mut ArrayBase<S, Ix1>)
    where
        S: DataMut<Elem = Self::Elem>,
    {
        let b = self.apply(a);
        azip!((a in a, &b in &b) *a = b);
    }

    /// Apply operator with move
    fn apply_into<S>(&self, mut a: ArrayBase<S, Ix1>) -> ArrayBase<S, Ix1>
    where
        S: DataOwned<Elem = Self::Elem> + DataMut,
    {
        self.apply_mut(&mut a);
        a
    }

    /// Apply operator to matrix out-place
    fn apply2<S>(&self, a: &ArrayBase<S, Ix2>) -> Array2<S::Elem>
    where
        S: Data<Elem = Self::Elem>,
    {
        let cols: Vec<_> = a.axis_iter(Axis(1)).map(|col| self.apply(&col)).collect();
        hstack(&cols).unwrap()
    }

    /// Apply operator to matrix in-place
    fn apply2_mut<S>(&self, a: &mut ArrayBase<S, Ix2>)
    where
        S: DataMut<Elem = Self::Elem>,
    {
        for mut col in a.axis_iter_mut(Axis(1)) {
            self.apply_mut(&mut col)
        }
    }

    /// Apply operator to matrix with move
    fn apply2_into<S>(&self, mut a: ArrayBase<S, Ix2>) -> ArrayBase<S, Ix2>
    where
        S: DataOwned<Elem = Self::Elem> + DataMut,
    {
        self.apply2_mut(&mut a);
        a
    }
}

impl<A, Sa> LinearOperator for ArrayBase<Sa, Ix2>
where
    A: Scalar,
    Sa: Data<Elem = A>,
{
    type Elem = A;

    fn apply<S>(&self, a: &ArrayBase<S, Ix1>) -> Array1<A>
    where
        S: Data<Elem = A>,
    {
        self.dot(a)
    }
}