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 r2rs_nmath::{rng::MarsagliaMulticarry, traits::RNG};
use strafe_testing::{
    r::{RString, RTester},
    r_assert_relative_equal,
};

use crate::funcs::opt::{brent::Brent, nmmin::NelderMead, Optim, BFGS, CG, L_BFGS_B, SANN};

pub fn brent_test_inner(
    p: f64,
    lower: f64,
    upper: f64,
    f_func: &dyn Fn(f64) -> f64,
    f_string: String,
) {
    let mut brent = Brent::default();
    brent.lower_bound = lower.into();
    brent.upper_bound = upper.into();
    let rust_ret = brent.optim(p, &f_func).expect("Error in Rust");

    let r_funcs = RString::from_str(&f_string).unwrap();

    let r_ret_par = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};lower={};upper={};{}",
            RString::from_f64(p),
            RString::from_f64(lower),
            RString::from_f64(upper),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, lower=lower, upper=upper, method='Brent')$par")
                .unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    r_assert_relative_equal!(r_ret_par, rust_ret.parameters, 1e-4);

    let r_ret_value = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};lower={};upper={};{}",
            RString::from_f64(p),
            RString::from_f64(lower),
            RString::from_f64(upper),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, lower=lower, upper=upper, method='Brent')$value")
                .unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    r_assert_relative_equal!(r_ret_value, rust_ret.value, 1e-4);

    let r_ret_counts = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};lower={};upper={};{}",
            RString::from_f64(p),
            RString::from_f64(lower),
            RString::from_f64(upper),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, lower=lower, upper=upper, method='Brent')$counts")
                .unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    assert_eq!(r_ret_counts[0] as usize, rust_ret.function_count);
    assert_eq!(r_ret_counts[1] as usize, rust_ret.gradient_count);

    let r_ret_conv = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};lower={};upper={};{}",
            RString::from_f64(p),
            RString::from_f64(lower),
            RString::from_f64(upper),
            r_funcs,
        )))
        .set_display(
            &RString::from_str(
                "optim(p, f_func, lower=lower, upper=upper, method='Brent')$convergence",
            )
            .unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    assert_eq!(r_ret_conv == 0.0, rust_ret.convergence);
}

pub fn cgmin_test_inner(
    p: &[f64],
    f_func: &dyn Fn(Vec<f64>) -> f64,
    f_string: String,
    g_func: Option<Box<dyn Fn(Vec<f64>) -> Vec<f64>>>,
    g_string: Option<String>,
) {
    let mut cg_min = CG::default();
    cg_min.grad_func = g_func;
    let rust_ret = cg_min.optim(p.to_vec(), f_func).expect("Error in Rust");

    let r_funcs =
        RString::from_str(&(f_string + g_string.unwrap_or("g_func=NULL;".to_string()).as_str()))
            .unwrap();

    let r_ret_par = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='CG')$par").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    for (r_ret, &rust_ret) in r_ret_par.into_iter().zip(rust_ret.parameters.iter()) {
        r_assert_relative_equal!(r_ret, rust_ret, 1e-4);
    }

    let r_ret_value = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='CG')$value").unwrap())
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    r_assert_relative_equal!(r_ret_value, rust_ret.value, 1e-4);

    let r_ret_counts = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='CG')$counts").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    assert_eq!(r_ret_counts[0] as usize, rust_ret.function_count);
    assert_eq!(r_ret_counts[1] as usize, rust_ret.gradient_count);

    let r_ret_conv = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, g_func, method='CG')$convergence").unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    assert_eq!(r_ret_conv == 0.0, rust_ret.convergence);
}

pub fn lbfgsb_test_inner(
    p: &[f64],
    f_func: &dyn Fn(Vec<f64>) -> f64,
    f_string: String,
    g_func: Option<Box<dyn Fn(Vec<f64>) -> Vec<f64>>>,
    g_string: Option<String>,
) {
    let mut lbfgsb = L_BFGS_B::default();
    lbfgsb.grad_func = g_func;
    let rust_ret = lbfgsb.optim(p.to_vec(), &f_func).expect("Error in Rust");

    let r_funcs =
        RString::from_str(&(f_string + g_string.unwrap_or("g_func=NULL;".to_string()).as_str()))
            .unwrap();

    let r_ret_par = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='L-BFGS-B')$par").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    for (r_ret, &rust_ret) in r_ret_par.into_iter().zip(rust_ret.parameters.iter()) {
        r_assert_relative_equal!(r_ret, rust_ret, 1e-4);
    }

    let r_ret_value = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, g_func, method='L-BFGS-B')$value").unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    // R can give somewhat different final values, but that's ok as long as it's not too different
    r_assert_relative_equal!(r_ret_value, rust_ret.value, 1e1);

    // let r_ret_counts = RTester::new()
    //     .set_script(&RString::from_string(format!(
    //         "p={};{}",
    //         RString::from_f64_slice(&p),
    //         r_funcs,
    //     )))
    //     .set_display(
    //         &RString::from_str("optim(p, f_func, g_func, method='L-BFGS-B')$counts").unwrap(),
    //     )
    //     .run()
    //     .expect("R running error")
    //     .as_f64_vec()
    //     .expect("R running error");

    // R doesn't keep track of the function and gradient runs in the same way as Rust does
    // assert_eq!(r_ret_counts[0] as usize, rust_ret.function_count);
    // assert_eq!(r_ret_counts[1] as usize, rust_ret.gradient_count);

    let r_ret_conv = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, g_func, method='L-BFGS-B')$convergence").unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    assert_eq!(r_ret_conv == 0.0, rust_ret.convergence);
}

