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, r_assert_relative_equal_result,
};
use strafe_trait::{Model, ModelBuilder};
use strafe_type::ModelMatrix;

use crate::regression::glm::{
    family::{binomial::Binomial, core::Family},
    link::core::Link,
    GeneralizedLinearRegressionBuilder,
};

pub fn family_option_string(family_option: &str) -> String {
    match family_option {
        "Binomial" => "binomial()".to_string(),
        "Gamma" => "Gamma()".to_string(),
        "Gaussian" => "gaussian()".to_string(),
        "InverseGaussian" => "inverse.gaussian()".to_string(),
        "Poisson" => "poisson()".to_string(),
        "Quasi(Constant)" => "quasi(variance = 'constant')".to_string(),
        "Quasi(Mu)" => "quasi(variance = 'mu')".to_string(),
        "Quasi(MuSqrd)" => "quasi(variance = 'mu^2')".to_string(),
        "Quasi(MuCubed)" => "quasi(variance = 'mu^3')".to_string(),
        "Quasi(MuOneMinusMu)" => "quasi(variance = 'mu(1-mu)')".to_string(),
        "QuasiBinomial" => "quasibinomial()".to_string(),
        "QuasiPoisson" => "quasipoisson()".to_string(),
        _ => unimplemented!(),
    }
}

pub fn binomial_logit_data(
    seed: u16,
    num_features: f64,
    num_rows: f64,
    fuzz_amount: f64,
) -> (DMatrix<f64>, DMatrix<f64>, DMatrix<f64>) {
    let num_features = num_features.to_usize().unwrap();
    let num_rows = num_rows.to_usize().unwrap();
    let fuzz_amount = fuzz_amount.to_usize().unwrap();

    let mut features = RGenerator::new()
        .set_seed(seed)
        .generate_uniform(num_features + 1, (-100, 100));
    features.iter_mut().for_each(|f| {
        if f.is_nan() {
            *f = 1.0
        }
    });

    let mut x = RGenerator::new()
        .set_seed(seed + 1)
        .generate_uniform_matrix(num_rows, num_features, (-10, 10));

    let mut weights = RGenerator::new()
        .set_seed(seed + 2)
        .generate_uniform_matrix(num_rows, 1, (0, 10));

    let mut errors =
        RGenerator::new()
            .set_seed(seed + 2)
            .generate_errors(num_rows, (0, 1), fuzz_amount);

    let pi = RTester::new()
        .set_seed(seed + 4)
        .set_display(&RString::from_string(format!(
            "exp(as.vector({0}%*%t(cbind(1,{1})))+{2})/(1+exp(as.vector({0}%*%t(cbind(1,{1})))+{2}))",
            RString::from_f64_slice(&features),
            RString::from_f64_matrix(&x),
            RString::from_f64_slice(&errors)
        )))
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R array error");

    let mut y = RTester::new()
        .set_seed(seed + 5)
        .set_display(&RString::from_string(format!(
            "matrix(rbinom(10000000, 1, prob={0}), nrow={1}, ncol={2})",
            RString::from_f64_slice(&pi),
            num_rows,
            1,
        )))
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("R array error");
    y.iter_mut().for_each(|y_i| {
        if y_i.is_nan() {
            *y_i = 1.0
        }
    });

    (x, y, weights)
}

