argmin_testfunctions/
schaffer.rs1use num::{Float, FromPrimitive};
34
35pub fn schaffer_n2<T>(param: &[T; 2]) -> T
47where
48 T: Float + FromPrimitive,
49{
50 let [x1, x2] = *param;
51
52 let n0_001 = T::from_f64(0.001).unwrap();
53 let n0_5 = T::from_f64(0.5).unwrap();
54 let n1 = T::from_f64(1.0).unwrap();
55
56 n0_5 + ((x1.powi(2) - x2.powi(2)).sin().powi(2) - n0_5)
57 / (n1 + n0_001 * (x1.powi(2) + x2.powi(2))).powi(2)
58}
59
60pub fn schaffer_n2_derivative<T>(param: &[T; 2]) -> [T; 2]
62where
63 T: Float + FromPrimitive,
64{
65 let [x1, x2] = *param;
66
67 let n0_001 = T::from_f64(0.001).unwrap();
68 let n0_004 = T::from_f64(0.004).unwrap();
69 let n0_5 = T::from_f64(0.5).unwrap();
70 let n1 = T::from_f64(1.0).unwrap();
71 let n4 = T::from_f64(4.0).unwrap();
72
73 let x1spx2s = x1.powi(2) + x2.powi(2);
74 let x1smx2s = x1.powi(2) - x2.powi(2);
75 let x2smx1s = x2.powi(2) - x1.powi(2);
76 let tmp = n0_001 * x1spx2s + n1;
77 let denom2 = tmp.powi(2);
78 let denom3 = tmp.powi(3);
79 let a = x1smx2s.sin() * x1smx2s.cos();
80 let a2 = x2smx1s.sin() * x2smx1s.cos();
81 let b = n0_004 * (x1smx2s.sin().powi(2) - n0_5);
82
83 [
84 (n4 * x1 * a) / denom2 - (x1 * b) / denom3,
85 (n4 * x2 * a2) / denom2 - (x2 * b) / denom3,
86 ]
87}
88
89pub fn schaffer_n2_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
91where
92 T: Float + FromPrimitive,
93{
94 let [x1, x2] = *param;
95
96 let n0_001 = T::from_f64(0.001).unwrap();
97 let n0_004 = T::from_f64(0.004).unwrap();
98 let n0_006 = T::from_f64(0.006).unwrap();
99 let n0_032 = T::from_f64(0.032).unwrap();
100 let n0_5 = T::from_f64(0.5).unwrap();
101 let n1 = T::from_f64(1.0).unwrap();
102 let n4 = T::from_f64(4.0).unwrap();
103 let n8 = T::from_f64(8.0).unwrap();
104
105 let x1s = x1.powi(2);
106 let x2s = x2.powi(2);
107 let x1spx2s = x1s + x2s;
108 let x1smx2s = x1s - x2s;
109 let x2smx1s = x2s - x1s;
110 let x1smx2ssin2 = x1smx2s.sin().powi(2);
111 let x2smx1ssin2 = x2smx1s.sin().powi(2);
112 let x1smx2scos2 = x1smx2s.cos().powi(2);
113 let x2smx1scos2 = x2smx1s.cos().powi(2);
114 let tmp = n0_001 * x1spx2s + n1;
115 let denom2 = tmp.powi(2);
116 let denom3 = tmp.powi(3);
117 let denom4 = tmp.powi(4);
118 let a = x1smx2s.sin() * x1smx2s.cos();
119 let a2 = x2smx1s.sin() * x2smx1s.cos();
120 let b = n0_004 * (x1smx2ssin2 - n0_5);
121
122 let offdiag =
123 (n8 * x1 * x2 * (x1smx2ssin2 - x1smx2scos2)) / denom2 + (n0_006 * x1 * x2 * b) / denom4;
124
125 [
126 [
127 (-(n8 * x1s * x1smx2ssin2) + (n4 * a) + (n8 * x1s * x1smx2scos2)) / denom2
128 - (b + (n0_032 * x1s * a)) / denom3
129 + (n0_006 * x1s * b) / denom4,
130 offdiag,
131 ],
132 [
133 offdiag,
134 (-(n8 * x2s * x2smx1ssin2) + (n4 * a2) + (n8 * x2s * x2smx1scos2)) / denom2
135 - (b + (n0_032 * x2s * a2)) / denom3
136 + (n0_006 * x2s * b) / denom4,
137 ],
138 ]
139}
140
141pub fn schaffer_n4<T>(param: &[T; 2]) -> T
153where
154 T: Float + FromPrimitive,
155{
156 let [x1, x2] = *param;
157 let n05 = T::from_f64(0.5).unwrap();
158 let n1 = T::from_f64(1.0).unwrap();
159 let n0001 = T::from_f64(0.001).unwrap();
160 n05 + ((x1.powi(2) - x2.powi(2)).abs().sin().cos().powi(2) - n05)
161 / (n1 + n0001 * (x1.powi(2) + x2.powi(2))).powi(2)
162}
163
164pub fn schaffer_n4_derivative<T>(param: &[T; 2]) -> [T; 2]
166where
167 T: Float + FromPrimitive,
168{
169 let [x1, x2] = *param;
170 let n0_5 = T::from_f64(0.5).unwrap();
171 let n1 = T::from_f64(1.0).unwrap();
172 let n0_001 = T::from_f64(0.001).unwrap();
173 let n0_004 = T::from_f64(0.004).unwrap();
174 let n4 = T::from_f64(4.0).unwrap();
175
176 let x1smx2s = x1.powi(2) - x2.powi(2);
177 let x2smx1s = x2.powi(2) - x1.powi(2);
178 let x1spx2s = x1.powi(2) + x2.powi(2);
179 let x1smx2scos = x1smx2s.cos();
180 let x2smx1scos = x2smx1s.cos();
181 let x1smx2sabs = x1smx2s.abs();
182 let x1smx2sabssin = x1smx2sabs.sin();
183 let x2smx1sabssin = x1smx2sabs.sin();
184 let x1smx2sabssincos = x1smx2sabssin.cos();
185 let x1smx2sabssinsin = x1smx2sabssin.sin();
186 let x1smx2sabssincos2 = x1smx2sabssin.cos().powi(2);
187 let x2smx1sabssincos = x2smx1sabssin.cos();
188 let x2smx1sabssinsin = x2smx1sabssin.sin();
189 let x2smx1sabssincos2 = x2smx1sabssin.cos().powi(2);
190 let denom_a = (n0_001 * x1spx2s + n1).powi(2);
191 let denom_b = (n0_001 * x1spx2s + n1).powi(3);
192
193 [
194 -(n4 * x1 * x1smx2s * x1smx2scos * x1smx2sabssincos * x1smx2sabssinsin)
195 / (denom_a * x1smx2sabs)
196 - (n0_004 * x1 * (x1smx2sabssincos2 - n0_5)) / denom_b,
197 -(n4 * x2 * x2smx1s * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin)
198 / (denom_a * x1smx2sabs)
199 - (n0_004 * x2 * (x2smx1sabssincos2 - n0_5)) / denom_b,
200 ]
201}
202
203pub fn schaffer_n4_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
205where
206 T: Float + FromPrimitive,
207{
208 let [x1, x2] = *param;
209
210 let eps = T::epsilon();
211 let n0_5 = T::from_f64(0.5).unwrap();
212 let n0 = T::from_f64(0.0).unwrap();
213 let n1 = T::from_f64(1.0).unwrap();
214 let n0_001 = T::from_f64(0.001).unwrap();
215 let n0_004 = T::from_f64(0.004).unwrap();
216 let n0_006 = T::from_f64(0.006).unwrap();
217 let n0_016 = T::from_f64(0.016).unwrap();
218 let n0_032 = T::from_f64(0.032).unwrap();
219 let n4 = T::from_f64(4.0).unwrap();
220 let n8 = T::from_f64(8.0).unwrap();
221
222 let x1s = x1.powi(2);
223 let x2s = x2.powi(2);
224 let x1smx2s = x1s - x2s;
225 let x2smx1s = x2s - x1s;
226 let x1spx2s = x1s + x2s;
227 let x1smx2scos = x1smx2s.cos();
228 let x1smx2ssin = x1smx2s.sin();
229 let x2smx1ssin = x2smx1s.sin();
230 let x2smx1scos = x2smx1s.cos();
231 let x1smx2scos2 = x1smx2scos.powi(2);
232 let x2smx1scos2 = x2smx1scos.powi(2);
233 let x1smx2sabs = x1smx2s.abs();
234 let x1smx2sabssin = x1smx2sabs.sin();
235 let x2smx1sabssin = x1smx2sabs.sin();
236 let x1smx2sabssincos = x1smx2sabssin.cos();
237 let x1smx2sabssinsin = x1smx2sabssin.sin();
238 let x2smx1sabssinsin = x2smx1sabssin.sin();
239 let x1smx2sabssinsin2 = x1smx2sabssinsin.powi(2);
240 let x2smx1sabssinsin2 = x2smx1sabssinsin.powi(2);
241 let x1smx2sabssincos2 = x1smx2sabssin.cos().powi(2);
242 let x2smx1sabssincos = x2smx1sabssin.cos();
243 let x2smx1sabssinsin = x2smx1sabssin.sin();
244 let x2smx1sabssincos2 = x2smx1sabssin.cos().powi(2);
245 let denom_a = (n0_001 * x1spx2s + n1).powi(2);
246 let denom_b = (n0_001 * x1spx2s + n1).powi(3);
247 let denom_c = (n0_001 * x1spx2s + n1).powi(4);
248
249 if x1smx2sabs <= eps {
250 [[n0, n0], [n0, n0]]
251 } else {
252 let a = (n8 * x1s * x1smx2scos2 * x1smx2sabssinsin2
253 - n8 * x1s * x1smx2scos2 * x1smx2sabssincos2)
254 / denom_a
255 + ((n8 * x1s * x1smx2s * x1smx2ssin * x1smx2sabssincos * x1smx2sabssinsin)
256 - (n4 * x1smx2s * x1smx2scos * x1smx2sabssincos * x1smx2sabssinsin))
257 / (denom_a * x1smx2sabs)
258 + (n0_032 * x1s * x1smx2s * x1smx2scos * x1smx2sabssincos * x1smx2sabssinsin)
259 / (denom_b * x1smx2sabs)
260 + (-n0_004 * (x1smx2sabssincos2 - n0_5)) / denom_b
261 + (n0_006 * n0_004 * x1s * (x1smx2sabssincos2 - n0_5)) / denom_c;
262
263 let b = (n8 * x2s * x2smx1scos2 * x2smx1sabssinsin2
264 - n8 * x2s * x2smx1scos2 * x2smx1sabssincos2)
265 / denom_a
266 + ((n8 * x2s * x2smx1s * x2smx1ssin * x2smx1sabssincos * x2smx1sabssinsin)
267 - (n4 * x2smx1s * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin))
268 / (denom_a * x1smx2sabs)
269 + (n0_032 * x2s * x2smx1s * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin)
270 / (denom_b * x1smx2sabs)
271 + (-n0_004 * (x2smx1sabssincos2 - n0_5)) / denom_b
272 + (n0_006 * n0_004 * x2s * (x2smx1sabssincos2 - n0_5)) / denom_c;
273
274 let offdiag = (n8 * x1 * x2 * x1smx2s * x2smx1scos2 * x2smx1sabssinsin2)
275 / (x2smx1s * denom_a)
276 + (n8 * x1 * x2 * x1smx2s * x2smx1ssin * x2smx1sabssincos * x2smx1sabssinsin)
277 / (x1smx2sabs * denom_a)
278 + (n8 * x1 * x2 * x1smx2s * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin)
279 / (x2smx1s * denom_a * x1smx2sabs)
280 + (n8 * x1 * x2 * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin)
281 / (denom_a * x1smx2sabs)
282 + (n0_016 * x1 * x2 * x2smx1s * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin)
283 / (denom_b * x1smx2sabs)
284 + (n0_016 * x1 * x2 * x1smx2s * x2smx1scos * x2smx1sabssincos * x2smx1sabssinsin)
285 / (denom_b * x1smx2sabs)
286 + (-n8 * x1 * x2 * x1smx2s * x2smx1scos2 * x2smx1sabssincos2) / (denom_a * x2smx1s)
287 + (n0_006 * n0_004 * x1 * x2 * (x2smx1sabssincos2 - n0_5)) / denom_c;
288
289 [[a, offdiag], [offdiag, b]]
290 }
291}
292
293#[cfg(test)]
294mod tests {
295 use super::*;
296 use approx::assert_relative_eq;
297 use finitediff::FiniteDiff;
298 use proptest::prelude::*;
299 use std::{f32, f64};
300
301 #[test]
302 fn test_schaffer_n2_optimum() {
303 assert_relative_eq!(schaffer_n2(&[0_f32, 0_f32]), 0.0, epsilon = f32::EPSILON);
304 assert_relative_eq!(schaffer_n2(&[0_f64, 0_f64]), 0.0, epsilon = f64::EPSILON);
305
306 let deriv = schaffer_n2_derivative(&[0.0, 0.0]);
307 for i in 0..2 {
308 assert_relative_eq!(deriv[i], 0.0, epsilon = f64::EPSILON);
309 }
310 }
311
312 proptest! {
313 #[test]
314 fn test_schaffer_n2_derivative_finitediff(a in -10.0..10.0, b in -10.0..10.0) {
315 let param = [a, b];
318 let derivative = schaffer_n2_derivative(¶m);
319 let derivative_fd = Vec::from(param).central_diff(&|x| schaffer_n2(&[x[0], x[1]]));
320 for i in 0..derivative.len() {
323 assert_relative_eq!(
324 derivative[i],
325 derivative_fd[i],
326 epsilon = 1e-3,
327 max_relative = 1e-2
328 );
329 }
330 }
331 }
332
333 proptest! {
334 #[test]
335 fn test_schaffer_n2_hessian_finitediff(a in -10.0..10.0, b in -10.0..10.0) {
336 let param = [a, b];
339 let hessian = schaffer_n2_hessian(¶m);
340 let hessian_fd =
341 Vec::from(param).forward_hessian(&|x| schaffer_n2_derivative(&[x[0], x[1]]).to_vec());
342 let n = hessian.len();
343 for i in 0..n {
346 assert_eq!(hessian[i].len(), n);
347 for j in 0..n {
348 assert_relative_eq!(
349 hessian[i][j],
350 hessian_fd[i][j],
351 epsilon = 1e-3,
352 max_relative = 1e-2
353 );
354 }
355 }
356 }
357 }
358
359 #[test]
360 fn test_schaffer_n4_optimum() {
361 assert_relative_eq!(
362 schaffer_n4(&[0_f32, 1.25313_f32]),
363 0.29257864,
364 epsilon = f32::EPSILON
365 );
366
367 let deriv = schaffer_n4_derivative(&[0.0, 1.25313]);
368 for i in 0..2 {
369 assert_relative_eq!(deriv[i], 0.0, epsilon = 1e-4);
370 }
371 }
372
373 proptest! {
374 #[test]
375 fn test_schaffer_n4_derivative_finitediff(a in -100.0..100.0, b in -100.0..100.0) {
376 let param = [a, b];
377 let derivative = schaffer_n4_derivative(¶m);
378 let derivative_fd = Vec::from(param).central_diff(&|x| schaffer_n4(&[x[0], x[1]]));
379 for i in 0..derivative.len() {
382 assert_relative_eq!(
383 derivative[i],
384 derivative_fd[i],
385 epsilon = 1e-3,
386 max_relative = 1e-2,
387 );
388 }
389 }
390 }
391
392 proptest! {
393 #[test]
394 fn test_schaffer_n4_hessian_finitediff(a in -10.0..10.0, b in -10.0..10.0) {
395 let param = [a, b];
398 let hessian = schaffer_n4_hessian(¶m);
399 let hessian_fd =
400 Vec::from(param).forward_hessian(&|x| schaffer_n4_derivative(&[x[0], x[1]]).to_vec());
401 let n = hessian.len();
402 for i in 0..n {
405 assert_eq!(hessian[i].len(), n);
406 for j in 0..n {
407 if hessian_fd[i][j].is_finite() {
408 assert_relative_eq!(
409 hessian[i][j],
410 hessian_fd[i][j],
411 epsilon = 1e-3,
412 max_relative = 1e-2
413 );
414 }
415 }
416 }
417 }
418 }
419}