Added quality function

src/history.rs

@@ -4,6 +4,7 @@ use crate::{
     agent::{self, Agent},
     batch::{self, Batch},
     gaussian::Gaussian,
+    matrix::Matrix,
     player::Player,
     sort_time, tuple_gt, tuple_max, Index, BETA, GAMMA, MU, P_DRAW, SIGMA,
 };
@@ -417,6 +418,67 @@ impl History {
 
         self.size += n;
     }
+
+    pub fn quality(&self, rating_groups: &[&[Gaussian]]) -> f64 {
+        let flatten_ratings = rating_groups
+            .iter()
+            .flat_map(|group| group.iter())
+            .collect::<Vec<_>>();
+
+        let flatten_weights = vec![1.0; flatten_ratings.len()].into_boxed_slice();
+
+        let length = flatten_ratings.len();
+
+        let mut mean_matrix = Matrix::new(length, 1);
+
+        for (i, rating) in flatten_ratings.iter().enumerate() {
+            mean_matrix[(i, 0)] = rating.mu;
+        }
+
+        let mut variance_matrix = Matrix::new(length, length);
+
+        for (i, rating) in flatten_ratings.iter().enumerate() {
+            variance_matrix[(i, i)] = rating.sigma.powi(2);
+        }
+
+        let mut rotated_a_matrix = Matrix::new(rating_groups.len() - 1, length);
+
+        let mut t = 0;
+        let mut x = 0;
+
+        for (row, group) in rating_groups.windows(2).enumerate() {
+            let current = group[0];
+            let next = group[1];
+
+            for n in t..t + current.len() {
+                rotated_a_matrix[(row, n)] = flatten_weights[n];
+
+                x += 1;
+            }
+
+            t += current.len();
+
+            for n in x..x + next.len() {
+                rotated_a_matrix[(row, n)] = -flatten_weights[n];
+            }
+
+            x += next.len();
+        }
+
+        let a_matrix = rotated_a_matrix.transpose();
+
+        let ata = self.beta.powi(2) * &rotated_a_matrix * &a_matrix;
+        let atsa = &rotated_a_matrix * &variance_matrix * &a_matrix;
+
+        let start = mean_matrix.transpose() * &a_matrix;
+        let middle = &ata + &atsa;
+        let end = &rotated_a_matrix * &mean_matrix;
+
+        let e_arg = (-0.5 * &start * &middle.inverse() * &end).determinant();
+        let s_arg = ata.determinant() / middle.determinant();
+
+        e_arg.exp() * s_arg.sqrt()
+    }
 }
 
 #[cfg(test)]

src/lib.rs

@@ -10,8 +10,8 @@ mod approx;
 pub mod batch;
 mod game;
 pub mod gaussian;
-// mod gaussian2;
 mod history;
+mod matrix;
 mod message;
 pub mod player;
 

src/matrix.rs

