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

//! Simplex vertex storage with best/worst tracking.
//!
//! Replaces SimplexParameters.h/.cxx. Stores N+1 vertices (each a parameter
//! vector + function value) and tracks indices of the best (lowest) and
//! worst (highest) function values.

#[derive(Debug, Clone)]
pub struct SimplexParameters {
    /// (function_value, parameter_vector) pairs for each vertex.
    params: Vec<(f64, Vec<f64>)>,
    /// Index of vertex with lowest function value.
    jlow: usize,
    /// Index of vertex with highest function value.
    jhigh: usize,
}

impl SimplexParameters {
    pub fn new(params: Vec<(f64, Vec<f64>)>) -> Self {
        let (jlow, jhigh) = Self::find_extremes(&params);
        Self {
            params,
            jlow,
            jhigh,
        }
    }

    fn find_extremes(params: &[(f64, Vec<f64>)]) -> (usize, usize) {
        let mut jlow = 0;
        let mut jhigh = 0;
        for (i, (fval, _)) in params.iter().enumerate() {
            if *fval < params[jlow].0 {
                jlow = i;
            }
            if *fval > params[jhigh].0 {
                jhigh = i;
            }
        }
        (jlow, jhigh)
    }

    /// Update a vertex and recompute extremes.
    pub fn update(&mut self, index: usize, fval: f64, vec: Vec<f64>) {
        self.params[index] = (fval, vec);
        let (jlow, jhigh) = Self::find_extremes(&self.params);
        self.jlow = jlow;
        self.jhigh = jhigh;
    }

    pub fn params(&self) -> &[(f64, Vec<f64>)] {
        &self.params
    }

    /// Index of best (lowest) vertex.
    pub fn jlow(&self) -> usize {
        self.jlow
    }

    pub fn jl(&self) -> usize {
        self.jlow()
    }

    /// Index of worst (highest) vertex.
    pub fn jhigh(&self) -> usize {
        self.jhigh
    }

    pub fn jh(&self) -> usize {
        self.jhigh()
    }

    /// Function value at best vertex.
    pub fn fval_best(&self) -> f64 {
        self.params[self.jlow].0
    }

    /// Function value at worst vertex.
    pub fn fval_worst(&self) -> f64 {
        self.params[self.jhigh].0
    }

    /// Parameter vector at best vertex.
    pub fn best(&self) -> &[f64] {
        &self.params[self.jlow].1
    }

    /// Compatibility helper for the current direction-like span in simplex space.
    pub fn dirin(&self) -> Vec<f64> {
        if self.params.is_empty() {
            return Vec::new();
        }
        let best = &self.params[self.jlow].1;
        let worst = &self.params[self.jhigh].1;
        best.iter().zip(worst).map(|(b, w)| w - b).collect()
    }

    /// Number of vertices (N+1).
    pub fn len(&self) -> usize {
        self.params.len()
    }

    /// Whether the simplex has no vertices.
    pub fn is_empty(&self) -> bool {
        self.params.is_empty()
    }

    /// EDM estimate: difference between worst and best function values.
    pub fn edm(&self) -> f64 {
        self.fval_worst() - self.fval_best()
    }

    /// Index of second-worst vertex.
    pub fn jsecond_high(&self) -> usize {
        let mut jsec = if self.jhigh == 0 { 1 } else { 0 };
        for (i, (fval, _)) in self.params.iter().enumerate() {
            if i != self.jhigh && *fval > self.params[jsec].0 {
                jsec = i;
            }
        }
        jsec
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn extremes() {
        let params = vec![
            (3.0, vec![1.0, 0.0]),
            (1.0, vec![0.0, 1.0]),
            (5.0, vec![1.0, 1.0]),
        ];
        let sp = SimplexParameters::new(params);
        assert_eq!(sp.jlow(), 1);
        assert_eq!(sp.jhigh(), 2);
        assert!((sp.fval_best() - 1.0).abs() < 1e-15);
        assert!((sp.fval_worst() - 5.0).abs() < 1e-15);
    }

    #[test]
    fn update_recomputes() {
        let params = vec![(3.0, vec![1.0]), (1.0, vec![2.0]), (5.0, vec![3.0])];
        let mut sp = SimplexParameters::new(params);
        sp.update(2, 0.5, vec![4.0]); // vertex 2 is now best
        assert_eq!(sp.jlow(), 2);
        assert!((sp.fval_best() - 0.5).abs() < 1e-15);
    }
}