pub fn nmmin_test_inner(p: &[f64], f_func: &dyn Fn(Vec<f64>) -> f64, f_string: String) {
    let mut nelder_mead = NelderMead::default();
    let rust_ret = nelder_mead
        .optim(p.to_vec(), f_func)
        .expect("Error in Rust");

    let r_funcs = RString::from_str(&f_string).unwrap();

    let r_ret_par = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, method='Nelder-Mead')$par").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    for (r_ret, &rust_ret) in r_ret_par.into_iter().zip(rust_ret.parameters.iter()) {
        r_assert_relative_equal!(r_ret, rust_ret, 1e-4);
    }

    let r_ret_value = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, method='Nelder-Mead')$value").unwrap())
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    r_assert_relative_equal!(r_ret_value, rust_ret.value, 1e-4);

    let r_ret_counts = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, method='Nelder-Mead')$counts").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    assert_eq!(r_ret_counts[0] as usize, rust_ret.function_count);
    assert_eq!(r_ret_counts[1] as usize, rust_ret.gradient_count);

    let r_ret_conv = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, method='Nelder-Mead')$convergence").unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    assert_eq!(r_ret_conv == 0.0, rust_ret.convergence);
}

pub fn samin_test_inner(
    p: &[f64],
    seed: u16,
    f_func: &dyn Fn(Vec<f64>) -> f64,
    f_string: String,
    g_func: Option<Box<dyn Fn(Vec<f64>) -> Vec<f64>>>,
    g_string: Option<String>,
) {
    let mut sann = SANN::<MarsagliaMulticarry>::default();
    sann.rng.set_seed(seed as u32);
    sann.grad_func = g_func;
    let rust_ret = sann.optim(p.to_vec(), f_func).expect("Error in Rust");

    let r_funcs =
        RString::from_str(&(f_string + g_string.unwrap_or("g_func=NULL;".to_string()).as_str()))
            .unwrap();

    let r_ret_par = RTester::new()
        .set_rand_type("Marsaglia-Multicarry")
        .set_seed(seed)
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='SANN')$par").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    for (r_ret, &rust_ret) in r_ret_par.into_iter().zip(rust_ret.parameters.iter()) {
        r_assert_relative_equal!(r_ret, rust_ret, 1e-4);
    }

    let r_ret_value = RTester::new()
        .set_rand_type("Marsaglia-Multicarry")
        .set_seed(seed)
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='SANN')$value").unwrap())
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    r_assert_relative_equal!(r_ret_value, rust_ret.value, 1e-2);

    let r_ret_counts = RTester::new()
        .set_rand_type("Marsaglia-Multicarry")
        .set_seed(seed)
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='SANN')$counts").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    assert_eq!(r_ret_counts[0] as usize, rust_ret.function_count);
    assert_eq!(r_ret_counts[1] as usize, rust_ret.gradient_count);

    let r_ret_conv = RTester::new()
        .set_rand_type("Marsaglia-Multicarry")
        .set_seed(seed)
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, g_func, method='SANN')$convergence").unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    assert_eq!(r_ret_conv == 0.0, rust_ret.convergence);
}

pub fn vmmin_test_inner(
    p: &[f64],
    f_func: &dyn Fn(Vec<f64>) -> f64,
    f_string: String,
    g_func: Option<Box<dyn Fn(Vec<f64>) -> Vec<f64>>>,
    g_string: Option<String>,
) {
    let mut bfgs = BFGS::default();
    bfgs.grad_func = g_func;
    let rust_ret = bfgs.optim(p.to_vec(), f_func).expect("Error in Rust");

    let r_funcs =
        RString::from_str(&(f_string + g_string.unwrap_or("g_func=NULL;".to_string()).as_str()))
            .unwrap();

    let r_ret_par = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='BFGS')$par").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    for (r_ret, &rust_ret) in r_ret_par.into_iter().zip(rust_ret.parameters.iter()) {
        r_assert_relative_equal!(r_ret, rust_ret, 1e-4);
    }

    let r_ret_value = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='BFGS')$value").unwrap())
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    r_assert_relative_equal!(r_ret_value, rust_ret.value, 1e-4);

    let r_ret_counts = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(&RString::from_str("optim(p, f_func, g_func, method='BFGS')$counts").unwrap())
        .run()
        .expect("R running error")
        .as_f64_vec()
        .expect("R running error");

    assert_eq!(r_ret_counts[0] as usize, rust_ret.function_count);
    assert_eq!(r_ret_counts[1] as usize, rust_ret.gradient_count);

    let r_ret_conv = RTester::new()
        .set_script(&RString::from_string(format!(
            "p={};{}",
            RString::from_f64_slice(&p),
            r_funcs,
        )))
        .set_display(
            &RString::from_str("optim(p, f_func, g_func, method='BFGS')$convergence").unwrap(),
        )
        .run()
        .expect("R running error")
        .as_f64()
        .expect("R running error");

    assert_eq!(r_ret_conv == 0.0, rust_ret.convergence);
}