metalforge 0.1.1

//! Dynamically Dimensioned Search (Tolson & Shoemaker 2007, WRR 43, W01413).
//!
//! A greedy single-solution global optimizer that scales the *number* of
//! perturbed dimensions down as the budget is consumed: early on it explores
//! many dimensions, late on it fine-tunes one at a time. Parsimonious — a
//! single control parameter `r` — and a workhorse for hydrological calibration,
//! which is why it is absent from generic optimization libraries but central to
//! forge. Migrated from the implementation validated in `rainflow-core` against
//! `airGR::Calibration_Michel`.

use super::{clamp, sample, Evaluator, Optimizer};
use crate::problem::Problem;
use crate::rng::Rng;
use crate::solution::{Report, Solution};
use crate::termination::Termination;

/// DDS configuration.
#[derive(Debug, Clone, Copy)]
pub struct Dds {
    /// Neighborhood perturbation size, as a fraction of each variable's range.
    /// The paper's robust default is `0.2`.
    pub r: f64,
    /// RNG seed; same seed + same problem + same budget ⇒ same result.
    pub seed: u64,
}

impl Default for Dds {
    fn default() -> Self {
        Dds { r: 0.2, seed: 42 }
    }
}

impl Optimizer for Dds {
    /// Minimizes `problem` within its bounds using DDS, starting from a uniform
    /// random point.
    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
        self.optimize_from(problem, term, None)
    }

    fn with_seed(&self, seed: u64) -> Self {
        Dds { seed, ..*self }
    }
}

impl Dds {
    /// Like [`optimize`](Optimizer::optimize) but with an explicit starting
    /// point. If `init` has the wrong length it is ignored and a random start
    /// is used.
    pub fn optimize_from(
        &self,
        problem: &dyn Problem,
        term: &Termination,
        init: Option<&[f64]>,
    ) -> Report {
        let bounds = problem.bounds();
        let dim = bounds.len();
        let mut rng = Rng::new(self.seed);

        let start: Vec<f64> = match init {
            Some(x0) if x0.len() == dim => x0.to_vec(),
            _ => sample(bounds, &mut rng),
        };
        let value = problem.objective(&start);
        let mut ev = Evaluator::new(
            problem,
            term,
            Solution {
                x: start,
                value: if value.is_finite() {
                    value
                } else {
                    f64::INFINITY
                },
            },
        );

        let m = (term.max_evaluations.max(2)) as f64;
        let mut i = 1usize;
        let mut candidate = vec![0.0; dim];
        while !ev.done() {
            // P(perturb a dimension) decays from ~1 toward 1/m over the budget.
            let p = 1.0 - (i as f64).ln() / m.ln();

            candidate.copy_from_slice(&ev.best.x);
            let mut perturbed = 0;
            for (j, &(lo, hi)) in bounds.iter().enumerate() {
                if rng.uniform() < p {
                    candidate[j] = perturb(ev.best.x[j], lo, hi, self.r, &mut rng);
                    perturbed += 1;
                }
            }
            if perturbed == 0 {
                // Always perturb at least one randomly chosen dimension.
                let j = rng.index(dim);
                let (lo, hi) = bounds[j];
                candidate[j] = perturb(ev.best.x[j], lo, hi, self.r, &mut rng);
            }

            // `eval` accepts greedily (best updates iff strictly better), which
            // is exactly DDS's acceptance rule.
            ev.eval(&candidate);
            i += 1;
        }

        ev.finish()
    }
}

/// One-dimensional DDS neighborhood move with boundary reflection
/// (Tolson & Shoemaker 2007, eq. 4): reflect once at the violated bound, and if
/// still outside, clamp to that bound.
fn perturb(x: f64, lo: f64, hi: f64, r: f64, rng: &mut Rng) -> f64 {
    let range = hi - lo;
    let mut xn = x + range * r * rng.normal();
    if xn < lo {
        xn = lo + (lo - xn);
        if xn > hi {
            xn = lo;
        }
    } else if xn > hi {
        xn = hi - (xn - hi);
        if xn < lo {
            xn = hi;
        }
    }
    clamp(xn, lo, hi)
}