argmin_testfunctions/
beale.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//! # Beale test function
9//!
10//! Defined as
11//!
12//! $$
13//! f(x_1, x_2) = (1.5 - x_1 + x_1 x_2)^2 + (2.25 - x_1 + x_1 x_2^2)^2 +
14//!               (2.625 - x_1 + x_1 x_2^3)^2
15//! $$
16//!
17//! where $x_i \in [-4.5,\\,4.5]$.
18//!
19//! The global minimum is at $f(x_1, x_2) = f(3, 0.5) = 0$.
20
21use num::{Float, FromPrimitive};
22
23/// Beale test function
24///
25/// Defined as
26///
27/// $$
28/// f(x_1, x_2) = (1.5 - x_1 + x_1 x_2)^2 + (2.25 - x_1 + x_1 x_2^2)^2 +
29///               (2.625 - x_1 + x_1 x_2^3)^2
30/// $$
31///
32/// where $x_i \in [-4.5,\\,4.5]$.
33///
34/// The global minimum is at $f(x_1, x_2) = f(3, 0.5) = 0$.
35pub fn beale<T>(param: &[T; 2]) -> T
36where
37    T: Float + FromPrimitive,
38{
39    let [x1, x2] = *param;
40    (T::from_f64(1.5).unwrap() - x1 + x1 * x2).powi(2)
41        + (T::from_f64(2.25).unwrap() - x1 + x1 * (x2.powi(2))).powi(2)
42        + (T::from_f64(2.625).unwrap() - x1 + x1 * (x2.powi(3))).powi(2)
43}
44
45/// Derivative of Beale test function
46pub fn beale_derivative<T>(param: &[T; 2]) -> [T; 2]
47where
48    T: Float + FromPrimitive,
49{
50    let n0_5 = T::from_f64(0.5).unwrap();
51    let n1 = T::from_f64(1.0).unwrap();
52    let n1_5 = T::from_f64(1.5).unwrap();
53    let n2 = T::from_f64(2.0).unwrap();
54    let n2_625 = T::from_f64(2.625).unwrap();
55    let n3 = T::from_f64(3.0).unwrap();
56    let n4_5 = T::from_f64(4.5).unwrap();
57    let n5_25 = T::from_f64(5.25).unwrap();
58    let n6 = T::from_f64(6.0).unwrap();
59    let n12_75 = T::from_f64(12.75).unwrap();
60
61    let [x1, x2] = *param;
62
63    let x2p2 = x2.powi(2);
64    let x2p3 = x2.powi(3);
65    let x2p4 = x2.powi(4);
66    let x2p5 = x2.powi(5);
67    let x2p6 = x2.powi(6);
68
69    [
70        n5_25 * x2p3
71            + n4_5 * x2p2
72            + n3 * x2
73            + n2 * x1 * (x2p6 + x2p4 - n2 * x2p3 - x2p2 - n2 * x2 + n3)
74            - n12_75,
75        n6 * x1
76            * (n2_625 * x2p2
77                + n1_5 * x2
78                + x1 * (x2p5 + (n2 / n3) * x2p3 - x2p2 - (n1 / n3) * (x2 + n1))
79                + n0_5),
80    ]
81}
82
83/// Hessian of Beale test function
84pub fn beale_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
85where
86    T: Float + FromPrimitive,
87{
88    let n2 = T::from_f64(2.0).unwrap();
89    let n3 = T::from_f64(3.0).unwrap();
90    let n4 = T::from_f64(4.0).unwrap();
91    let n8 = T::from_f64(8.0).unwrap();
92    let n9 = T::from_f64(9.0).unwrap();
93    let n12 = T::from_f64(12.0).unwrap();
94    let n15_75 = T::from_f64(15.75).unwrap();
95    let n30 = T::from_f64(30.0).unwrap();
96    let n31_5 = T::from_f64(31.5).unwrap();
97
98    let [x1, x2] = *param;
99
100    let x1p2 = x1.powi(2);
101    let x2p2 = x2.powi(2);
102    let x2p3 = x2.powi(3);
103    let x2p4 = x2.powi(4);
104    let x2p5 = x2.powi(5);
105    let x2p6 = x2.powi(6);
106
107    let offdiag = n12 * x1 * x2p5 + n8 * x1 * x2p3 - n12 * x1 * x2p2 + n15_75 * x2p2 - n4 * x1 * x2
108        + n9 * x2
109        - n4 * x1
110        + n3;
111
112    [
113        [
114            n2 * (x2p6 + x2p4 - n2 * x2p3 - x2p2 - n2 * x2 + n3),
115            offdiag,
116        ],
117        [
118            offdiag,
119            n30 * x1p2 * x2p4 + n12 * x1p2 * x2p2 - n12 * x1p2 * x2 + n31_5 * x1 * x2 - n2 * x1p2
120                + n9 * x1,
121        ],
122    ]
123}
124
125#[cfg(test)]
126mod tests {
127    use super::*;
128    use approx::assert_relative_eq;
129    use finitediff::FiniteDiff;
130    use proptest::prelude::*;
131    use std::{f32, f64};
132
133    #[test]
134    fn test_beale_optimum() {
135        assert_relative_eq!(beale(&[3.0_f32, 0.5_f32]), 0.0, epsilon = f32::EPSILON);
136        assert_relative_eq!(beale(&[3.0_f64, 0.5_f64]), 0.0, epsilon = f64::EPSILON);
137    }
138
139    proptest! {
140        #[test]
141        fn test_beale_derivative(a in -4.5..4.5, b in -4.5..4.5) {
142            let param = [a, b];
143            let derivative = beale_derivative(&param);
144            let derivative_fd = Vec::from(param).central_diff(&|x| beale(&[x[0], x[1]]));
145            for i in 0..derivative.len() {
146                assert_relative_eq!(derivative[i], derivative_fd[i], epsilon = 1e-2);
147            }
148        }
149    }
150
151    proptest! {
152        #[test]
153        fn test_beale_hessian_finitediff(a in -4.5..4.5, b in -4.5..4.5) {
154            let param = [a, b];
155            let hessian = beale_hessian(&param);
156            let hessian_fd =
157                Vec::from(param).forward_hessian(&|x| beale_derivative(&[x[0], x[1]]).to_vec());
158            let n = hessian.len();
159            for i in 0..n {
160                assert_eq!(hessian[i].len(), n);
161                for j in 0..n {
162                    assert_relative_eq!(
163                        hessian[i][j],
164                        hessian_fd[i][j],
165                        epsilon = 1e-5,
166                        max_relative = 1e-2
167                    );
168                }
169            }
170        }
171    }
172}