From 062c9d3765ed51183b70d28a01420facc8ac899c Mon Sep 17 00:00:00 2001 From: Anders Olsson Date: Thu, 26 Oct 2023 11:09:30 +0200 Subject: [PATCH] Added quality function --- src/history.rs | 62 ++++++++++++++ src/lib.rs | 2 +- src/matrix.rs | 213 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 276 insertions(+), 1 deletion(-) create mode 100644 src/matrix.rs diff --git a/src/history.rs b/src/history.rs index 4a09dc4..0d26b8c 100644 --- a/src/history.rs +++ b/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::>(); + + 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)] diff --git a/src/lib.rs b/src/lib.rs index 4d9192a..76fb17b 100644 --- a/src/lib.rs +++ b/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; diff --git a/src/matrix.rs b/src/matrix.rs new file mode 100644 index 0000000..cf74061 --- /dev/null +++ b/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::>(); + + 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 + } +}