argmin_testfunctions/
bukin.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//! # Bukin test function No. 6
9//!
10//! Defined as
11//!
12//! $$
13//! f(x_1, x_2) = 100\sqrt{\left|x_2 - 0.01x_1^2\right|} + 0.01\left|x_1 + 10\right|
14//! $$
15//!
16//! where $x_1 \in [-15,\\,-5]$ and $x_2 \in [-3,\\,3]$.
17//!
18//! The global minimum is at $f(x_1, x_2) = f(-10, 1) = 0$.
19
20use num::{Float, FromPrimitive};
21
22/// Bukin test function No. 6
23///
24/// Defined as
25///
26/// $$
27/// f(x_1, x_2) = 100\sqrt{\left|x_2 - 0.01x_1^2\right|} + 0.01\left|x_1 + 10\right|
28/// $$
29///
30/// where $x_1 \in [-15,\\,-5]$ and $x_2 \in [-3,\\,3]$.
31///
32/// The global minimum is at $f(x_1, x_2) = f(-10, 1) = 0$.
33pub fn bukin_n6<T>(param: &[T; 2]) -> T
34where
35    T: Float + FromPrimitive,
36{
37    let [x1, x2] = *param;
38    let n001 = T::from_f64(0.01).unwrap();
39    let n10 = T::from_f64(10.0).unwrap();
40    let n100 = T::from_f64(100.0).unwrap();
41    n100 * (x2 - n001 * x1.powi(2)).abs().sqrt() + n001 * (x1 + n10).abs()
42}
43
44/// Derivative of Bukin test function No. 6
45pub fn bukin_n6_derivative<T>(param: &[T; 2]) -> [T; 2]
46where
47    T: Float + FromPrimitive,
48{
49    let [x1, x2] = *param;
50
51    let n0 = T::from_f64(0.0).unwrap();
52    let n0_01 = T::from_f64(0.01).unwrap();
53    let n10 = T::from_f64(10.0).unwrap();
54    let n50 = T::from_f64(50.0).unwrap();
55    let eps = T::epsilon();
56
57    let denom = (x2 - n0_01 * x1.powi(2)).abs().powi(3).sqrt();
58    let tmp = x2 - n0_01 * x1.powi(2);
59
60    if denom.abs() <= eps {
61        // Deriviative is actually not defined at optimum. Therefore, as soon as we get close,
62        // we'll set the derivative to 0
63        [n0, n0]
64    } else {
65        [
66            n0_01 * (x1 + n10).signum() - x1 * tmp / denom,
67            n50 * tmp / denom,
68        ]
69    }
70}
71
72/// Hessian of Bukin test function No. 6
73pub fn bukin_n6_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
74where
75    T: Float + FromPrimitive,
76{
77    let [x1, x2] = *param;
78
79    let n0 = T::from_f64(0.0).unwrap();
80    let n0_01 = T::from_f64(0.01).unwrap();
81    let n0_02 = T::from_f64(0.02).unwrap();
82    let n0_0001 = T::from_f64(0.0001).unwrap();
83    let n0_5 = T::from_f64(0.5).unwrap();
84    let n25 = T::from_f64(25.0).unwrap();
85    let eps = T::epsilon() * T::from_f64(1e-4).unwrap();
86
87    let tmp = x2 - n0_01 * x1.powi(2);
88    let denom = tmp.abs().powi(7).sqrt();
89
90    if denom.abs() <= eps {
91        // Hessian is actually not defined at optimum. Therefore, as soon as we get close,
92        // we'll set the Hessian to 0
93        [[n0, n0], [n0, n0]]
94    } else {
95        let offdiag = n0_5 * x1 * tmp.powi(2) / denom;
96        [
97            [
98                x2 * (-n0_0001 * x1.powi(4) + n0_02 * x2 * x1.powi(2) - x2.powi(2)) / denom,
99                offdiag,
100            ],
101            [offdiag, -n25 * tmp.powi(2) / denom],
102        ]
103    }
104}
105
106#[cfg(test)]
107mod tests {
108    use super::*;
109    use approx::assert_relative_eq;
110    use finitediff::FiniteDiff;
111    use proptest::prelude::*;
112    use std::{f32, f64};
113
114    #[test]
115    fn test_bukin_n6_optimum() {
116        assert_relative_eq!(bukin_n6(&[-10_f32, 1_f32]), 0.0, epsilon = f32::EPSILON);
117        assert_relative_eq!(bukin_n6(&[-10_f64, 1_f64]), 0.0, epsilon = f64::EPSILON);
118
119        let deriv = bukin_n6_derivative(&[-10_f64, 1_f64]);
120        for i in 0..2 {
121            assert_relative_eq!(deriv[i], 0.0, epsilon = f64::EPSILON);
122        }
123
124        let hessian = bukin_n6_hessian(&[-10_f64, 1_f64]);
125        for i in 0..2 {
126            for j in 0..2 {
127                assert_relative_eq!(hessian[i][j], 0.0, epsilon = f64::EPSILON);
128            }
129        }
130    }
131
132    proptest! {
133        #[test]
134        fn test_bukin_n6_derivative_finitediff(a in -15.0..-5.0, b in -3.0..3.0) {
135            let param = [a, b];
136            let derivative = bukin_n6_derivative(&param);
137            let derivative_fd = Vec::from(param).central_diff(&|x| bukin_n6(&[x[0], x[1]]));
138            for i in 0..derivative.len() {
139                assert_relative_eq!(
140                    derivative[i],
141                    derivative_fd[i],
142                    epsilon = 1e-5,
143                    max_relative = 1e-2
144                );
145            }
146        }
147    }
148
149    proptest! {
150        #[test]
151        fn test_bukin_n6_hessian_finitediff(a in -15.0..-5.0, b in -3.0..3.0) {
152            let param = [a, b];
153            let hessian = bukin_n6_hessian(&param);
154            let hessian_fd = Vec::from(param).central_hessian(&|x| bukin_n6_derivative(&[x[0], x[1]]).to_vec());
155            let n = hessian.len();
156            // println!("1: {a}/{b} {hessian:?}");
157            // println!("2: {a}/{b} {hessian_fd:?}");
158            for i in 0..n {
159                assert_eq!(hessian[i].len(), n);
160                for j in 0..n {
161                    assert_relative_eq!(
162                        hessian[i][j],
163                        hessian_fd[i][j],
164                        epsilon = 1e-5,
165                        max_relative = 1e-2,
166                    );
167                }
168            }
169        }
170    }
171}