gosh-fire 0.1.1

FIRE algorithm for geometry optimization
Documentation
// [[file:../fire.note::*docs][docs:2]]
/// Represents an optimization problem with cache for avoiding unnecessary
/// function re-evaluations.
///
/// # Examples
///
/// ```ignore
/// let mut x = vec![0.0; 5];
/// let mut prb = CachedProblem::new(&x, f);
/// let d = [0.2; 5];
/// prb.take_line_step(&d, 1.0);
/// let fx = prb.value();
/// let gx = prb.gradient();
/// ```
// docs:2 ends here

// [[file:../fire.note::*imports][imports:1]]
use crate::common::*;
use vecfx::*;
// imports:1 ends here

// [[file:../fire.note::*cached][cached:1]]
#[derive(Clone, Debug)]
pub struct CachedProblem<E>
where
    E: FnMut(&[f64], &mut [f64]) -> f64,
{
    f: E,
    x: Vec<f64>,
    fx: Option<f64>,
    gx: Option<Vec<f64>>,

    epsilon: f64,
    neval: usize,

    // cache previous point
    x_prev: Vec<f64>,
    fx_prev: Option<f64>,
    gx_prev: Option<Vec<f64>>,
}

impl<E> CachedProblem<E>
where
    E: FnMut(&[f64], &mut [f64]) -> f64,
{
    /// Construct a CachedProblem
    ///
    /// # Parameters
    ///
    /// * x: initial position
    /// * f: a closure for function evaluation of value and gradient.
    pub fn new(x: &[f64], f: E) -> Self {
        CachedProblem {
            f,
            epsilon: 1e-8,
            x: x.to_vec(),
            fx: None,
            gx: None,
            neval: 0,
            x_prev: x.to_vec(),
            fx_prev: None,
            gx_prev: None,
        }
    }

    /// The number of function calls
    pub fn ncalls(&self) -> usize {
        self.neval
    }

    /// evaluate function value and gradient at current position
    fn eval(&mut self) -> (f64, &[f64]) {
        let n = self.x.len();

        let gx: &mut [f64] = self.gx.get_or_insert(vec![0.0; n]);
        let v = (self.f)(&self.x, gx);
        self.fx = Some(v);
        self.neval += 1;

        // update previous point
        self.fx_prev = self.fx;
        self.x_prev = self.x.clone();
        self.gx_prev = Some(gx.to_vec());

        (v, gx)
    }

    /// Return function value at current position.
    ///
    /// The function will be evaluated when necessary.
    pub fn value(&mut self) -> f64 {
        match self.fx {
            // found cached value.
            Some(v) => v,
            // first time calculation
            None => {
                let (fx, _) = self.eval();
                fx
            }
        }
    }

    /// Return function value at previous point
    pub fn value_prev(&self) -> f64 {
        match self.fx_prev {
            Some(fx) => fx,
            None => panic!("not evaluated yet"),
        }
    }

    pub fn gradient_prev(&self) -> &[f64] {
        match self.gx_prev {
            Some(ref gx) => gx,
            None => panic!("not evaluated yet"),
        }
    }

    /// Return a reference to gradient at current position.
    ///
    /// The function will be evaluated when necessary.
    pub fn gradient(&mut self) -> &[f64] {
        match self.gx {
            // found cached value.
            Some(ref gx) => gx,
            // first time calculation
            None => {
                let (_, gx) = self.eval();
                gx
            }
        }
    }

    /// Return a reference to current position vector.
    pub fn position(&self) -> &[f64] {
        &self.x
    }

    /// Update position `x` at a prescribed displacement and step size.
    ///
    /// x += step * displ
    pub fn take_line_step(&mut self, displ: &[f64], step: f64) {
        // position changed
        if step * displ.vec2norm() > self.epsilon {
            // update position vector with displacement
            self.x.vecadd(displ, step);
            self.fx = None;
            self.gx = None;
        }
    }

    /// Revert to previous point
    pub fn revert(&mut self) {
        self.fx = self.fx_prev;
        self.x.veccpy(&self.x_prev);
        self.gx = self.gx_prev.clone();
    }
}
// cached:1 ends here