@@ -0,0 +1,213 @@
+use std::ops;
+
+fn det(m: &[f64], x: usize) -> f64 {
+    if x == 1 {
+        m[0]
+    } else if x == 2 {
+        m[0] * m[3] - m[1] * m[2]
+    } else {
+        let mut d = 0.0;
+
+        for n in 0..x {
+            let ms = m
+                .iter()
+                .enumerate()
+                .skip(x)
+                .filter(|(i, _)| (i % x) != n)
+                .map(|(_, v)| *v)
+                .collect::<Vec<_>>();
+
+            d += (-1.0f64).powi(n as i32) * m[n] * det(&ms, x - 1);
+        }
+
+        d
+    }
+}
+
+#[derive(Clone, Debug)]
+pub struct Matrix {
+    data: Box<[f64]>,
+    height: usize,
+    width: usize,
+}
+
+impl Matrix {
+    pub fn new(height: usize, width: usize) -> Matrix {
+        Matrix {
+            data: vec![0.0; height * width].into_boxed_slice(),
+            height,
+            width,
+        }
+    }
+
+    pub fn transpose(&self) -> Matrix {
+        let mut matrix = Matrix::new(self.width, self.height);
+
+        for c in 0..self.width {
+            for r in 0..self.height {
+                matrix[(c, r)] = self[(r, c)];
+            }
+        }
+
+        matrix
+    }
+
+    pub fn minor(&self, row_n: usize, col_n: usize) -> Matrix {
+        let mut matrix = Matrix::new(self.height - 1, self.width - 1);
+
+        let mut nr = 0;
+
+        for r in 0..self.height {
+            if r == row_n {
+                continue;
+            }
+
+            let mut nc = 0;
+
+            for c in 0..self.width {
+                if c == col_n {
+                    continue;
+                }
+
+                matrix[(nr, nc)] = self[(r, c)];
+
+                nc += 1;
+            }
+
+            nr += 1;
+        }
+
+        matrix
+    }
+
+    pub fn determinant(&self) -> f64 {
+        debug_assert!(self.width == self.height);
+
+        det(&self.data, self.width)
+    }
+
+    pub fn adjugate(&self) -> Matrix {
+        debug_assert!(self.width == self.height);
+
+        let mut matrix = Matrix::new(self.height, self.width);
+
+        if matrix.height == 2 {
+            matrix[(0, 0)] = self[(1, 1)];
+            matrix[(0, 1)] = -self[(0, 1)];
+            matrix[(1, 0)] = -self[(1, 0)];
+            matrix[(1, 1)] = self[(0, 0)];
+        } else {
+            for r in 0..matrix.height {
+                for c in 0..matrix.width {
+                    let sign = if (r + c) % 2 == 0 { 1.0 } else { -1.0 };
+
+                    matrix[(r, c)] = self.minor(r, c).determinant() * sign;
+                }
+            }
+        }
+
+        matrix
+    }
+
+    pub fn inverse(&self) -> Matrix {
+        let mut matrix = Matrix::new(self.width, self.height);
+
+        if self.height == self.width && self.height == 1 {
+            matrix[(0, 0)] = 1.0 / self[(0, 0)];
+        } else {
+            panic!("eh, okey")
+        }
+
+        matrix
+    }
+}
+
+impl ops::Index<(usize, usize)> for Matrix {
+    type Output = f64;
+
+    fn index(&self, pos: (usize, usize)) -> &Self::Output {
+        &self.data[(self.width * pos.0) + pos.1]
+    }
+}
+
+impl ops::IndexMut<(usize, usize)> for Matrix {
+    fn index_mut(&mut self, pos: (usize, usize)) -> &mut Self::Output {
+        &mut self.data[(self.width * pos.0) + pos.1]
+    }
+}
+
+impl<'a> ops::Mul<&'a Matrix> for f64 {
+    type Output = Matrix;
+
+    fn mul(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(rhs.height, rhs.width);
+
+        for r in 0..rhs.height {
+            for c in 0..rhs.width {
+                matrix[(r, c)] = self * rhs[(r, c)];
+            }
+        }
+
+        matrix
+    }
+}
+
+impl<'a> ops::Mul<&'a Matrix> for Matrix {
+    type Output = Matrix;
+
+    fn mul(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(self.height, rhs.width);
+
+        for r in 0..matrix.height {
+            for c in 0..matrix.width {
+                let mut value = 0.0;
+
+                for x in 0..self.width {
+                    value += self[(r, x)] * rhs[(x, c)];
+                }
+
+                matrix[(r, c)] = value;
+            }
+        }
+
+        matrix
+    }
+}
+
+impl<'a> ops::Mul<&'a Matrix> for &'a Matrix {
+    type Output = Matrix;
+
+    fn mul(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(self.height, rhs.width);
+
+        for r in 0..matrix.height {
+            for c in 0..matrix.width {
+                let mut value = 0.0;
+
+                for x in 0..self.width {
+                    value += self[(r, x)] * rhs[(x, c)];
+                }
+
+                matrix[(r, c)] = value;
+            }
+        }
+
+        matrix
+    }
+}
+
+impl<'a> ops::Add<&'a Matrix> for &'a Matrix {
+    type Output = Matrix;
+
+    fn add(self, rhs: &'a Matrix) -> Matrix {
+        let mut matrix = Matrix::new(self.height, self.width);
+
+        for r in 0..matrix.height {
+            for c in 0..matrix.width {
+                matrix[(r, c)] = self[(r, c)] + rhs[(r, c)];
+            }
+        }
+
+        matrix
+    }
+}