minuit2 0.4.1

Pure Rust port of CERN Minuit2 parameter optimization engine
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249

//! Two-point numerical gradient calculator.
//!
//! Replaces Numerical2PGradientCalculator.cxx. Computes the gradient using
//! central differences: g_i = (f(x+h) - f(x-h)) / 2h, with adaptive step
//! sizing refined over multiple cycles. Also computes g2 (second derivative
//! estimate) and gstep (optimal step size).

use nalgebra::DVector;

use crate::minimum::gradient::FunctionGradient;
use crate::minimum::parameters::MinimumParameters;
use crate::mn_fcn::MnFcn;
use crate::strategy::MnStrategy;
use crate::user_transformation::MnUserTransformation;

pub struct Numerical2PGradientCalculator {
    strategy: MnStrategy,
}

impl Numerical2PGradientCalculator {
    pub fn new(strategy: MnStrategy) -> Self {
        Self { strategy }
    }

    /// Compute gradient from scratch (no previous gradient available).
    /// Uses the heuristic gradient's gstep as initial step sizes.
    pub fn compute(
        &self,
        fcn: &MnFcn,
        params: &MinimumParameters,
        trafo: &MnUserTransformation,
        initial_gradient: &FunctionGradient,
    ) -> FunctionGradient {
        let n = trafo.variable_parameters();
        let eps2 = trafo.precision().eps2();
        let fcnmin = params.fval();
        let dfmin = 8.0 * eps2 * (fcnmin.abs() + fcn.up());
        let vrysml = 8.0 * eps2 * eps2;

        let x = params.vec();
        let ncycles = self.strategy.grad_ncycles();
        let step_tol = self.strategy.grad_step_tol();
        let grad_tol = self.strategy.grad_tol();

        let mut grad = DVector::zeros(n);
        let mut g2 = DVector::zeros(n);
        let mut gstep = DVector::zeros(n);

        for i in 0..n {
            let ext_idx = trafo.ext_of_int(i);
            let xi = x[i];
            let p = &trafo.parameters()[ext_idx];
            let has_limits = p.has_limits() || p.has_lower_limit() || p.has_upper_limit();

            // Initial step from heuristic gradient
            let mut gstepi = initial_gradient.gstep()[i].max(vrysml);
            let mut g2i = initial_gradient.g2()[i];

            // Ncycles of refinement
            for _cycle in 0..ncycles {
                // Optimal step: balance truncation vs roundoff error
                let optstp = (dfmin / (g2i.abs() + eps2)).sqrt();
                let mut step = optstp.max(0.1 * gstepi.abs());

                // Bounded parameter: cap step at 0.5
                if has_limits {
                    step = step.min(0.5);
                }

                // Clamp step
                let stpmax = 10.0 * gstepi.abs();
                let stpmin = vrysml.max(8.0 * eps2 * xi.abs());
                step = step.clamp(stpmin, stpmax);

                let stepb4 = gstepi;
                let grdb4 = grad[i];

                gstepi = step;

                // Central differences: f(x+h) - f(x-h)
                let mut xp = x.clone();
                let mut xm = x.clone();
                xp[i] = xi + step;
                xm[i] = xi - step;

                let fp = fcn.call(xp.as_slice());
                let fm = fcn.call(xm.as_slice());

                let grdi = 0.5 * (fp - fm) / step;
                let g2i_new = (fp + fm - 2.0 * fcnmin) / (step * step);

                grad[i] = grdi;
                g2[i] = g2i_new;
                gstep[i] = gstepi;
                g2i = g2i_new;

                // Check convergence: step stabilization
                if _cycle > 0 {
                    let step_change = (gstepi - stepb4).abs() / gstepi.abs();
                    if step_change < step_tol {
                        break;
                    }

                    // Gradient stabilization
                    let grad_change = (grdi - grdb4).abs() / (grdi.abs() + dfmin / step);
                    if grad_change < grad_tol {
                        break;
                    }
                }
            }
        }

        FunctionGradient::new(grad, g2, gstep)
    }

