dsmc 0.1.2

Digitalservo Motion Control Library
Documentation 
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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

use std::ops::{AddAssign, MulAssign};

use num_traits::Float;

pub struct GaussianProcessRegression<T> {
    pub x_sample: Vec<T>,
    pub y_sample: Vec<T>,
    pub x_max: T,
    pub x_min: T,
    kernel: fn(T, T) -> T,
    sense_variance: T,
    inv_cov: Vec<Vec<T>>,
    sample: usize,
}

#[derive(Debug)]
pub struct PredictedValue<T> {
    pub mean: T,
    pub stdev: T,
}

impl<T: Float + AddAssign + MulAssign> GaussianProcessRegression<T> {
    pub fn new(kernel: fn(T, T) -> T, sigma: T) -> Self {
        Self {
            x_sample: vec![],
            y_sample: vec![],
            x_max: T::zero(),
            x_min: T::zero(),
            kernel,
            sense_variance: sigma,
            inv_cov: vec![vec![]],
            sample: 0,
        }
    }

    pub fn add(&mut self, x: T, y: T) {
        self.x_sample.push(x);
        self.y_sample.push(y);

        if self.sample == 0 {
            self.x_max = x;
            self.x_min = x;
        } else {
            if self.x_max < x {
                self.x_max = x;
            }
            if self.x_min > x {
                self.x_min = x;
            }
        }

        self.sample += 1;
    }

    pub fn predict(&mut self, x: T) -> PredictedValue<T> {
        if self.inv_cov.len() != self.sample {
            let mut buffer: Vec<Vec<T>> = vec![vec![T::zero(); self.sample]; self.sample];
            for i in 0..self.sample {
                for j in 0..self.sample {
                    buffer[i][j] = (self.kernel)(self.x_sample[i], self.x_sample[j]);
                }
            }
            self.inv_cov = inverse(&buffer).unwrap();
        }

        let mut k: Vec<T> = vec![T::zero(); self.sample];
        for i in 0..self.sample {
            k[i] = (self.kernel)(self.x_sample[i], x);
        }

        let mut buffer1: Vec<T> = vec![T::zero(); self.sample];
        let mut buffer2: T = T::zero();

        let mut mean: T = T::zero();
        for i in 0..self.sample {
            for j in 0..self.sample {
                buffer1[i] += self.inv_cov[i][j] * self.y_sample[j];
            }
            mean += k[i] * buffer1[i];
        }

        for i in 0..self.sample {
            buffer1[i] = T::zero();
        }

        for i in 0..self.sample {
            for j in 0..self.sample {
                buffer1[i] += self.inv_cov[i][j] * k[j];
            }
            buffer2 += k[i] * buffer1[i];
        }

        let stdev: T = ((self.kernel)(x, x) - buffer2 + self.sense_variance)
            .abs()
            .sqrt();

        PredictedValue { mean, stdev }
    }
}

fn inverse<T: Float + MulAssign>(m: &Vec<Vec<T>>) -> Option<Vec<Vec<T>>> {
    let vlen: usize = m.len();

    let mut m1: Vec<Vec<T>> = m.clone();
    let mut m2: Vec<Vec<T>> = vec![vec![T::zero(); vlen]; vlen];
    for i in 0..vlen {
        m2[i][i] = T::one();
    }

    for i in 0..vlen {
        let mut max_row_option: usize = 0;
        let mut max_value_option: T = T::zero();
        for j in i..vlen {
            if m1[j][i].abs() > max_value_option {
                max_value_option = m1[j][i].abs();
                max_row_option = j;
            }
        }

        /* Error Handling */
        if max_value_option == T::zero() {
            return None;
        }

        let m1_buffer: Vec<T> = m1[i].clone();
        let m2_buffer: Vec<T> = m2[i].clone();
        m1[i] = m1[max_row_option].clone();
        m2[i] = m2[max_row_option].clone();
        m1[max_row_option] = m1_buffer;
        m2[max_row_option] = m2_buffer;

        let scaler: T = T::one() / m1[i][i];
        for j in 0..vlen {
            m1[i][j] *= scaler;
            m2[i][j] *= scaler;
        }

        for j in 0..vlen {
            if i != j {
                let scaler: T = m1[j][i];
                for k in 0..vlen {
                    m1[j][k] = m1[j][k] - scaler * m1[i][k];
                    m2[j][k] = m2[j][k] - scaler * m2[i][k];
                }
            }
        }
    }

    Some(m2)
}