minuit2 0.4.2

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

//! Migrad seed generator.
//!
//! Replaces MnSeedGenerator.cxx (Migrad path). Creates the initial MinimumSeed
//! by evaluating the FCN, computing a numerical gradient (not just heuristic),
//! and building Vā‚€ = diag(1/g2_i).

use nalgebra::{DMatrix, DVector};

use crate::fcn::FCNGradient;
use crate::gradient::{
    AnalyticalGradientCalculator, InitialGradientCalculator, Numerical2PGradientCalculator,
};
use crate::minimum::error::MinimumError;
use crate::minimum::parameters::MinimumParameters;
use crate::minimum::seed::MinimumSeed;
use crate::minimum::state::MinimumState;
use crate::mn_fcn::MnFcn;
use crate::strategy::MnStrategy;
use crate::user_transformation::MnUserTransformation;

pub struct MigradSeedGenerator;

impl MigradSeedGenerator {
    /// Generate seed using numerical gradients (central differences).
    pub fn generate(
        fcn: &MnFcn,
        trafo: &MnUserTransformation,
        strategy: &MnStrategy,
    ) -> MinimumSeed {
        let n = trafo.variable_parameters();
        let eps2 = trafo.precision().eps2();

        // 1. Get initial internal parameter values
        let int_values = trafo.initial_internal_values();
        let int_vec = DVector::from_vec(int_values.clone());

        // 2. Evaluate FCN at starting point
        let fval = fcn.call(&int_values);
        let params = MinimumParameters::new(int_vec, fval);

        // 3. Compute heuristic gradient (no FCN calls ā€” just from step sizes)
        let heuristic_calc = InitialGradientCalculator::new(*strategy);
        let heuristic_grad = heuristic_calc.compute(fcn, &params, trafo);

        // 4. Compute numerical gradient (2-point central differences)
        let numerical_calc = Numerical2PGradientCalculator::new(*strategy);
        let gradient = numerical_calc.compute(fcn, &params, trafo, &heuristic_grad);

        // 5. Build Vā‚€ = diag(1/g2_i), fallback to 1.0 for non-positive g2
        let mut v0 = DMatrix::zeros(n, n);
        for i in 0..n {
            let g2i = gradient.g2()[i];
            v0[(i, i)] = if g2i > eps2 { 1.0 / g2i } else { 1.0 };
        }

        let dcovar = 1.0; // approximate: initial V is rough
        let error = MinimumError::new(v0, dcovar);

        // 6. EDM = 0.5 * g^T * V * g
        let edm = {
            let g = gradient.grad();
            let e = error.matrix();
            0.5 * g.dot(&(e * g))
        };

        let state = MinimumState::new(params, error, gradient, edm, fcn.num_of_calls());

        MinimumSeed::new(state, trafo.clone())
    }

    /// Generate seed using analytical gradients from user.
    pub fn generate_with_gradient(
        fcn: &dyn FCNGradient,
        trafo: &MnUserTransformation,
        _strategy: &MnStrategy,
    ) -> MinimumSeed {
        let n = trafo.variable_parameters();
        let eps2 = trafo.precision().eps2();

        // 1. Get initial internal parameter values
        let int_values = trafo.initial_internal_values();
        let int_vec = DVector::from_vec(int_values.clone());

        // 2. Evaluate FCN at starting point (needed for EDM calculation)
        let fval = fcn.value(&trafo.transform(&int_values));
        let params = MinimumParameters::new(int_vec, fval);

        // 3. Compute analytical gradient (user-provided, with g2/gstep heuristics)
        let gradient = AnalyticalGradientCalculator::compute(fcn, trafo, &params);

        // 4. Build Vā‚€ = diag(1/g2_i), fallback to 1.0 for non-positive g2
        let mut v0 = DMatrix::zeros(n, n);
        for i in 0..n {
            let g2i = gradient.g2()[i];
            v0[(i, i)] = if g2i > eps2 { 1.0 / g2i } else { 1.0 };
        }

        let dcovar = 1.0; // approximate: initial V is rough
        let error = MinimumError::new(v0, dcovar);

        // 5. EDM = 0.5 * g^T * V * g
        let edm = {
            let g = gradient.grad();
            let e = error.matrix();
            0.5 * g.dot(&(e * g))
        };

        // Note: no MnFcn call counter here since analytical gradient doesn't eval FCN
        // We'll use state with nfcn=1 (for the initial FCN eval only)
        let state = MinimumState::new(params, error, gradient, edm, 1);

        MinimumSeed::new(state, trafo.clone())
    }

    pub fn call_with_analytical_gradient_calculator(
        fcn: &dyn FCNGradient,
        trafo: &MnUserTransformation,
        strategy: &MnStrategy,
    ) -> MinimumSeed {
        Self::generate_with_gradient(fcn, trafo, strategy)
    }
}