r2rs-stats 0.1.1

// "Whatever you do, work at it with all your heart, as working for the Lord,
// not for human masters, since you know that you will receive an inheritance
// from the Lord as a reward. It is the Lord Christ you are serving."
// (Col 3:23-24)

use std::str::FromStr;

use nalgebra::DMatrix;
use num_traits::ToPrimitive;
use strafe_testing::{
    r::{RGenerator, RString, RTester},
    r_assert_relative_equal,
};

use crate::regression::glm::link::{
    cauchit::Cauchit, cloglog::CLogLog, core::Link, identity::Identity, inverse::Inverse, log::Log,
    logit::Logit, one_div_mu_sqrd::OneDivMuSqrd, power::Power, probit::Probit, sqrt::Sqrt,
};

pub enum LinkOptions {
    Cauchit,
    CLogLog,
    Identity,
    Inverse,
    Log,
    Logit,
    OneDivMuSqrd,
    Power(f64),
    Probit,
    Sqrt,
}

impl ToString for LinkOptions {
    fn to_string(&self) -> String {
        match self {
            LinkOptions::Cauchit => "make.link('cauchit')".to_string(),
            LinkOptions::CLogLog => "make.link('cloglog')".to_string(),
            LinkOptions::Identity => "make.link('identity')".to_string(),
            LinkOptions::Inverse => "make.link('inverse')".to_string(),
            LinkOptions::Log => "make.link('log')".to_string(),
            LinkOptions::Logit => "make.link('logit')".to_string(),
            LinkOptions::OneDivMuSqrd => "make.link('1/mu^2')".to_string(),
            LinkOptions::Power(lambda) => format!("power(lambda={lambda})"),
            LinkOptions::Probit => "make.link('probit')".to_string(),
            LinkOptions::Sqrt => "make.link('sqrt')".to_string(),
        }
    }
}

impl LinkOptions {
    fn to_link(&self) -> Box<dyn Link> {
        match self {
            LinkOptions::Cauchit => Box::new(Cauchit::new()),
            LinkOptions::CLogLog => Box::new(CLogLog::new()),
            LinkOptions::Identity => Box::new(Identity::new()),
            LinkOptions::Inverse => Box::new(Inverse::new()),
            LinkOptions::Log => Box::new(Log::new()),
            LinkOptions::Logit => Box::new(Logit::new()),
            LinkOptions::OneDivMuSqrd => Box::new(OneDivMuSqrd::new()),
            LinkOptions::Power(lambda) => Box::new(Power::new().with_lambda(*lambda)),
            LinkOptions::Probit => Box::new(Probit::new()),
            LinkOptions::Sqrt => Box::new(Sqrt::new()),
        }
    }
}

pub fn link_test_inner(seed: u16, num_rows: f64, link: LinkOptions) {
    let num_rows = num_rows.to_usize().unwrap();

    let eta_vec = RGenerator::new()
        .set_seed(seed)
        .generate_uniform(num_rows, (-100, 100));
    let etas = DMatrix::from_iterator(eta_vec.len(), 1, eta_vec.iter().copied());

    let mu_vec = RGenerator::new()
        .set_seed(seed + 1)
        .generate_uniform(num_rows, (0, 1));
    let mus = DMatrix::from_iterator(mu_vec.len(), 1, mu_vec.iter().copied());

    let library = RString::from_library_name("stats");

    // Link Function
    let r_link_func = RTester::new()
        .add_library(&library)
        .set_script(&RString::from_string(format!(
            "mu={};link={};",
            RString::from_f64_matrix(&mus),
            RString::from_str(&link.to_string()).unwrap(),
        )))
        .set_display(&RString::from_str("link$linkfun(mu)").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("Could not read matrix");

    let rust_link_func = link.to_link().link_func(&mus);

    for (&r_val, &rust_val) in r_link_func.iter().zip(rust_link_func.iter()) {
        r_assert_relative_equal!(r_val, rust_val, 1e-4);
    }

    // Link Inverse
    let r_link_inverse = RTester::new()
        .add_library(&library)
        .set_script(&RString::from_string(format!(
            "eta={};link={};",
            RString::from_f64_matrix(&etas),
            RString::from_str(&link.to_string()).unwrap(),
        )))
        .set_display(&RString::from_str("link$linkinv(eta)").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("Could not read matrix");

    let rust_link_inverse = link.to_link().link_inverse(&etas);

    for (&r_val, &rust_val) in r_link_inverse.iter().zip(rust_link_inverse.iter()) {
        r_assert_relative_equal!(r_val, rust_val, 1e-4);
    }

    // Mu-Eta
    let r_mu_eta = RTester::new()
        .add_library(&library)
        .set_script(&RString::from_string(format!(
            "mu={};link={};",
            RString::from_f64_matrix(&mus),
            RString::from_str(&link.to_string()).unwrap(),
        )))
        .set_display(&RString::from_str("link$mu.eta(mu)").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("Could not read matrix");

    let rust_mu_eta = link.to_link().mu_eta(&mus);

    for (&r_val, &rust_val) in r_mu_eta.iter().zip(rust_mu_eta.iter()) {
        r_assert_relative_equal!(r_val, rust_val, 1e-4);
    }

    // Valid Eta
    let r_mu_eta = RTester::new()
        .add_library(&library)
        .set_script(&RString::from_string(format!(
            "eta={};link={};",
            RString::from_f64_matrix(&etas),
            RString::from_str(&link.to_string()).unwrap(),
        )))
        .set_display(&RString::from_str("link$valideta(eta)").unwrap())
        .run()
        .expect("R error")
        .as_bool()
        .expect("Could not read bool");

    let rust_mu_eta = link.to_link().valid_eta(&etas);

    assert_eq!(r_mu_eta, rust_mu_eta);
}