swoop 0.1.0

Simple, lightweight optimisation algorithms in pure Rust
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

//! Algorithms for scalar univariate function optimisation

mod bounded;
pub use bounded::bounded;

mod brent;
pub use brent::brent;

mod golden;
use crate::SwoopErrors;
pub use golden::golden;

/// Struct to represent the result of a scalar univariate function optimisation
#[derive(Debug, Clone)]
pub struct ScalarOptimisationResult {
    /// Value of the objective function
    pub fun: f64,
    /// Number of evaluations of the objective function
    pub nfev: usize,
    /// Whether the optimisation was successful or not
    pub success: bool,
    /// The solution of the optimization
    pub x: f64,
}

/// Trait to implement for a scalar univarite objective function
pub trait ScalarObjectiveFunction {
    /// Method to implement the objective function that will be used for evaluation when
    /// optimising
    fn evaluate(&self, x: f64) -> f64;
}

pub(crate) struct BracketResult {
    pub(crate) bracket: (f64, f64, f64),
    pub(crate) fun_calls: usize,
}

/// Bracket the minimum of the function.
///
/// Given a function and distinct initial points, search in the
/// downhill direction (as defined by the initial points) and return
/// new points `xa, xb, xc` that bracket the minimum of the function
/// `f(xa) > f(xb) < f(xc)`.
///
/// It doesn't always mean that obtained solution will satisfy `xa<=x<=xb`.
///
/// # Parameters
/// * `objective_function` - Objective function to minimize
/// * `xa` - bottom of the bracketing interval
/// * `xb` - top of the bracketing interval
/// * `grow_limit` - Maximum grow limit
/// * `maxiter` - Maximum number of iterations to perform.
pub(crate) fn bracket<T: ScalarObjectiveFunction>(
    objective_function: &T,
    mut xa: f64,
    mut xb: f64,
    grow_limit: f64,
    maxiter: usize,
) -> Result<BracketResult, SwoopErrors> {
    // the golden ratio
    let gold = 1.618_034;

    let very_small_number = 1e-21;
    let mut fa = objective_function.evaluate(xa);
    let mut fb = objective_function.evaluate(xb);

    if fa < fb {
        (xa, xb) = (xb, xa);
        (fa, fb) = (fb, fa);
    }

    let mut xc = xb + gold * (xb - xa);
    let mut fc = objective_function.evaluate(xc);
    let mut fun_calls = 3usize;
    let mut iter = 0usize;

    let mut tmp1: f64;
    let mut tmp2: f64;
    let mut val: f64;

    let mut denom: f64;
    let mut w: f64;
    let mut wlim: f64;

    let mut fw: f64;

    while fc < fb {
        tmp1 = (xb - xa) * (fb - fc);
        tmp2 = (xb - xc) * (fb - fa);
        val = tmp2 - tmp1;

        if val.abs() < very_small_number {
            denom = 2.0 * very_small_number;
        } else {
            denom = 2.0 * val;
        }

        w = xb - ((xb - xc) * tmp2 - (xb - xa) * tmp1) / denom;
        wlim = xb + grow_limit * (xc - xb);

        if iter > maxiter {
            return Err(SwoopErrors::MaxIterExceeded);
        }
        iter += 1usize;

        if (w - xc) * (xb - w) > 0.0 {
            fw = objective_function.evaluate(w);
            fun_calls += 1usize;
            if fw < fc {
                xa = xb;
                xb = w;
                //fa = fb;
                //fb = fw;
                return Ok(BracketResult {
                    bracket: (xa, xb, xc),
                    fun_calls,
                });
            } else if fw > fb {
                xc = w;
                //fc = fw;
                return Ok(BracketResult {
                    bracket: (xa, xb, xc),
                    fun_calls,
                });
            }

            w = xc + gold * (xc - xb);
            fw = objective_function.evaluate(w);
            fun_calls += 1usize;
        } else if (w - wlim) * (wlim - xc) >= 0.0 {
            w = wlim;
            fw = objective_function.evaluate(w);
            fun_calls += 1usize;
        } else if (w - wlim) * (xc - w) > 0.0 {
            fw = objective_function.evaluate(w);
            fun_calls += 1usize;
            if fw < fc {
                xb = xc;
                xc = w;
                w = xc + gold * (xc - xb);
                fb = fc;
                fc = fw;
                fw = objective_function.evaluate(w);
                fun_calls += 1usize;
            }
        } else {
            w = xc + gold * (xc - xb);
            fw = objective_function.evaluate(w);
            fun_calls += 1usize;
        }
        xa = xb;
        xb = xc;
        xc = w;
        fa = fb;
        fb = fc;
        fc = fw;
    }

    Ok(BracketResult {
        bracket: (xa, xb, xc),
        fun_calls,
    })
}

/// Bracket the minimum of the function with default arguments
///
/// # Parameters
/// * `objective_function` - Objective function to minimize
pub(crate) fn bracket_default<T: ScalarObjectiveFunction>(
    objective_function: &T,
) -> Result<BracketResult, SwoopErrors> {
    bracket(objective_function, 0.0f64, 1.0f64, 110.0f64, 1000usize)
}