argmin_testfunctions/
eggholder.rs

1// Copyright 2018-2024 argmin developers
2//
3// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
4// http://apache.org/licenses/LICENSE-2.0> or the MIT license <LICENSE-MIT or
5// http://opensource.org/licenses/MIT>, at your option. This file may not be
6// copied, modified, or distributed except according to those terms.
7
8//! # Eggholder test function.
9//!
10//! Defined as
11//!
12//! $$
13//! f(x_1, x_2) = -(x_2 + 47) \cdot \sin\left( \sqrt{\left| x_2 + \frac{x_1}{2} + 47 \right|} \right) -
14//!                x_1 \cdot \sin\left(\sqrt{\left|x_1 - (x_2 + 47)\right|}\right)
15//! $$
16//!
17//! where $x_i \in [-512,\\,512]$.
18//!
19//! The global minimum is at $f(x_1,\\,x_2) = f(512,\\,404.2319) = -959.6407$.
20
21use num::{Float, FromPrimitive};
22
23/// Eggholder test function.
24///
25/// Defined as
26///
27/// $$
28/// f(x_1, x_2) = -(x_2 + 47) \cdot \sin\left( \sqrt{\left| x_2 + \frac{x_1}{2} + 47 \right|} \right) -
29///                x_1 \cdot \sin\left(\sqrt{\left|x_1 - (x_2 + 47)\right|}\right)
30/// $$
31///
32/// where $x_i \in [-512,\\,512]$.
33///
34/// The global minimum is at $f(x_1,\\,x_2) = f(512,\\,404.2319) = -959.6407$.
35pub fn eggholder<T>(param: &[T; 2]) -> T
36where
37    T: Float + FromPrimitive,
38{
39    let [x1, x2] = *param;
40    let n47 = T::from_f64(47.0).unwrap();
41    -(x2 + n47)
42        * (x2 + x1 / T::from_f64(2.0).unwrap() + n47)
43            .abs()
44            .sqrt()
45            .sin()
46        - x1 * (x1 - (x2 + n47)).abs().sqrt().sin()
47}
48
49/// Derivative of Eggholder test function.
50pub fn eggholder_derivative<T>(param: &[T; 2]) -> [T; 2]
51where
52    T: Float + FromPrimitive,
53{
54    let [x1, x2] = *param;
55
56    let eps = T::epsilon();
57    let n0 = T::from_f64(0.0).unwrap();
58    let n2 = T::from_f64(2.0).unwrap();
59    let n4 = T::from_f64(4.0).unwrap();
60    let n47 = T::from_f64(47.0).unwrap();
61
62    let x1mx2m47 = x1 - x2 - n47;
63    let x1mx2m47abs = x1mx2m47.abs();
64    let x1mx2m47abssqrt = x1mx2m47abs.sqrt();
65    let x1mx2m47abssqrtsin = x1mx2m47abssqrt.sin();
66    let x1mx2m47abssqrtcos = x1mx2m47abssqrt.cos();
67    let x1hpx2p47 = x1 / n2 + x2 + n47;
68    let x1hpx2p47abs = x1hpx2p47.abs();
69    let x1hpx2p47abssqrt = x1hpx2p47abs.sqrt();
70    let x1hpx2p47abssqrtsin = x1hpx2p47abssqrt.sin();
71    let x1hpx2p47abssqrtcos = x1hpx2p47abssqrt.cos();
72    let x2mx1p47 = x2 - x1 + n47;
73    let x2mx1p47abs = x2mx1p47.abs();
74    let x2mx1p47abssqrt = x2mx1p47abs.sqrt();
75    let x2mx1p47abssqrtcos = x2mx1p47abssqrt.cos();
76
77    [
78        -x1mx2m47abssqrtsin
79            - if x1mx2m47abs <= eps {
80                n0
81            } else {
82                (x1 * x1mx2m47 * x1mx2m47abssqrtcos) / (n2 * x1mx2m47abssqrt.powi(3))
83            }
84            - if x1hpx2p47abs <= eps {
85                n0
86            } else {
87                ((x2 + n47) * x1hpx2p47abssqrtcos * x1hpx2p47) / (n4 * x1hpx2p47abssqrt.powi(3))
88            },
89        -x1hpx2p47abssqrtsin
90            - if x1hpx2p47abs <= eps {
91                n0
92            } else {
93                ((x2 + n47) * x1hpx2p47 * x1hpx2p47abssqrtcos) / (n2 * x1hpx2p47abssqrt.powi(3))
94            }
95            - if x2mx1p47abs <= eps {
96                n0
97            } else {
98                (x1 * x2mx1p47 * x2mx1p47abssqrtcos) / (n2 * x2mx1p47abssqrt.powi(3))
99            },
100    ]
101}
102
103/// Hessian of Eggholder test function.
104///
105/// This function can return [`NaN`](num::Float::nan)-valued elements under the following conditions:
106///
107/// * $|x_1 - x_2 - 47| <= \epsilon$ and $x_1 \neq 0$
108/// * $|x_2 - x_1 + 47| <= \epsilon$ and $x_1 \neq 0$
109/// * $\left|\frac{x_1}{2} + x_2 + 47\right| <= \epsilon$ and $|x_2 + 47| \neq 0$
110///
111/// where $\epsilon$ is `T`'s [epsilon](num::Float::epsilon).
112pub fn eggholder_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
113where
114    T: Float + FromPrimitive,
115{
116    let [x1, x2] = *param;
117
118    let eps = T::epsilon();
119    let n0 = T::from_f64(0.0).unwrap();
120    let n2 = T::from_f64(2.0).unwrap();
121    let n3 = T::from_f64(3.0).unwrap();
122    let n4 = T::from_f64(4.0).unwrap();
123    let n8 = T::from_f64(8.0).unwrap();
124    let n16 = T::from_f64(16.0).unwrap();
125    let n47 = T::from_f64(47.0).unwrap();
126
127    let x1mx2m47 = x1 - x2 - n47;
128    let x1mx2m47abs = x1mx2m47.abs();
129    let x1mx2m47abssqrt = x1mx2m47abs.sqrt();
130    let x1mx2m47abssqrtsin = x1mx2m47abssqrt.sin();
131    let x1mx2m47abssqrtcos = x1mx2m47abssqrt.cos();
132    let x1hpx2p47 = x1 / n2 + x2 + n47;
133    let x1hpx2p47abs = x1hpx2p47.abs();
134    let x1hpx2p47abssqrt = x1hpx2p47abs.sqrt();
135    let x1hpx2p47abssqrtsin = x1hpx2p47abssqrt.sin();
136    let x1hpx2p47abssqrtcos = x1hpx2p47abssqrt.cos();
137    let x2mx1p47 = x2 - x1 + n47;
138    let x2mx1p47abs = x2mx1p47.abs();
139    let x2mx1p47abssqrt = x2mx1p47abs.sqrt();
140    let x2mx1p47abssqrtcos = x2mx1p47abssqrt.cos();
141    let x2mx1p47abssqrtsin = x2mx1p47abssqrt.sin();
142
143    let a = if x1mx2m47abs <= eps {
144        n0
145    } else {
146        (x1 * x1mx2m47abssqrtsin) / (n4 * x1mx2m47abs)
147            - (x1mx2m47 * x1mx2m47abssqrtcos) / (x1mx2m47abssqrt.powi(3))
148            + (n3 * x1 * x1mx2m47.powi(2) * x1mx2m47abssqrtcos) / (n4 * x1mx2m47abssqrt.powi(7))
149    } + if x1mx2m47abs <= eps && x1.abs() <= eps {
150        n0
151    } else {
152        -(x1 * x1mx2m47abssqrtcos) / (n2 * x1mx2m47abssqrt.powi(3))
153    } + if x1hpx2p47abs <= eps {
154        n0
155    } else {
156        (n3 * (x2 + n47) * x1hpx2p47abssqrtcos * x1hpx2p47.powi(2))
157            / (n16 * x1hpx2p47abssqrt.powi(7))
158            + ((x2 + n47) * x1hpx2p47abssqrtsin) / (n16 * x1hpx2p47abs)
159    } - if x1hpx2p47abs <= eps && (x1 + n47).abs() <= eps {
160        n0
161    } else {
162        ((x2 + n47) * x1hpx2p47abssqrtcos) / (n8 * x1hpx2p47abssqrt.powi(3))
163    };
164
165    let b = if x1hpx2p47abs <= eps {
166        n0
167    } else {
168        ((x2 + n47) * x1hpx2p47abssqrtsin) / (n4 * x1hpx2p47abs)
169            - (x1hpx2p47 * x1hpx2p47abssqrtcos) / (x1hpx2p47abssqrt.powi(3))
170            + (n3 * (x2 + n47) * x1hpx2p47.powi(2) * x1hpx2p47abssqrtcos)
171                / (n4 * x1hpx2p47abssqrt.powi(7))
172    } - if x1hpx2p47abs <= eps && (x2 + n47).abs() <= eps {
173        n0
174    } else {
175        ((x2 + n47) * x1hpx2p47abssqrtcos) / (n2 * x1hpx2p47abssqrt.powi(3))
176    } + if x2mx1p47abs <= eps {
177        n0
178    } else {
179        (x1 * x2mx1p47abssqrtsin) / (n4 * x2mx1p47abs)
180            + (n3 * x1 * x2mx1p47.powi(2) * x2mx1p47abssqrtcos) / (n4 * x2mx1p47abssqrt.powi(7))
181    } - if x2mx1p47abs <= eps && x1.abs() <= eps {
182        n0
183    } else {
184        (x1 * x2mx1p47abssqrtcos) / (n2 * x2mx1p47abssqrt.powi(3))
185    };
186
187    let offdiag = if x1hpx2p47abs <= eps {
188        n0
189    } else {
190        ((x2 + n47) * x1hpx2p47abssqrtsin) / (n8 * x1hpx2p47abs)
191            - (x1hpx2p47 * x1hpx2p47abssqrtcos) / (n4 * x1hpx2p47abssqrt.powi(3))
192            + (n3 * (x2 + n47) * x1hpx2p47.powi(2) * x1hpx2p47abssqrtcos)
193                / (n8 * x1hpx2p47abssqrt.powi(7))
194    } - if x1hpx2p47abs <= eps && (x2 + n47).abs() <= eps {
195        n0
196    } else {
197        ((x2 + n47) * x1hpx2p47abssqrtcos) / (n4 * x1hpx2p47abssqrt.powi(3))
198    } + if x2mx1p47abs <= eps {
199        n0
200    } else {
201        -(x1 * x2mx1p47 * x2mx1p47abssqrtsin) / (n4 * x2mx1p47 * x2mx1p47abs)
202            - (x2mx1p47 * x2mx1p47abssqrtcos) / (n2 * x2mx1p47abssqrt.powi(3))
203            - (n3 * x1 * x2mx1p47.powi(2) * x2mx1p47abssqrtcos) / (n4 * x2mx1p47abssqrt.powi(7))
204    } + if x2mx1p47abs <= eps && x1.abs() <= eps {
205        n0
206    } else {
207        (x1 * x2mx1p47abssqrtcos) / (n2 * x2mx1p47abssqrt.powi(3))
208    };
209
210    [[a, offdiag], [offdiag, b]]
211}
212
213#[cfg(test)]
214mod tests {
215    use super::*;
216    use approx::assert_relative_eq;
217    use finitediff::FiniteDiff;
218    use proptest::prelude::*;
219    use std::{f32, f64};
220
221    #[test]
222    fn test_eggholder_optimum() {
223        assert_relative_eq!(
224            eggholder(&[512.0_f32, 404.2319_f32]),
225            -959.6407_f32,
226            epsilon = f32::EPSILON
227        );
228        assert_relative_eq!(
229            eggholder(&[512.0_f64, 404.2319_f64]),
230            -959.6406627106155_f64,
231            epsilon = f64::EPSILON
232        );
233    }
234
235    proptest! {
236        #[test]
237        fn test_eggholder_derivative_finitediff(a in -512.0..512.0, b in -512.0..512.0) {
238            let param = [a, b];
239            let derivative = eggholder_derivative(&param);
240            let derivative_fd = Vec::from(param).central_diff(&|x| eggholder(&[x[0], x[1]]));
241            for i in 0..derivative.len() {
242                assert_relative_eq!(
243                    derivative[i],
244                    derivative_fd[i],
245                    epsilon = 1e-4,
246                    max_relative = 1e-2
247                );
248            }
249        }
250    }
251
252    proptest! {
253        #[test]
254        fn test_eggholder_hessian_finitediff(a in -512.0..512.0, b in -512.0..512.0) {
255            let param = [a, b];
256            let hessian = eggholder_hessian(&param);
257            let hessian_fd =
258                Vec::from(param).central_hessian(&|x| eggholder_derivative(&[x[0], x[1]]).to_vec());
259            let n = hessian.len();
260            // println!("1: {hessian:?} at {a}/{b}");
261            // println!("2: {hessian_fd:?} at {a}/{b}");
262            for i in 0..n {
263                assert_eq!(hessian[i].len(), n);
264                for j in 0..n {
265                    if hessian_fd[i][j].is_finite() {
266                        assert_relative_eq!(
267                            hessian[i][j],
268                            hessian_fd[i][j],
269                            epsilon = 1e-4,
270                            max_relative = 1e-2
271                        );
272                    }
273                }
274            }
275        }
276    }
277}