    /// Compute gradient using previous gradient's step sizes as starting point.
    /// More efficient than `compute()` since step sizes are already tuned.
    pub fn compute_with_previous(
        &self,
        fcn: &MnFcn,
        params: &MinimumParameters,
        trafo: &MnUserTransformation,
        previous: &FunctionGradient,
    ) -> FunctionGradient {
        let n = trafo.variable_parameters();
        let eps2 = trafo.precision().eps2();
        let fcnmin = params.fval();
        let dfmin = 8.0 * eps2 * (fcnmin.abs() + fcn.up());
        let vrysml = 8.0 * eps2 * eps2;

        let x = params.vec();
        let ncycles = self.strategy.grad_ncycles();
        let step_tol = self.strategy.grad_step_tol();
        let grad_tol = self.strategy.grad_tol();

        let mut grad = DVector::zeros(n);
        let mut g2 = DVector::zeros(n);
        let mut gstep = DVector::zeros(n);

        for i in 0..n {
            let ext_idx = trafo.ext_of_int(i);
            let xi = x[i];
            let p = &trafo.parameters()[ext_idx];
            let has_limits = p.has_limits() || p.has_lower_limit() || p.has_upper_limit();

            // Start from previous step sizes
            let mut gstepi = previous.gstep()[i].max(vrysml);
            let mut g2i = previous.g2()[i];

            for _cycle in 0..ncycles {
                let optstp = (dfmin / (g2i.abs() + eps2)).sqrt();
                let mut step = optstp.max(0.1 * gstepi.abs());

                if has_limits {
                    step = step.min(0.5);
                }

                let stpmax = 10.0 * gstepi.abs();
                let stpmin = vrysml.max(8.0 * eps2 * xi.abs());
                step = step.clamp(stpmin, stpmax);

                let stepb4 = gstepi;
                let grdb4 = grad[i];

                gstepi = step;

                let mut xp = x.clone();
                let mut xm = x.clone();
                xp[i] = xi + step;
                xm[i] = xi - step;

                let fp = fcn.call(xp.as_slice());
                let fm = fcn.call(xm.as_slice());

                let grdi = 0.5 * (fp - fm) / step;
                let g2i_new = (fp + fm - 2.0 * fcnmin) / (step * step);

                grad[i] = grdi;
                g2[i] = g2i_new;
                gstep[i] = gstepi;
                g2i = g2i_new;

                if _cycle > 0 {
                    let step_change = (gstepi - stepb4).abs() / gstepi.abs();
                    if step_change < step_tol {
                        break;
                    }

                    let grad_change = (grdi - grdb4).abs() / (grdi.abs() + dfmin / step);
                    if grad_change < grad_tol {
                        break;
                    }
                }
            }
        }

        FunctionGradient::new(grad, g2, gstep)
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::fcn::FCN;
    use crate::parameter::MinuitParameter;
    use crate::user_transformation::MnUserTransformation;

    struct Quadratic;
    impl FCN for Quadratic {
        fn value(&self, p: &[f64]) -> f64 {
            // f(x,y) = x² + 4y²
            p[0] * p[0] + 4.0 * p[1] * p[1]
        }
    }

    #[test]
    fn numerical_gradient_quadratic() {
        let params = vec![
            MinuitParameter::new(0, "x", 3.0, 0.1),
            MinuitParameter::new(1, "y", 2.0, 0.1),
        ];
        let trafo = MnUserTransformation::new(params);
        let fcn = MnFcn::new(&Quadratic, &trafo);
        let strategy = MnStrategy::default();

        // Evaluate at (3, 2) → f = 9 + 16 = 25
        let x = DVector::from_vec(vec![3.0, 2.0]);
        let min_params = MinimumParameters::new(x, 25.0);

        // Heuristic initial gradient (just for step sizes)
        let ig = crate::gradient::InitialGradientCalculator::new(strategy);
        let init_grad = ig.compute(&fcn, &min_params, &trafo);

        let calc = Numerical2PGradientCalculator::new(strategy);
        let grad = calc.compute(&fcn, &min_params, &trafo, &init_grad);

        // df/dx = 2x = 6, df/dy = 8y = 16
        assert!(
            (grad.grad()[0] - 6.0).abs() < 0.01,
            "dfdx should be ~6.0, got {}",
            grad.grad()[0]
        );
        assert!(
            (grad.grad()[1] - 16.0).abs() < 0.1,
            "dfdy should be ~16.0, got {}",
            grad.grad()[1]
        );
    }
}