argmin_testfunctions/
eggholder.rs1use num::{Float, FromPrimitive};
22
23pub 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
49pub 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
103pub 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(¶m);
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(¶m);
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 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}