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::{error::Error, ops::DerefMut};

use num_traits::Float;
use r2rs_nmath::{rng::NMathRNG, traits::RNG};

use crate::funcs::opt::{fminfn::fminfn, Optim, OptimOptions, OptimResponse};

pub struct SANN<R: RNG> {
    pub tmax: usize,
    pub ti: f64,
    pub rng: Box<R>,
    pub grad_func: Option<Box<dyn Fn(Vec<f64>) -> Vec<f64>>>,
    pub optim_options: OptimOptions,
}

impl<R: RNG + Default> Default for SANN<R> {
    fn default() -> Self {
        Self {
            tmax: 10,
            ti: 10.0,
            rng: Box::<R>::default(),
            grad_func: None,
            optim_options: OptimOptions {
                parscale: None,
                fnscale: 1.0,
                abstol: f64::neg_infinity(),
                reltol: f64::epsilon().sqrt(),
                maxit: 10_000,
            },
        }
    }
}

impl<R: RNG> Optim for SANN<R> {
    type ParameterType = Vec<f64>;
    fn optim(
        &mut self,
        par: Self::ParameterType,
        min_func: &dyn Fn(Self::ParameterType) -> f64,
    ) -> Result<OptimResponse<Self::ParameterType>, Box<dyn Error>> {
        let npar = par.len();
        let parscale = if let Some(parscale) = &self.optim_options.parscale {
            parscale.clone()
        } else {
            vec![1.0; npar]
        };
        let mut bvec = par
            .iter()
            .zip(parscale.iter())
            .map(|(par, scale)| par / scale)
            .collect::<Vec<_>>();

        let mut grcount = 0;
        let mut fncount = 0;
        let mut min = f64::nan();
        let mut par_out = Vec::new();

        min = samin(
            npar,
            &mut bvec,
            min_func,
            self.grad_func.as_deref(),
            self.optim_options.maxit,
            self.tmax,
            &parscale,
            self.optim_options.fnscale,
            self.ti,
            self.rng.deref_mut(),
        );
        fncount = if npar > 0 {
            self.optim_options.maxit
        } else {
            1
        };
        par_out = bvec.clone();

        for i in 0..npar {
            par_out[i] *= parscale[i];
        }

        Ok(OptimResponse {
            parameters: par_out,
            value: min,
            function_count: fncount,
            gradient_count: grcount,
            convergence: true,
        })
    }
}

const BIG: f64 = 1.0e35;
const E1: f64 = 1.7182818; /* exp(1.0)-1.0 */

/* Given a starting point pb[0..n-1], simulated annealing minimization
   is performed on the function fminfn. The starting temperature
   is input as ti. To make sann work silently set trace to zero.
   sann makes in total maxit function evaluations, tmax
   evaluations at each temperature. Returned quantities are pb
   (the location of the minimum), and yb (the minimum value of
   the function func).  Author: Adrian Trapletti
*/
pub fn samin(
    n: usize,
    pb: &mut [f64],
    min_func: &dyn Fn(Vec<f64>) -> f64,
    gen_func_opt: Option<&dyn Fn(Vec<f64>) -> Vec<f64>>,
    maxit: usize,
    tmax: usize,
    parscale: &[f64],
    fnscale: f64,
    ti: f64,
    rng: &mut impl RNG,
) -> f64 {
    if n == 0 {
        // don't even attempt to optimize
        fminfn(n, pb, parscale, fnscale, min_func)
    } else {
        let mut p = vec![0.0; n];
        let mut ptry = vec![0.0; n];
        let mut yb: f64 = fminfn(n, pb, parscale, fnscale, min_func);
        if !yb.is_finite() {
            yb = BIG;
        }
        p[..n].copy_from_slice(&pb[..n]);
        let mut y = yb;
        let scale = 1.0 / ti;
        let mut its = 1;
        while its < maxit {
            // cool down system
            let t = ti / (its as f64 + E1).ln(); // temperature annealing schedule
            let mut k = 1;
            while (k <= tmax) && (its < maxit) {
                // iterate at constant temperature
                genptry(n, &p, &mut ptry, scale * t, parscale, rng, gen_func_opt);
                let mut ytry: f64 = fminfn(n, &ptry, parscale, fnscale, min_func);
                if !ytry.is_finite() {
                    ytry = BIG;
                }
                let dy = ytry - y;
                if dy <= 0.0 || (rng.unif_rand() < (-dy / t).exp()) {
                    // accept new point?
                    p[..n].copy_from_slice(&ptry[..n]);
                    y = ytry; // update system state p, y
                    if y <= yb {
                        // if system state is best, then update best system state pb, *yb
                        pb[..n].copy_from_slice(&p[..n]);
                        yb = y;
                    }
                }
                its += 1;
                k += 1;
            }
        }
        yb
    }
}

fn genptry(
    n: usize,
    p: &[f64],
    ptry: &mut [f64],
    scale: f64,
    parscale: &[f64],
    rng: &mut impl RNG,
    gen_func_opt: Option<&dyn Fn(Vec<f64>) -> Vec<f64>>,
) {
    if let Some(gen_func) = gen_func_opt {
        // user defined generation of candidate point
        let mut x = vec![0.0; n];
        for i in 0..n {
            x[i] = p[i] * parscale[i];
        }
        let s = gen_func(x);
        for i in 0..n {
            ptry[i] = s[i] / parscale[i];
        }
    } else {
        // default Gaussian Markov kernel
        for i in 0..n {
            ptry[i] = p[i] + scale * rng.norm_rand(); // new candidate point
        }
    }
}