diff --git a/Cargo.toml b/Cargo.toml index 44fa244..5fc43b6 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -7,10 +7,11 @@ license = "MIT OR Apache-2.0" [dependencies] lapacke = "0.5" -ndarray = { version = "0.15", features = ["approx"] } +ndarray = { version = "0.15", features = ["approx", "blas"] } [dev-dependencies] approx = "0.4" -intel-mkl-src = "0.5" rand = "0.8" time = "0.2" +blas-src = { version = "0.8", features = ["openblas"] } +openblas-src = { version = "0.10", features = ["cblas", "system"] } diff --git a/examples/abcdef.rs b/examples/abcdef.rs index cce7778..01726a0 100644 --- a/examples/abcdef.rs +++ b/examples/abcdef.rs @@ -1,4 +1,4 @@ -extern crate intel_mkl_src; +extern crate blas_src; use kickscore as ks; diff --git a/examples/kickscore-basics.rs b/examples/kickscore-basics.rs index 32c84dd..e3d39a5 100644 --- a/examples/kickscore-basics.rs +++ b/examples/kickscore-basics.rs @@ -1,4 +1,4 @@ -extern crate intel_mkl_src; +extern crate blas_src; use kickscore as ks; diff --git a/examples/nba-history.rs b/examples/nba-history.rs index ba023ee..ff005e4 100644 --- a/examples/nba-history.rs +++ b/examples/nba-history.rs @@ -1,4 +1,4 @@ -extern crate intel_mkl_src; +extern crate blas_src; use std::cmp::Ordering; use std::collections::HashSet; diff --git a/src/kernel.rs b/src/kernel.rs index 29ebbaf..544d217 100644 --- a/src/kernel.rs +++ b/src/kernel.rs @@ -217,7 +217,7 @@ impl Kernel for Vec> { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/kernel/affine.rs b/src/kernel/affine.rs index 94f1c6c..258b093 100644 --- a/src/kernel/affine.rs +++ b/src/kernel/affine.rs @@ -77,7 +77,7 @@ impl Kernel for Affine { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/kernel/constant.rs b/src/kernel/constant.rs index b06111d..1340e1f 100644 --- a/src/kernel/constant.rs +++ b/src/kernel/constant.rs @@ -59,7 +59,7 @@ impl Kernel for Constant { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/kernel/exponential.rs b/src/kernel/exponential.rs index 9bd05c5..cba50b0 100644 --- a/src/kernel/exponential.rs +++ b/src/kernel/exponential.rs @@ -62,7 +62,7 @@ impl Kernel for Exponential { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/kernel/matern32.rs b/src/kernel/matern32.rs index bf3d468..cf63391 100644 --- a/src/kernel/matern32.rs +++ b/src/kernel/matern32.rs @@ -89,7 +89,7 @@ impl Kernel for Matern32 { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/kernel/matern52.rs b/src/kernel/matern52.rs index 280def1..d4533e0 100644 --- a/src/kernel/matern52.rs +++ b/src/kernel/matern52.rs @@ -117,7 +117,7 @@ impl Kernel for Matern52 { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/kernel/wiener.rs b/src/kernel/wiener.rs index 2d84abc..1a86d0a 100644 --- a/src/kernel/wiener.rs +++ b/src/kernel/wiener.rs @@ -68,7 +68,7 @@ impl Kernel for Wiener { #[cfg(test)] mod tests { - extern crate intel_mkl_src; + extern crate blas_src; use approx::assert_abs_diff_eq; use rand::{distributions::Standard, thread_rng, Rng}; diff --git a/src/observation/ordinal.rs b/src/observation/ordinal.rs index 9301397..9d15bc2 100644 --- a/src/observation/ordinal.rs +++ b/src/observation/ordinal.rs @@ -1,7 +1,9 @@ -use std::f64::consts::TAU; +use std::f64::consts::{SQRT_2, TAU}; use crate::storage::Storage; -use crate::utils::{logphi, normcdf, normpdf, ROOTS_HERMITENORM_WS, ROOTS_HERMITENORM_XS}; +use crate::utils::{ + logphi, logsumexp2, normcdf, normpdf, ROOTS_HERMITENORM_WS, ROOTS_HERMITENORM_XS, +}; use super::{f_params, Core, Observation}; @@ -52,6 +54,66 @@ fn ll_probit_tie(x: f64, margin: f64) -> f64 { } } +fn lambdas() -> [f64; 5] { + [ + 0.44 * SQRT_2, + 0.41 * SQRT_2, + 0.40 * SQRT_2, + 0.39 * SQRT_2, + 0.36 * SQRT_2, + ] +} + +const CS: [f64; 5] = [ + 1.146480988574439e+02, + -1.508871030070582e+03, + 2.676085036831241e+03, + -1.356294962039222e+03, + 7.543285642111850e+01, +]; + +fn mm_logit_win(mean_cav: f64, cov_cav: f64) { + let mut arr1 = [0.0; 5]; + let mut arr2 = [0.0; 5]; + let mut arr3 = [0.0; 5]; + + for (i, x) in lambdas().iter().enumerate() { + let (a, b, c) = mm_probit_win(x * mean_cav, x * x * cov_cav); + + arr1[i] = a; + arr2[i] = b; + arr3[i] = c; + } + + let logpart1 = logsumexp2(arr1, CS); + + /* + dlogpart1 = (np.dot(np.exp(arr1) * arr2, CS * LAMBDAS) + / np.dot(np.exp(arr1), CS)) + d2logpart1 = (np.dot(np.exp(arr1) * (arr2 * arr2 + arr3), + CS * LAMBDAS * LAMBDAS) + / np.dot(np.exp(arr1), CS)) - (dlogpart1 * dlogpart1) + */ +} + +fn mm_logit_tie(x: f64, margin: f64) { + // +} + +fn ll_logit_win(x: f64, margin: f64) -> f64 { + let z = x - margin; + + if z > 0.0 { + -((-z).exp().ln_1p()) + } else { + z - z.exp().ln_1p() + } +} + +fn ll_logit_tie(x: f64, margin: f64) -> f64 { + ll_logit_win(x, margin) + ll_logit_win(-x, margin) + (2.0 * margin).exp_m1().ln() +} + fn cvi_expectations(mean: f64, var: f64, ll_fct: F) -> (f64, f64, f64) where F: Fn(f64) -> f64, @@ -153,16 +215,20 @@ impl Observation for ProbitTieObservation { pub struct LogitWinObservation { core: Core, - _margin: f64, + margin: f64, } impl LogitWinObservation { pub fn new(storage: &mut Storage, elems: &[(usize, f64)], t: f64, margin: f64) -> Self { Self { core: Core::new(storage, elems, t), - _margin: margin, + margin, } } + + pub fn probability(_storage: &Storage, _elems: &[(usize, f64)], _t: f64, _margin: f64) -> f64 { + todo!() + } } impl Observation for LogitWinObservation { @@ -171,22 +237,30 @@ impl Observation for LogitWinObservation { } fn kl_update(&mut self, storage: &mut Storage, lr: f64) -> f64 { - self.core.kl_update(storage, lr, |mean, var| todo!()) + let margin = self.margin; + + self.core.kl_update(storage, lr, |mean, var| { + cvi_expectations(mean, var, |x| ll_logit_win(x, margin)) + }) } } pub struct LogitTieObservation { core: Core, - _margin: f64, + margin: f64, } impl LogitTieObservation { pub fn new(storage: &mut Storage, elems: &[(usize, f64)], t: f64, margin: f64) -> Self { Self { core: Core::new(storage, elems, t), - _margin: margin, + margin, } } + + pub fn probability(_storage: &Storage, _elems: &[(usize, f64)], _t: f64, _margin: f64) -> f64 { + todo!() + } } impl Observation for LogitTieObservation { @@ -195,7 +269,11 @@ impl Observation for LogitTieObservation { } fn kl_update(&mut self, storage: &mut Storage, lr: f64) -> f64 { - self.core.kl_update(storage, lr, |mean, var| todo!()) + let margin = self.margin; + + self.core.kl_update(storage, lr, |mean, var| { + cvi_expectations(mean, var, |x| ll_logit_tie(x, margin)) + }) } } diff --git a/src/utils.rs b/src/utils.rs index 9b4ee72..719e4b9 100644 --- a/src/utils.rs +++ b/src/utils.rs @@ -1,3 +1,4 @@ +use std::cmp::Ordering; use std::f64::consts::{PI, SQRT_2, TAU}; use crate::math::erfc; @@ -88,6 +89,17 @@ pub fn logphi(z: f64) -> (f64, f64) { } } +pub fn logsumexp2(xs: [f64; 5], bs: [f64; 5]) -> f64 { + let a = xs.iter().max_by(|a, b| a.partial_cmp(b).unwrap()).unwrap(); + + a + bs + .iter() + .zip(xs.iter().map(|x| (x - a).exp())) + .map(|(b, x)| b * x) + .sum::() + .ln() +} + pub const ROOTS_HERMITENORM_XS: [f64; 30] = [ -9.706236, -8.68083772, diff --git a/tests/binary-1.rs b/tests/binary-1.rs index 0e5dd1e..ec84cc9 100644 --- a/tests/binary-1.rs +++ b/tests/binary-1.rs @@ -1,4 +1,4 @@ -extern crate intel_mkl_src; +extern crate blas_src; use kickscore as ks; diff --git a/tests/kickscore-basics.rs b/tests/kickscore-basics.rs index 8b6a1a2..2836b18 100644 --- a/tests/kickscore-basics.rs +++ b/tests/kickscore-basics.rs @@ -1,4 +1,4 @@ -extern crate intel_mkl_src; +extern crate blas_src; use approx::assert_abs_diff_eq; use kickscore as ks; @@ -41,5 +41,5 @@ fn kickscore_basic() { let (p_win, _p_los) = model.probabilities(&[&"Jerry"], &[&"Tom"], -1.0); - assert_abs_diff_eq!(p_win, 0.903756079972532, epsilon = f64::EPSILON); + assert_abs_diff_eq!(p_win, 0.9037560799725326, epsilon = f64::EPSILON); } diff --git a/tests/nba-history.rs b/tests/nba-history.rs index f21c57c..f32a126 100644 --- a/tests/nba-history.rs +++ b/tests/nba-history.rs @@ -1,4 +1,4 @@ -extern crate intel_mkl_src; +extern crate blas_src; use std::cmp::Ordering; use std::collections::HashSet;