argmin_testfunctions/
mccorminck.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//! # McCormick test function
9//!
10//! Defined as
11//!
12//! $$
13//! f(x_1,\\,x_2) = \sin(x_1 + x_2) + (x_1 - x_2)^2 - 1.5x_1 + 2.5x_2 + 1
14//! $$
15//!
16//! where $x_1 \in [-1.5,\\,4]$ and $x_2 \in [-3,\\,4]$.
17//!
18//! The global minimum is at $f(x_1,\\,x_2) = f(-0.54719,\\,-1.54719) = -1.913228$.
19
20use num::{Float, FromPrimitive};
21
22/// McCormick test function.
23///
24/// Defined as
25///
26/// $$
27/// f(x_1,\\,x_2) = \sin(x_1 + x_2) + (x_1 - x_2)^2 - 1.5x_1 + 2.5x_2 + 1
28/// $$
29///
30/// where $x_1 \in [-1.5,\\,4]$ and $x_2 \in [-3,\\,4]$.
31///
32/// The global minimum is at $f(x_1,\\,x_2) = f(-0.54719,\\,-1.54719) = -1.913228$.
33pub fn mccorminck<T>(param: &[T; 2]) -> T
34where
35    T: Float + FromPrimitive,
36{
37    let [x1, x2] = *param;
38    (x1 + x2).sin() + (x1 - x2).powi(2) - T::from_f64(1.5).unwrap() * x1
39        + T::from_f64(2.5).unwrap() * x2
40        + T::from_f64(1.0).unwrap()
41}
42
43/// Derivative of McCormick test function.
44pub fn mccorminck_derivative<T>(param: &[T; 2]) -> [T; 2]
45where
46    T: Float + FromPrimitive,
47{
48    let [x1, x2] = *param;
49
50    let n2 = T::from_f64(2.0).unwrap();
51    let n3 = T::from_f64(3.0).unwrap();
52    let n5 = T::from_f64(5.0).unwrap();
53
54    [
55        (x1 + x2).cos() + n2 * (x1 - x2) - n3 / n2,
56        (x1 + x2).cos() - n2 * (x1 - x2) + n5 / n2,
57    ]
58}
59
60/// Hessian of McCormick test function.
61pub fn mccorminck_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
62where
63    T: Float + FromPrimitive,
64{
65    let [x1, x2] = *param;
66
67    let n2 = T::from_f64(2.0).unwrap();
68
69    let a = (x1 + x2).sin();
70
71    let diag = n2 - a;
72    let offdiag = -n2 - a;
73
74    [[diag, offdiag], [offdiag, diag]]
75}
76
77#[cfg(test)]
78mod tests {
79    use super::*;
80    use approx::assert_relative_eq;
81    use finitediff::FiniteDiff;
82    use proptest::prelude::*;
83
84    #[test]
85    fn test_mccorminck_optimum() {
86        assert_relative_eq!(
87            mccorminck(&[-0.54719_f32, -1.54719_f32]),
88            -1.9132228,
89            epsilon = f32::EPSILON
90        );
91        assert_relative_eq!(
92            mccorminck(&[-0.54719_f64, -1.54719_f64]),
93            -1.9132229544882274,
94            epsilon = f32::EPSILON.into()
95        );
96
97        let deriv = mccorminck_derivative(&[-0.54719_f64, -1.54719_f64]);
98        println!("1: {deriv:?}");
99        for i in 0..2 {
100            assert_relative_eq!(deriv[i], 0.0, epsilon = 1e-4);
101        }
102    }
103
104    proptest! {
105        #[test]
106        fn test_mccorminck_derivative_finitediff(a in -1.5..4.0, b in -3.0..4.0) {
107            let param = [a, b];
108            let derivative = mccorminck_derivative(&param);
109            let derivative_fd = Vec::from(param).central_diff(&|x| mccorminck(&[x[0], x[1]]));
110            for i in 0..derivative.len() {
111                assert_relative_eq!(
112                    derivative[i],
113                    derivative_fd[i],
114                    epsilon = 1e-5,
115                    max_relative = 1e-2
116                );
117            }
118        }
119    }
120
121    proptest! {
122        #[test]
123        fn test_mccorminck_hessian_finitediff(a in -1.5..4.0, b in -3.0..4.0) {
124            let param = [a, b];
125            let hessian = mccorminck_hessian(&param);
126            let hessian_fd =
127                Vec::from(param).central_hessian(&|x| mccorminck_derivative(&[x[0], x[1]]).to_vec());
128            let n = hessian.len();
129            // println!("1: {hessian:?} at {a}/{b}");
130            // println!("2: {hessian_fd:?} at {a}/{b}");
131            for i in 0..n {
132                assert_eq!(hessian[i].len(), n);
133                for j in 0..n {
134                    if hessian_fd[i][j].is_finite() {
135                        assert_relative_eq!(
136                            hessian[i][j],
137                            hessian_fd[i][j],
138                            epsilon = 1e-5,
139                            max_relative = 1e-2
140                        );
141                    }
142                }
143            }
144        }
145    }
146}