pub fn fit_models<L: 'static + Link + Clone, F: 'static + Family<L> + Clone>(
    x: DMatrix<f64>,
    y: DMatrix<f64>,
    weights: DMatrix<f64>,
    family: F,
    family_string: &str,
) {
    let mut model = GeneralizedLinearRegressionBuilder::new(family)
        .with_x(&ModelMatrix::from(&x))
        .with_y(&ModelMatrix::from(&y))
        .with_weights(&ModelMatrix::from(&weights))
        .build();
    // let model_test = model.test(&()).unwrap();

    let r_regression = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};weights={};family={};",
            RString::from_f64_matrix(&y),
            RString::from_f64_matrix(&x),
            RString::from_f64_matrix(&weights),
            RString::from_str(family_string).unwrap(),
        )))
        .set_display(&RString::from_str("summary(glm(y~x,weights=weights,family=family))").unwrap())
        .run()
        .expect("R error")
        .as_glm_regression();

    for i in 0..model.parameters().unwrap().len() {
        // Estimates
        assert!(
            r_regression.coefs[i].estimate
                >= model.parameters().unwrap()[i].confidence_interval().0
        );
        assert!(
            r_regression.coefs[i].estimate
                <= model.parameters().unwrap()[i].confidence_interval().1
        );
        r_assert_relative_equal_result!(
            Ok(r_regression.coefs[i].estimate),
            &model.parameters().unwrap()[i].estimate(),
            1e-4
        );
        r_assert_relative_equal_result!(
            Ok(r_regression.coefs[i].estimate),
            &model.parameters().unwrap()[i].estimate(),
            1e-4
        );

        //     // Significance Test
        //     r_assert_relative_equal_result!(
        //         Ok(r_regression.coefs[i].z.precision_round_to(4)),
        //         model_test.significant_coef_tests[i]
        //             .statistic()
        //             .precision_round_to(4),
        //         1e-4
        //     );
        //
        //     if r_regression.coefs[i].p <= 0.05 {
        //         if model_test.significant_coef_tests[i].probability_value() > 0.05 {
        //             r_assert_relative_equal_result!(
        //                 Ok(r_regression.coefs[i].p),
        //                 model_test.significant_coef_tests[i].probability_value(),
        //                 1e-4
        //             );
        //         }
        //     } else {
        //         if model_test.significant_coef_tests[i].probability_value() <= 0.05 {
        //             r_assert_relative_equal_result!(
        //                 Ok(r_regression.coefs[i].p),
        //                 model_test.significant_coef_tests[i].probability_value(),
        //                 1e-4
        //             );
        //         }
        //     }
    }

    r_assert_relative_equal_result!(Ok(r_regression.null_dev), &model.null_deviance, 1e-4);
    r_assert_relative_equal_result!(Ok(r_regression.resid_dev), &model.residual_deviance, 1e-4);
    r_assert_relative_equal_result!(Ok(r_regression.aic), &model.aic, 1e-4);

    let r_res = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};weights={};family={};",
            RString::from_f64_matrix(&y),
            RString::from_f64_matrix(&x),
            RString::from_f64_matrix(&weights),
            RString::from_str(family_string).unwrap(),
        )))
        .set_display(
            &RString::from_str("unname(glm(y~x,weights=weights,family=family)$residuals)").unwrap(),
        )
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model.residuals().unwrap().matrix().iter().zip(r_res.iter()) {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-4);
    }

    let r_pred = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};weights={};family={};",
            RString::from_f64_matrix(&y),
            RString::from_f64_matrix(&x),
            RString::from_f64_matrix(&weights),
            RString::from_str(family_string).unwrap(),
        )))
        .set_display(
            &RString::from_str("unname(glm(y~x,weights=weights,family=family)$fitted.values)")
                .unwrap(),
        )
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model
        .predictions()
        .unwrap()
        .estimate()
        .matrix()
        .iter()
        .zip(r_pred.iter())
    {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-4);
    }
    let r_vcov = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};weights={};family={};",
            RString::from_f64_matrix(&y),
            RString::from_f64_matrix(&x),
            RString::from_f64_matrix(&weights),
            RString::from_str(family_string).unwrap(),
        )))
        .set_display(&RString::from_str("vcov(glm(y~x,weights=weights,family=family))").unwrap())
        .run()
        .expect("R error")
        .as_f64_matrix()
        .expect("R error");

    for (&rust_ret, &r_ret) in model.variance().unwrap().iter().zip(r_vcov.iter()) {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-4);
    }

    let r_res = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};weights={};family={};",
            RString::from_f64_matrix(&y),
            RString::from_f64_matrix(&x),
            RString::from_f64_matrix(&weights),
            RString::from_str(family_string).unwrap(),
        )))
        .set_display(
            &RString::from_str("unname(rstudent(glm(y~x,weights=weights,family=family)))").unwrap(),
        )
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model
        .studentized_residuals()
        .unwrap()
        .matrix()
        .iter()
        .zip(r_res.iter())
    {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-2);
    }

    let r_res = RTester::new()
        .set_script(&RString::from_string(format!(
            "y={};x={};weights={};family={};",
            RString::from_f64_matrix(&y),
            RString::from_f64_matrix(&x),
            RString::from_f64_matrix(&weights),
            RString::from_str(family_string).unwrap(),
        )))
        .set_display(
            &RString::from_str("unname(rstandard(glm(y~x,weights=weights,family=family)))")
                .unwrap(),
        )
        .run()
        .expect("R error")
        .as_f64_vec()
        .expect("R error");

    for (&rust_ret, &r_ret) in model
        .standardized_residuals()
        .unwrap()
        .matrix()
        .iter()
        .zip(r_res.iter())
    {
        r_assert_relative_equal!(rust_ret, r_ret, 1e-2);
    }
}

pub fn generalized_linear_regression_test_inner(
    seed: u16,
    num_features: f64,
    num_rows: f64,
    fuzz_amount: f64,
) {
    let (x, y, weights) = binomial_logit_data(seed, num_features, num_rows, fuzz_amount);
    fit_models(
        x,
        y,
        weights,
        Binomial::default(),
        &family_option_string("Binomial"),
    );
}