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)

/*  BFGS variable-metric method, based on Pascal code
in J.C. Nash, `Compact Numerical Methods for Computers', 2nd edition,
converted by p2c then re-crafted by B.D. Ripley */
use std::{convert::TryFrom, error::Error};

use nalgebra::DMatrix;
use num_traits::Float;

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

pub struct BFGS {
    pub nREPORT: i32,
    pub grad_func: Option<Box<dyn Fn(Vec<f64>) -> Vec<f64>>>,
    pub optim_options: OptimOptions,
}

impl Default for BFGS {
    fn default() -> Self {
        Self {
            nREPORT: 10,
            grad_func: None,
            optim_options: OptimOptions {
                parscale: None,
                fnscale: 1.0,
                abstol: f64::neg_infinity(),
                reltol: f64::epsilon().sqrt(),
                maxit: 100,
            },
        }
    }
}

impl Optim for BFGS {
    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 ndeps = vec![1e-3; npar];
        let bounds = None;

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

        let mask = vec![true; par.len()];
        (par_out, min, fncount, grcount, fail) = vmmin(
            npar,
            &bvec,
            min_func,
            self.grad_func.as_deref(),
            self.optim_options.maxit,
            &mask,
            &ndeps,
            &parscale,
            self.optim_options.fnscale,
            bounds,
            self.optim_options.abstol,
            self.optim_options.reltol,
            self.nREPORT,
        )
        .expect("Error calculating vmmin");

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

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

const RELTEST: f64 = 10.0;
const ACCTOL: f64 = 0.0001;
const STEPREDN: f64 = 0.2;

pub fn vmmin(
    n0: usize,
    b: &[f64],
    min_func: &dyn Fn(Vec<f64>) -> f64,
    grad_func: Option<&dyn Fn(Vec<f64>) -> Vec<f64>>,
    maxit: usize,
    mask: &[bool],
    ndeps: &[f64],
    parscale: &[f64],
    fnscale: f64,
    bounds: Option<(Vec<f64>, Vec<f64>)>,
    abstol: f64,
    reltol: f64,
    nREPORT: i32,
) -> Result<(Vec<f64>, f64, usize, usize, bool), Box<dyn Error>> {
    let mut fail = false;
    if maxit == 0 {
        let mut b = b.to_vec();
        let min = fminfn(n0, &b, parscale, fnscale, min_func);
        Ok((b, min, 0, 0, fail))
    } else {
        if nREPORT <= 0 {
            return Err(Box::try_from("nREPORT must be > 0").unwrap());
        }
        let mut l = vec![0; n0];

        let mut n = 0;
        for (i, &mask_i) in mask.iter().take(n0).enumerate() {
            if mask_i {
                l[n] = i;
                n += 1;
            }
        }
        let mut t = vec![0.0; n];
        let mut x = vec![0.0; n];
        let mut c = vec![0.0; n];
        let mut B = DMatrix::<f64>::zeros(n, n);
        let mut f: f64 = fminfn(n0, b, parscale, fnscale, min_func);
        if !f.is_finite() {
            return Err(Box::try_from("initial value in vmmin is not finite").unwrap());
        }
        let mut fmin = f;
        let mut funcount = 1;
        let mut gradcount = 1;
        let (mut b, mut g) = fmingr(
            n0,
            b,
            ndeps,
            parscale,
            fnscale,
            bounds.clone(),
            min_func,
            grad_func,
        );
        let mut iter = 1;
        let mut ilast = gradcount;

        let mut count = 0;
        loop {
            if ilast == gradcount {
                for i in 0..n {
                    for j in 0..i {
                        B[(i, j)] = 0.0;
                    }
                    B[(i, i)] = 1.0;
                }
            }

            for i in 0..n {
                x[i] = b[l[i]];
                c[i] = g[l[i]];
            }

            let mut gradproj = 0.0;

            for i in 0..n {
                let mut s = 0.0;
                for j in 0..=i {
                    s -= B[(i, j)] * g[l[j]];
                }
                for j in (i + 1)..n {
                    s -= B[(j, i)] * g[l[j]];
                }
                t[i] = s;
                gradproj += s * g[l[i]];
            }

            if gradproj < 0.0 {
                // search direction is downhill
                let mut steplength = 1.0;
                let mut accpoint = false;

                loop {
                    count = 0;

                    for i in 0..n {
                        b[l[i]] = x[i] + steplength * t[i];
                        if RELTEST + x[i] == RELTEST + b[l[i]] {
                            // no change
                            count += 1;
                        }
                    }

                    if count < n {
                        f = fminfn(n0, &b, parscale, fnscale, min_func);
                        funcount += 1;
                        accpoint = f.is_finite() && f <= fmin + gradproj * steplength * ACCTOL;
                        if !accpoint {
                            steplength *= STEPREDN;
                        }
                    }

                    if count == n || accpoint {
                        break;
                    }
                }

                let enough = f > abstol && (f - fmin).abs() > reltol * (fmin.abs() + reltol);

                if !enough {
                    // stop if value is small or if relative change is low
                    count = n;
                    fmin = f;
                }

                if count < n {
                    // making progress
                    fmin = f;
                    (b, g) = fmingr(
                        n0,
                        &b,
                        ndeps,
                        parscale,
                        fnscale,
                        bounds.clone(),
                        min_func,
                        grad_func,
                    );
                    gradcount += 1;
                    iter += 1;

                    let mut D1 = 0.0;
                    for i in 0..n {
                        t[i] *= steplength;
                        c[i] = g[l[i]] - c[i];
                        D1 += t[i] * c[i];
                    }

                    if D1 > 0.0 {
                        let mut D2 = 0.0;
                        for i in 0..n {
                            let mut s = 0.0;
                            for j in 0..=i {
                                s += B[(i, j)] * c[j];
                            }
                            for j in (i + 1)..n {
                                s += B[(j, i)] * c[j];
                            }
                            x[i] = s;
                            D2 += s * c[i];
                        }
                        D2 = 1.0 + D2 / D1;
                        for i in 0..n {
                            for j in 0..=i {
                                B[(i, j)] += (D2 * t[i] * t[j] - x[i] * t[j] - t[i] * x[j]) / D1;
                            }
                        }
                    } else {
                        // D1 < 0.0
                        ilast = gradcount;
                    }
                } else {
                    // no progress
                    if ilast < gradcount {
                        count = 0;
                        ilast = gradcount;
                    }
                }
            } else {
                // uphill search
                count = 0;
                if ilast == gradcount {
                    count = n;
                } else {
                    ilast = gradcount;
                }
                // Resets unless has just been reset
            }

            if iter >= maxit {
                fail = true;
                break;
            }

            if gradcount - ilast > 2 * n {
                // periodic restart
                ilast = gradcount;
            }

            if count == n && ilast == gradcount {
                break;
            }
        }
        Ok((b, fmin, funcount, gradcount, fail))
    }
}