forge_core/algo/
dds.rs

1//! Dynamically Dimensioned Search (Tolson & Shoemaker 2007, WRR 43, W01413).
2//!
3//! A greedy single-solution global optimizer that scales the *number* of
4//! perturbed dimensions down as the budget is consumed: early on it explores
5//! many dimensions, late on it fine-tunes one at a time. Parsimonious — a
6//! single control parameter `r` — and a workhorse for hydrological calibration,
7//! which is why it is absent from generic optimization libraries but central to
8//! forge. Migrated from the implementation validated in `rainflow-core` against
9//! `airGR::Calibration_Michel`.
10
11use super::{clamp, sample, Evaluator, Optimizer};
12use crate::problem::Problem;
13use crate::rng::Rng;
14use crate::solution::{Report, Solution};
15use crate::termination::Termination;
16
17/// DDS configuration.
18#[derive(Debug, Clone, Copy)]
19pub struct Dds {
20    /// Neighborhood perturbation size, as a fraction of each variable's range.
21    /// The paper's robust default is `0.2`.
22    pub r: f64,
23    /// RNG seed; same seed + same problem + same budget ⇒ same result.
24    pub seed: u64,
25}
26
27impl Default for Dds {
28    fn default() -> Self {
29        Dds { r: 0.2, seed: 42 }
30    }
31}
32
33impl Optimizer for Dds {
34    /// Minimizes `problem` within its bounds using DDS, starting from a uniform
35    /// random point.
36    fn optimize(&self, problem: &dyn Problem, term: &Termination) -> Report {
37        self.optimize_from(problem, term, None)
38    }
39
40    fn with_seed(&self, seed: u64) -> Self {
41        Dds { seed, ..*self }
42    }
43}
44
45impl Dds {
46    /// Like [`optimize`](Optimizer::optimize) but with an explicit starting
47    /// point.
48    ///
49    /// # Panics
50    /// If the problem is invalid (see [`crate::problem::validate`]) or `init`
51    /// has the wrong length — a shape mismatch is a caller bug that must not
52    /// be silently papered over with a random start.
53    pub fn optimize_from(
54        &self,
55        problem: &dyn Problem,
56        term: &Termination,
57        init: Option<&[f64]>,
58    ) -> Report {
59        crate::problem::validate(problem).unwrap_or_else(|e| panic!("Dds: invalid problem: {e}"));
60        let bounds = problem.bounds();
61        let dim = bounds.len();
62        let mut rng = Rng::new(self.seed);
63
64        let start: Vec<f64> = match init {
65            Some(x0) => {
66                assert_eq!(
67                    x0.len(),
68                    dim,
69                    "Dds::optimize_from: init has length {} but the problem has {} variables",
70                    x0.len(),
71                    dim
72                );
73                x0.to_vec()
74            }
75            None => sample(bounds, &mut rng),
76        };
77        let value = problem.objective(&start);
78        let mut ev = Evaluator::new(
79            problem,
80            term,
81            Solution {
82                x: start,
83                value: if value.is_finite() {
84                    value
85                } else {
86                    f64::INFINITY
87                },
88            },
89        );
90
91        // The schedule runs over the *candidate* evaluations that remain after
92        // the initial point (Tolson & Shoemaker's m excludes initialization),
93        // so the final iteration reaches P ≈ 0 as published.
94        let m = (term.max_evaluations.saturating_sub(1)).max(2) as f64;
95        let mut i = 1usize;
96        let mut candidate = vec![0.0; dim];
97        while !ev.done() {
98            // P(perturb a dimension) decays from ~1 toward 1/m over the budget.
99            let p = 1.0 - (i as f64).ln() / m.ln();
100
101            candidate.copy_from_slice(&ev.best.x);
102            let mut perturbed = 0;
103            for (j, &(lo, hi)) in bounds.iter().enumerate() {
104                if rng.uniform() < p {
105                    candidate[j] = perturb(ev.best.x[j], lo, hi, self.r, &mut rng);
106                    perturbed += 1;
107                }
108            }
109            if perturbed == 0 {
110                // Always perturb at least one randomly chosen dimension.
111                let j = rng.index(dim);
112                let (lo, hi) = bounds[j];
113                candidate[j] = perturb(ev.best.x[j], lo, hi, self.r, &mut rng);
114            }
115
116            // `eval` accepts greedily (best updates on `<=` for finite values,
117            // the paper's rule — ties may drift along plateaus).
118            ev.eval(&candidate);
119            i += 1;
120        }
121
122        ev.finish()
123    }
124}
125
126/// One-dimensional DDS neighborhood move with boundary reflection
127/// (Tolson & Shoemaker 2007, eq. 4): reflect once at the violated bound, and if
128/// still outside, clamp to that bound.
129fn perturb(x: f64, lo: f64, hi: f64, r: f64, rng: &mut Rng) -> f64 {
130    let range = hi - lo;
131    let mut xn = x + range * r * rng.normal();
132    if xn < lo {
133        xn = lo + (lo - xn);
134        if xn > hi {
135            xn = lo;
136        }
137    } else if xn > hi {
138        xn = hi - (xn - hi);
139        if xn < lo {
140            xn = hi;
141        }
142    }
143    clamp(xn, lo, hi)
144}
forge_core/algo/dds.rs

forge_core/algo/
dds.rs
