From cb988ce13b59c42493fe8543f90b18cf019bb13c Mon Sep 17 00:00:00 2001 From: Dmitry Ustalov Date: Sat, 15 Jun 2024 18:42:59 +0200 Subject: [PATCH] Add translation from Python --- Cargo.toml | 15 ++++ pyproject.toml | 16 ++++ src/lib.rs | 198 +++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 229 insertions(+) create mode 100644 pyproject.toml create mode 100644 src/lib.rs diff --git a/Cargo.toml b/Cargo.toml index 34325a2..fcd7c69 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -4,4 +4,19 @@ version = "0.0.1" edition = "2021" license = "Apache-2.0" +[lib] +name = "evalica" +crate-type = ["cdylib"] + [dependencies] +pyo3 = "0.21.1" +ndarray = "0.15.6" +ndarray-linalg = "0.16.0" +num-traits = "0.2.19" +approx = "0.5.1" +rand = "0.8.5" +numpy = "0.21.0" + +[features] +extension-module = ["pyo3/extension-module"] +default = ["extension-module"] diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..6db5b84 --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,16 @@ +[build-system] +requires = ["maturin>=1.6,<2.0"] +build-backend = "maturin" + +[project] +name = "evalica" +requires-python = ">= 3.8" +classifiers = [ + "Programming Language :: Rust", + "Programming Language :: Python :: Implementation :: CPython", + "Programming Language :: Python :: Implementation :: PyPy", +] +dynamic = ["version"] + +[tool.maturin] +features = ["pyo3/extension-module"] diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..eb2358d --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,198 @@ +use numpy::PyArrayMethods; +use pyo3::prelude::*; + +use ndarray::prelude::*; +use ndarray::Array2; +use ndarray_linalg::Norm; +use numpy::{IntoPyArray, PyArray1, PyArray2}; +use rand::prelude::*; + +const EPS: f64 = 1e-8; + +pub fn bradley_terry(m: &Array2) -> (Array1, usize) { + let t = m.t().to_owned() + m; + + let active = t.mapv(|x| x > 0); + + let w: Array1 = m.sum_axis(Axis(1)); + + let mut z: Array2 = Array2::zeros(m.raw_dim()); + + let mut p: Array1 = Array1::ones(m.shape()[0]); + let mut p_new: Array1 = p.clone(); + + let mut converged = false; + let mut iterations = 0; + + while !converged { + iterations += 1; + + for ((i, j), &active_val) in active.indexed_iter() { + if active_val { + z[[i, j]] = t[[i, j]] as f64 / (p[i] + p[j]); + } + } + + p_new.fill(0.0); + + for i in 0..m.shape()[0] { + p_new[i] = w[i] as f64 / z.column(i).sum(); + } + + p_new /= p_new.sum(); + + let diff_norm = (&p_new - &p).norm(); + + converged = diff_norm < EPS; + + p.assign(&p_new); + } + + (p, iterations) +} + +fn compute_ties_and_wins(m: &Array2) -> (Array2, Array2) { + let mut t = m.clone(); + for ((i, j), t) in t.indexed_iter_mut() { + *t = std::cmp::min(m[[i, j]], m[[j, i]]); + } + let w = m - &t; + (t, w) +} + +pub fn newman(m: &Array2, seed: u64, tolerance: f64, limit: usize) -> (Array1, usize) { + let (t, w) = compute_ties_and_wins(m); + + let mut rng = StdRng::seed_from_u64(seed); + + let mut pi: Array1 = Array1::from_shape_fn(m.shape()[0], |_| rng.gen_range(0.0..1.0)); + let mut v: f64 = rng.gen_range(0.0..1.0); + + let mut converged = false; + let mut iterations = 0; + + while !converged && iterations < limit { + iterations += 1; + + let pi_broadcast = pi.broadcast((pi.len(), pi.len())).unwrap().to_owned(); + let pi_broadcast_t = pi_broadcast.t().to_owned(); + let pi_sum = &pi_broadcast + &pi_broadcast_t; + let sqrt_pi_product = (pi_broadcast.clone() * pi_broadcast_t.clone()).mapv(f64::sqrt); + + let denominator_common = &pi_sum + 2.0 * v * &sqrt_pi_product; + + let v_numerator = + (&t.mapv(|x| x as f64) * (&pi_broadcast + &pi_broadcast_t) / &denominator_common).sum() + / 2.0; + + let v_denominator = + (&w.mapv(|x| x as f64) * (2.0 * &sqrt_pi_product) / &denominator_common).sum(); + + v = v_numerator / v_denominator; + + if v.is_nan() { + v = tolerance; + } + + let pi_old = pi.clone(); + + let pi_numerator = ((w.mapv(|x| x as f64) + t.mapv(|x| x as f64) / 2.0) + * (&pi_broadcast + v * &sqrt_pi_product) + / (&pi_sum + 2.0 + v * &sqrt_pi_product)) + .sum_axis(Axis(1)); + + let pi_denominator = ((w.mapv(|x| x as f64) + t.mapv(|x| x as f64) / 2.0) + * (1.0 + v * &sqrt_pi_product) + / (&pi_sum + 2.0 + v * &sqrt_pi_product)) + .sum_axis(Axis(0)); + + pi = &pi_numerator / &pi_denominator; + + pi.iter_mut().for_each(|x| { + if x.is_nan() { + *x = tolerance; + } + }); + + converged = pi.iter().zip(pi_old.iter()).all(|(p_new, p_old)| { + (p_new / (p_new + 1.0) - p_old / (p_old + 1.0)).abs() <= tolerance + }); + } + + (pi, iterations) +} + +#[cfg(test)] +mod tests { + use super::*; + use ndarray::array; + + #[test] + fn test_bradley_terry() { + let m: Array2 = array![ + [0, 1, 2, 0, 1], + [2, 0, 2, 1, 0], + [1, 2, 0, 0, 1], + [1, 2, 1, 0, 2], + [2, 0, 1, 3, 0] + ]; + + let (p, iterations) = bradley_terry(&m); + + assert_eq!(p.len(), m.shape()[0]); + assert_ne!(iterations, 0); + + let expected_p = array![0.12151104, 0.15699947, 0.11594851, 0.31022851, 0.29531247]; + + for (a, b) in p.iter().zip(expected_p.iter()) { + assert!((a - b).abs() < EPS); + } + } + + #[test] + fn test_newman() { + let m = array![ + [0, 1, 2, 0, 1], + [2, 0, 2, 1, 0], + [1, 2, 0, 0, 1], + [1, 2, 1, 0, 2], + [2, 0, 1, 3, 0] + ]; + + let seed = 0; + let tolerance = 1e-6; + let limit = 100; + + let (pi, iterations) = newman(&m, seed, tolerance, limit); + + assert_eq!(pi.len(), m.shape()[0]); + assert_ne!(iterations, 0); + } +} + +#[pyfunction] +fn py_bradley_terry(py: Python, m: &Bound>) -> PyResult<(Py>, usize)> { + let m = unsafe { m.as_array().to_owned() }; + let (pi, iterations) = bradley_terry(&m); + Ok((pi.into_pyarray_bound(py).unbind(), iterations)) +} + +#[pyfunction] +fn py_newman( + py: Python, + m: &Bound>, + seed: u64, + tolerance: f64, + limit: usize, +) -> PyResult<(Py>, usize)> { + let m = unsafe { m.as_array().to_owned() }; + let (pi, iterations) = newman(&m, seed, tolerance, limit); + Ok((pi.into_pyarray_bound(py).unbind(), iterations)) +} + +#[pymodule] +fn evalica(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(py_bradley_terry, m)?)?; + m.add_function(wrap_pyfunction!(py_newman, m)?)?; + Ok(()) +}