argmin_testfunctions/
schaffer.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//! # Schaffer test function No. 2
9//!
10//! Defined as
11//!
12//! $$
13//! f(x_1,\\,x_2) = 0.5 + \frac{\sin^2(x_1^2 - x_2^2) - 0.5}{\left[1 + 0.001(x_1^2 + x_2^2)\right]^2}
14//! $$
15//!
16//! where $x_i \in [-100, 100]$.
17//!
18//! The global minimum is at $f(x_1,\\,x_2) = f(0,\\,0) = 0$.
19//!
20//! # Schaffer test function No. 4
21//!
22
23//! Defined as
24//!
25//! $$
26//! f(x_1,\\,x_2) = 0.5 + \frac{\cos^2(\sin(\left|x_1^2 - x_2^2\right|)) - 0.5}{\left[1 + 0.001(x_1^2 + x_2^2)\right]^2}
27//! $$
28//!
29//! where $x_i \in [-100, 100]$.
30//!
31//! The global minimum is at $f(x_1,\\,x_2) = f(0,\\,1.25313) = 0.291992$.
32
33use num::{Float, FromPrimitive};
34
35/// Schaffer test function No. 2.
36///
37/// Defined as
38///
39/// $$
40/// f(x_1,\\,x_2) = 0.5 + \frac{\sin^2(x_1^2 - x_2^2) - 0.5}{\left[1 + 0.001(x_1^2 + x_2^2)\right]^2}
41/// $$
42///
43/// where $x_i \in [-100, 100]$.
44///
45/// The global minimum is at $f(x_1,\\,x_2) = f(0,\\,0) = 0$.
46pub 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
60/// Derivative of Schaffer test function No. 2
61pub 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
89/// Hessian of Schaffer test function No. 2
90pub 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
141/// Schaffer test function No. 4
142///
143/// Defined as
144///
145/// $$
146/// f(x_1,\\,x_2) = 0.5 + \frac{\cos^2(\sin(\left|x_1^2 - x_2^2\right|)) - 0.5}{\left[1 + 0.001(x_1^2 + x_2^2)\right]^2}
147/// $$
148///
149/// where $x_i \in [-100, 100]$.
150///
151/// The global minimum is at $f(x_1,\\,x_2) = f(0,\\,1.25313) = 0.291992$.
152pub 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
164/// Derivative of Schaffer test function No. 4
165pub 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
203/// Hessian of Schaffer test function No. 4
204pub 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            // Note: prop testing range is smaller than range of function, simply because finitediff
316            // has huge errors far away from the optimum.
317            let param = [a, b];
318            let derivative = schaffer_n2_derivative(&param);
319            let derivative_fd = Vec::from(param).central_diff(&|x| schaffer_n2(&[x[0], x[1]]));
320            // println!("1: {derivative:?} at {a}/{b}");
321            // println!("2: {derivative_fd:?} at {a}/{b}");
322            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            // Note: prop testing range is smaller than range of function, simply because finitediff
337            // has huge errors far away from the optimum.
338            let param = [a, b];
339            let hessian = schaffer_n2_hessian(&param);
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            // println!("1: {hessian:?} at {a}/{b}");
344            // println!("2: {hessian_fd:?} at {a}/{b}");
345            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(&param);
378            let derivative_fd = Vec::from(param).central_diff(&|x| schaffer_n4(&[x[0], x[1]]));
379            // println!("1: {derivative:?} at {a}/{b}");
380            // println!("2: {derivative_fd:?} at {a}/{b}");
381            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            // Note: prop testing range is smaller than range of function, simply because finitediff
396            // has huge errors far away from the optimum.
397            let param = [a, b];
398            let hessian = schaffer_n4_hessian(&param);
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            // println!("1: {hessian:?} at {a}/{b}");
403            // println!("2: {hessian_fd:?} at {a}/{b}");
404            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}