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
use super::PositiveDefiniteKernel;
use crate::{KernelAdd, KernelError, KernelMul};
use rayon::prelude::*;
use std::{f64::consts::PI, ops::Add, ops::Mul};
#[derive(Clone, Debug)]
pub struct SpectralMixture {
p: usize,
q: usize,
}
impl SpectralMixture {
pub fn new(p: usize, q: usize) -> Self {
Self { p, q }
}
}
impl PositiveDefiniteKernel<Vec<f64>> for SpectralMixture {
fn params_len(&self) -> usize {
self.q + self.p * self.q + self.p * self.q
}
fn value(&self, params: &[f64], x: &Vec<f64>, xprime: &Vec<f64>) -> Result<f64, KernelError> {
if params.len() != self.params_len() {
return Err(KernelError::ParametersLengthMismatch.into());
}
if self.p != x.len() {
return Err(KernelError::ParametersLengthMismatch.into());
}
if x.len() != xprime.len() {
return Err(KernelError::InvalidArgument.into());
}
let w = ¶ms[0..self.q];
let v = ¶ms[self.q..self.q + self.p * self.q];
let mu = ¶ms[self.q + self.p * self.q..self.q + self.p * self.q + self.p * self.q];
let fx = (0..self.q)
.into_par_iter()
.map(|q| {
w[q] * (0..self.p)
.into_par_iter()
.map(|p| {
(-2.0 * PI.powi(2) * (x[p] - xprime[p]).powi(2) * v[self.p * p + q]).exp()
* (2.0 * PI * (x[p] - xprime[p]) * mu[self.p * p + q]).cos()
})
.product::<f64>()
})
.sum();
Ok(fx)
}
}
impl<R> Add<R> for SpectralMixture
where
R: PositiveDefiniteKernel<Vec<f64>>,
{
type Output = KernelAdd<Self, R, Vec<f64>>;
fn add(self, rhs: R) -> Self::Output {
Self::Output::new(self, rhs)
}
}
impl<R> Mul<R> for SpectralMixture
where
R: PositiveDefiniteKernel<Vec<f64>>,
{
type Output = KernelMul<Self, R, Vec<f64>>;
fn mul(self, rhs: R) -> Self::Output {
Self::Output::new(self, rhs)
}
}
#[cfg(test)]
mod tests {
use crate::*;
#[test]
fn it_works() {
let kernel = SpectralMixture::new(1, 2);
let test_value = kernel.value(&[1.0, 1.0], &vec![0.0, 0.0, 0.0], &vec![0.0, 0.0, 0.0]);
match test_value {
Err(KernelError::ParametersLengthMismatch) => (),
_ => panic!(),
};
}
}