argmin_testfunctions/
beale.rs1use num::{Float, FromPrimitive};
22
23pub 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
45pub 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
83pub 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(¶m);
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(¶m);
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}