argmin_testfunctions/
picheny.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//! # Picheny test function
9//!
10//! Variation of the Goldstein-Price test function.
11//!
12//! Defined as
13//!
14//! $$
15//! f(x_1,\\,x_2) = \frac{1}{2.427}\biggl[
16//!                    \log\Bigl[\bigl(1 + (\bar{x}_1 + \bar{x}_2 + 1)^2(19 - 14\bar{x}_1 +
17//!                    3\bar{x}_1^2 - 14\bar{x}_2 + 6\bar{x}_1\bar{x}_2 + 3\bar{x}_2^2)\bigr) \\\\
18//!                    \times\bigl(30 + (2\bar{x}_1 - 3\bar{x}_2)^2(18 - 32\bar{x}_1 + 12\bar{x}_1^2 +
19//!                    48\bar{x}_2 - 36\bar{x}_1\bar{x}_2 + 27\bar{x}_2^2)\bigr)\Bigr]
20//!                 - 8.693\biggr]
21//! $$
22//!
23//! where $\bar{x}_i = 4x_i - 2$ and $x_i \in [0, 1]$.
24//!
25//! The global minimum is at $f(x_1,\\,x_2) = f(0.5,\\,0.25) = 3.3851993182036826$.
26
27use num::{Float, FromPrimitive};
28
29/// Picheny test function.
30///
31/// Variation on the [Goldstein-Price test function](crate::goldsteinprice::goldsteinprice) introduced by [Picheny et al. (2012)](https://hal.science/hal-00658212).
32///
33/// Defined as
34///
35/// $$
36/// f(x_1,\\,x_2) = \frac{1}{2.427}\biggl[
37///                    \log\Bigl[\bigl(1 + (\bar{x}_1 + \bar{x}_2 + 1)^2(19 - 14\bar{x}_1 +
38///                    3\bar{x}_1^2 - 14\bar{x}_2 + 6\bar{x}_1\bar{x}_2 + 3\bar{x}_2^2)\bigr) \\\\
39///                    \times\bigl(30 + (2\bar{x}_1 - 3\bar{x}_2)^2(18 - 32\bar{x}_1 + 12\bar{x}_1^2 +
40///                    48\bar{x}_2 - 36\bar{x}_1\bar{x}_2 + 27\bar{x}_2^2)\bigr)\Bigr]
41///                 - 8.693\biggr]
42/// $$
43///
44/// where $\bar{x}_i = 4x_i - 2$ and $x_i \in [0, 1]$.
45///
46/// The global minimum is at $f(x_1,\\,x_2) = f(0.5,\\,0.25) = 3.3851993182036826$.
47pub fn picheny<T>(param: &[T; 2]) -> T
48where
49    T: Float + FromPrimitive,
50{
51    let n1 = T::from_f64(1.0).unwrap();
52    let n2 = T::from_f64(2.0).unwrap();
53    let n3 = T::from_f64(3.0).unwrap();
54    let n4 = T::from_f64(4.0).unwrap();
55    let n6 = T::from_f64(6.0).unwrap();
56    let n12 = T::from_f64(12.0).unwrap();
57    let n14 = T::from_f64(14.0).unwrap();
58    let n18 = T::from_f64(18.0).unwrap();
59    let n19 = T::from_f64(19.0).unwrap();
60    let n27 = T::from_f64(27.0).unwrap();
61    let n30 = T::from_f64(30.0).unwrap();
62    let n32 = T::from_f64(32.0).unwrap();
63    let n36 = T::from_f64(36.0).unwrap();
64    let n48 = T::from_f64(48.0).unwrap();
65
66    let [x1, x2] = param.map(|x| n4 * x - n2);
67
68    T::from_f64(1.0 / 2.427).unwrap()
69        * (((n1
70            + (x1 + x2 + n1).powi(2)
71                * (n19 - n14 * (x1 + x2) + n3 * (x1.powi(2) + x2.powi(2)) + n6 * x1 * x2))
72            * (n30
73                + (n2 * x1 - n3 * x2).powi(2)
74                    * (n18 - n32 * x1 + n12 * x1.powi(2) + n48 * x2 - n36 * x1 * x2
75                        + n27 * x2.powi(2))))
76        .log10()
77            - T::from_f64(8.693).unwrap())
78}
79/// Derivative of Picheny test function.
80pub fn picheny_derivative<T>(param: &[T; 2]) -> [T; 2]
81where
82    T: Float + FromPrimitive,
83{
84    let n1 = T::from_f64(1.0).unwrap();
85    let n2 = T::from_f64(2.0).unwrap();
86    let n3 = T::from_f64(3.0).unwrap();
87    let n4 = T::from_f64(4.0).unwrap();
88    let n6 = T::from_f64(6.0).unwrap();
89    let n8 = T::from_f64(8.0).unwrap();
90    let n10 = T::from_f64(10.0).unwrap();
91    let n12 = T::from_f64(12.0).unwrap();
92    let n14 = T::from_f64(14.0).unwrap();
93    let n16 = T::from_f64(16.0).unwrap();
94    let n18 = T::from_f64(18.0).unwrap();
95    let n19 = T::from_f64(19.0).unwrap();
96    let n24 = T::from_f64(24.0).unwrap();
97    let n27 = T::from_f64(27.0).unwrap();
98    let n30 = T::from_f64(30.0).unwrap();
99    let n32 = T::from_f64(32.0).unwrap();
100    let n36 = T::from_f64(36.0).unwrap();
101    let n48 = T::from_f64(48.0).unwrap();
102    let n56 = T::from_f64(56.0).unwrap();
103    let n96 = T::from_f64(96.0).unwrap();
104    let n128 = T::from_f64(128.0).unwrap();
105    let n144 = T::from_f64(144.0).unwrap();
106    let n192 = T::from_f64(192.0).unwrap();
107    let n216 = T::from_f64(216.0).unwrap();
108    let factor = T::from_f64(1.0 / 2.427).unwrap();
109
110    let [x1, x2] = *param;
111
112    let a = (factor
113        * ((n8
114            * (n4 * x1 + n4 * x2 - n3)
115            * (n3 * (n4 * x1 - n2).powi(2) + n6 * (n4 * x2 - n2) * (n4 * x1 - n2)
116                - n14 * (n4 * x1 - n2)
117                + n3 * (n4 * x2 - n2).powi(2)
118                - n14 * (n4 * x2 - n2)
119                + n19)
120            + (n4 * x1 + n4 * x2 - n3).powi(2)
121                * (n24 * (n4 * x1 - n2) + n24 * (n4 * x2 - n2) - n56))
122            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
123                * (n12 * (n4 * x1 - n2).powi(2)
124                    - n36 * (n4 * x2 - n2) * (n4 * x1 - n2)
125                    - n32 * (n4 * x1 - n2)
126                    + n27 * (n4 * x2 - n2).powi(2)
127                    + n48 * (n4 * x2 - n2)
128                    + n18)
129                + n30)
130            + ((n4 * x1 + n4 * x2 - n3).powi(2)
131                * (n3 * (n4 * x1 - n2).powi(2) + n6 * (n4 * x2 - n2) * (n4 * x1 - n2)
132                    - n14 * (n4 * x1 - n2)
133                    + n3 * (n4 * x2 - n2).powi(2)
134                    - n14 * (n4 * x2 - n2)
135                    + n19)
136                + n1)
137                * (n16
138                    * (n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2))
139                    * (n12 * (n4 * x1 - n2).powi(2)
140                        - n36 * (n4 * x2 - n2) * (n4 * x1 - n2)
141                        - n32 * (n4 * x1 - n2)
142                        + n27 * (n4 * x2 - n2).powi(2)
143                        + n48 * (n4 * x2 - n2)
144                        + n18)
145                    + (n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
146                        * (n96 * (n4 * x1 - n2) - n144 * (n4 * x2 - n2) - n128))))
147        / (n10.ln()
148            * ((n4 * x1 + n4 * x2 - n3).powi(2)
149                * (n3 * (n4 * x1 - n2).powi(2) + n6 * (n4 * x2 - n2) * (n4 * x1 - n2)
150                    - n14 * (n4 * x1 - n2)
151                    + n3 * (n4 * x2 - n2).powi(2)
152                    - n14 * (n4 * x2 - n2)
153                    + n19)
154                + n1)
155            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
156                * (n12 * (n4 * x1 - n2).powi(2)
157                    - n36 * (n4 * x2 - n2) * (n4 * x1 - n2)
158                    - n32 * (n4 * x1 - n2)
159                    + n27 * (n4 * x2 - n2).powi(2)
160                    + n48 * (n4 * x2 - n2)
161                    + n18)
162                + n30));
163
164    let b = (factor
165        * ((n8
166            * (n4 * x2 + n4 * x1 - n3)
167            * (n3 * (n4 * x2 - n2).powi(2) + n6 * (n4 * x1 - n2) * (n4 * x2 - n2)
168                - n14 * (n4 * x2 - n2)
169                + n3 * (n4 * x1 - n2).powi(2)
170                - n14 * (n4 * x1 - n2)
171                + n19)
172            + (n4 * x2 + n4 * x1 - n3).powi(2)
173                * (n24 * (n4 * x2 - n2) + n24 * (n4 * x1 - n2) - n56))
174            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
175                * (n27 * (n4 * x2 - n2).powi(2) - n36 * (n4 * x1 - n2) * (n4 * x2 - n2)
176                    + n48 * (n4 * x2 - n2)
177                    + n12 * (n4 * x1 - n2).powi(2)
178                    - n32 * (n4 * x1 - n2)
179                    + n18)
180                + n30)
181            + ((n4 * x2 + n4 * x1 - n3).powi(2)
182                * (n3 * (n4 * x2 - n2).powi(2) + n6 * (n4 * x1 - n2) * (n4 * x2 - n2)
183                    - n14 * (n4 * x2 - n2)
184                    + n3 * (n4 * x1 - n2).powi(2)
185                    - n14 * (n4 * x1 - n2)
186                    + n19)
187                + n1)
188                * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
189                    * (n216 * (n4 * x2 - n2) - n144 * (n4 * x1 - n2) + n192)
190                    - n24
191                        * (n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2))
192                        * (n27 * (n4 * x2 - n2).powi(2) - n36 * (n4 * x1 - n2) * (n4 * x2 - n2)
193                            + n48 * (n4 * x2 - n2)
194                            + n12 * (n4 * x1 - n2).powi(2)
195                            - n32 * (n4 * x1 - n2)
196                            + n18))))
197        / (n10.ln()
198            * ((n4 * x2 + n4 * x1 - n3).powi(2)
199                * (n3 * (n4 * x2 - n2).powi(2) + n6 * (n4 * x1 - n2) * (n4 * x2 - n2)
200                    - n14 * (n4 * x2 - n2)
201                    + n3 * (n4 * x1 - n2).powi(2)
202                    - n14 * (n4 * x1 - n2)
203                    + n19)
204                + n1)
205            * ((n2 * (n4 * x1 - n2) - n3 * (n4 * x2 - n2)).powi(2)
206                * (n27 * (n4 * x2 - n2).powi(2) - n36 * (n4 * x1 - n2) * (n4 * x2 - n2)
207                    + n48 * (n4 * x2 - n2)
208                    + n12 * (n4 * x1 - n2).powi(2)
209                    - n32 * (n4 * x1 - n2)
210                    + n18)
211                + n30));
212
213    [a, b]
214}
215
216/// Hessian of Picheny test function.
217pub fn picheny_hessian<T>(param: &[T; 2]) -> [[T; 2]; 2]
218where
219    T: Float + FromPrimitive,
220{
221    // oh dear...
222    let n2 = T::from_f64(2.0).unwrap();
223    let n3 = T::from_f64(3.0).unwrap();
224    let n4 = T::from_f64(4.0).unwrap();
225    let n5 = T::from_f64(5.0).unwrap();
226    let n6 = T::from_f64(6.0).unwrap();
227    let n10 = T::from_f64(10.0).unwrap();
228    let n11 = T::from_f64(11.0).unwrap();
229    let n96 = T::from_f64(96.0).unwrap();
230    let n144 = T::from_f64(144.0).unwrap();
231    let n192 = T::from_f64(192.0).unwrap();
232    let n277 = T::from_f64(277.0).unwrap();
233    let n432 = T::from_f64(432.0).unwrap();
234    let n576 = T::from_f64(576.0).unwrap();
235    let n768 = T::from_f64(768.0).unwrap();
236    let n864 = T::from_f64(864.0).unwrap();
237    let n896 = T::from_f64(896.0).unwrap();
238    let n1080 = T::from_f64(1080.0).unwrap();
239    let n1152 = T::from_f64(1152.0).unwrap();
240    let n1167 = T::from_f64(1167.0).unwrap();
241    let n1512 = T::from_f64(1512.0).unwrap();
242    let n1603 = T::from_f64(1603.0).unwrap();
243    let n2048 = T::from_f64(2048.0).unwrap();
244    let n2304 = T::from_f64(2304.0).unwrap();
245    let n2688 = T::from_f64(2688.0).unwrap();
246    let n3024 = T::from_f64(3024.0).unwrap();
247    let n5376 = T::from_f64(5376.0).unwrap();
248    let n5594 = T::from_f64(5594.0).unwrap();
249    let n5874 = T::from_f64(5874.0).unwrap();
250    let n6144 = T::from_f64(6144.0).unwrap();
251    let n6912 = T::from_f64(6912.0).unwrap();
252    let n9216 = T::from_f64(9216.0).unwrap();
253    let n13824 = T::from_f64(13824.0).unwrap();
254    let n20736 = T::from_f64(20736.0).unwrap();
255    let n21546 = T::from_f64(21546.0).unwrap();
256    let n24576 = T::from_f64(24576.0).unwrap();
257    let n27648 = T::from_f64(27648.0).unwrap();
258    let n31104 = T::from_f64(31104.0).unwrap();
259    let n36864 = T::from_f64(36864.0).unwrap();
260    let n77456 = T::from_f64(77456.0).unwrap();
261    let n82944 = T::from_f64(82944.0).unwrap();
262    let n85088 = T::from_f64(85088.0).unwrap();
263    let n98304 = T::from_f64(98304.0).unwrap();
264    let n118688 = T::from_f64(118688.0).unwrap();
265    let n124416 = T::from_f64(124416.0).unwrap();
266    let n127904 = T::from_f64(127904.0).unwrap();
267    let n155024 = T::from_f64(155024.0).unwrap();
268    let n163936 = T::from_f64(163936.0).unwrap();
269    let n165888 = T::from_f64(165888.0).unwrap();
270    let n248832 = T::from_f64(248832.0).unwrap();
271    let n255264 = T::from_f64(255264.0).unwrap();
272    let n255808 = T::from_f64(255808.0).unwrap();
273    let n327872 = T::from_f64(327872.0).unwrap();
274    let n331776 = T::from_f64(331776.0).unwrap();
275    let n393216 = T::from_f64(393216.0).unwrap();
276    let n442368 = T::from_f64(442368.0).unwrap();
277    let n497664 = T::from_f64(497664.0).unwrap();
278    let n516096 = T::from_f64(516096.0).unwrap();
279    let n679936 = T::from_f64(679936.0).unwrap();
280    let n688128 = T::from_f64(688128.0).unwrap();
281    let n705280 = T::from_f64(705280.0).unwrap();
282    let n995328 = T::from_f64(995328.0).unwrap();
283    let n1057536 = T::from_f64(1057536.0).unwrap();
284    let n1057920 = T::from_f64(1057920.0).unwrap();
285    let n1136640 = T::from_f64(1136640.0).unwrap();
286    let n1253376 = T::from_f64(1253376.0).unwrap();
287    let n1290240 = T::from_f64(1290240.0).unwrap();
288    let n1298688 = T::from_f64(1298688.0).unwrap();
289    let n1327104 = T::from_f64(1327104.0).unwrap();
290    let n1435136 = T::from_f64(1435136.0).unwrap();
291    let n1441280 = T::from_f64(1441280.0).unwrap();
292    let n1490944 = T::from_f64(1490944.0).unwrap();
293    let n1515520 = T::from_f64(1515520.0).unwrap();
294    let n1556736 = T::from_f64(1556736.0).unwrap();
295    let n1781760 = T::from_f64(1781760.0).unwrap();
296    let n1854464 = T::from_f64(1854464.0).unwrap();
297    let n1880064 = T::from_f64(1880064.0).unwrap();
298    let n1990656 = T::from_f64(1990656.0).unwrap();
299    let n2088960 = T::from_f64(2088960.0).unwrap();
300    let n2161920 = T::from_f64(2161920.0).unwrap();
301    let n2322432 = T::from_f64(2322432.0).unwrap();
302    let n2363520 = T::from_f64(2363520.0).unwrap();
303    let n2520320 = T::from_f64(2520320.0).unwrap();
304    let n2550112 = T::from_f64(2550112.0).unwrap();
305    let n2580480 = T::from_f64(2580480.0).unwrap();
306    let n2654208 = T::from_f64(2654208.0).unwrap();
307    let n2752512 = T::from_f64(2752512.0).unwrap();
308    let n2985984 = T::from_f64(2985984.0).unwrap();
309    let n3113472 = T::from_f64(3113472.0).unwrap();
310    let n3133440 = T::from_f64(3133440.0).unwrap();
311    let n3896064 = T::from_f64(3896064.0).unwrap();
312    let n4079616 = T::from_f64(4079616.0).unwrap();
313    let n4231680 = T::from_f64(4231680.0).unwrap();
314    let n4323840 = T::from_f64(4323840.0).unwrap();
315    let n4523520 = T::from_f64(4523520.0).unwrap();
316    let n4546560 = T::from_f64(4546560.0).unwrap();
317    let n5040640 = T::from_f64(5040640.0).unwrap();
318    let n5287680 = T::from_f64(5287680.0).unwrap();
319    let n5492736 = T::from_f64(5492736.0).unwrap();
320    let n5971968 = T::from_f64(5971968.0).unwrap();
321    let n6266880 = T::from_f64(6266880.0).unwrap();
322    let n6785280 = T::from_f64(6785280.0).unwrap();
323    let n6850560 = T::from_f64(6850560.0).unwrap();
324    let n7127040 = T::from_f64(7127040.0).unwrap();
325    let n7175680 = T::from_f64(7175680.0).unwrap();
326    let n7399680 = T::from_f64(7399680.0).unwrap();
327    let n7728640 = T::from_f64(7728640.0).unwrap();
328    let n8515584 = T::from_f64(8515584.0).unwrap();
329    let n9089280 = T::from_f64(9089280.0).unwrap();
330    let n9134080 = T::from_f64(9134080.0).unwrap();
331    let n9400320 = T::from_f64(9400320.0).unwrap();
332    let n9454080 = T::from_f64(9454080.0).unwrap();
333    let n10081280 = T::from_f64(10081280.0).unwrap();
334    let n13284864 = T::from_f64(13284864.0).unwrap();
335    let n13570560 = T::from_f64(13570560.0).unwrap();
336    let n13731840 = T::from_f64(13731840.0).unwrap();
337    let n13934592 = T::from_f64(13934592.0).unwrap();
338    let n14799360 = T::from_f64(14799360.0).unwrap();
339    let n23185920 = T::from_f64(23185920.0).unwrap();
340    let n27463680 = T::from_f64(27463680.0).unwrap();
341    let n29598720 = T::from_f64(29598720.0).unwrap();
342    let n27402240 = T::from_f64(27402240.0).unwrap();
343
344    let factor = T::from_f64(9.88875154511743).unwrap();
345
346    let [x_1, x_2] = *param;
347
348    let offdiag = -(factor
349        * (n768 * x_2.powi(3) + x_1 * (n2304 * x_2.powi(2) - n5376 * x_2 + n3024)
350            - n2688 * x_2.powi(2)
351            + x_1.powi(2) * (n2304 * x_2 - n2688)
352            + n3024 * x_2
353            + n768 * x_1.powi(3)
354            - n1080)
355        * (n331776 * x_2.powi(7) - n497664 * x_2.powi(6)
356            + x_1
357                * (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
358                    + n1441280 * x_2.powi(3)
359                    + n1556736 * x_2.powi(2)
360                    - n255808 * x_2
361                    + n5594)
362            - n1057536 * x_2.powi(5)
363            + x_1.powi(2)
364                * (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
365                    - n6785280 * x_2.powi(2)
366                    - n255264 * x_2
367                    + n77456)
368            + x_1.powi(3)
369                * (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
370                    + n5040640 * x_2.powi(2)
371                    + n4231680 * x_2
372                    - n118688)
373            + n2363520 * x_2.powi(4)
374            + x_1.powi(4)
375                * (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2) - n7175680 * x_2 - n705280)
376            - n1298688 * x_2.powi(3)
377            + n163936 * x_2.powi(2)
378            + x_1.powi(5) * (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464)
379            + n5874 * x_2
380            + x_1.powi(6) * (-n688128 * x_2 - n1490944)
381            + n393216 * x_1.powi(7)
382            - n1603))
383        / (n10.ln()
384            * (n192 * x_2.powi(4)
385                + x_1 * (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080)
386                - n896 * x_2.powi(3)
387                + x_1.powi(2) * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512)
388                + n1512 * x_2.powi(2)
389                + x_1.powi(3) * (n768 * x_2 - n896)
390                - n1080 * x_2
391                + n192 * x_1.powi(4)
392                + n277)
393                .powi(2)
394            * (n31104 * x_2.powi(4) - n6912 * x_2.powi(3)
395                + x_1 * (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96)
396                + x_1.powi(2) * (n82944 * x_2.powi(2) - n9216 * x_2 - n192)
397                - n432 * x_2.powi(2)
398                + n144 * x_2
399                + x_1.powi(3) * (n2048 - n36864 * x_2)
400                + n6144 * x_1.powi(4)
401                + n11))
402        - (factor
403            * (n124416 * x_2.powi(3) - n20736 * x_2.powi(2)
404                + x_1 * (-n248832 * x_2.powi(2) + n27648 * x_2 + n576)
405                + x_1.powi(2) * (n165888 * x_2 - n9216)
406                - n864 * x_2
407                - n36864 * x_1.powi(3)
408                + n144)
409            * (n331776 * x_2.powi(7) - n497664 * x_2.powi(6)
410                + x_1
411                    * (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
412                        + n1441280 * x_2.powi(3)
413                        + n1556736 * x_2.powi(2)
414                        - n255808 * x_2
415                        + n5594)
416                - n1057536 * x_2.powi(5)
417                + x_1.powi(2)
418                    * (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
419                        - n6785280 * x_2.powi(2)
420                        - n255264 * x_2
421                        + n77456)
422                + x_1.powi(3)
423                    * (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
424                        + n5040640 * x_2.powi(2)
425                        + n4231680 * x_2
426                        - n118688)
427                + n2363520 * x_2.powi(4)
428                + x_1.powi(4)
429                    * (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2)
430                        - n7175680 * x_2
431                        - n705280)
432                - n1298688 * x_2.powi(3)
433                + n163936 * x_2.powi(2)
434                + x_1.powi(5) * (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464)
435                + n5874 * x_2
436                + x_1.powi(6) * (-n688128 * x_2 - n1490944)
437                + n393216 * x_1.powi(7)
438                - n1603))
439            / (n10.ln()
440                * (n192 * x_2.powi(4)
441                    + x_1 * (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080)
442                    - n896 * x_2.powi(3)
443                    + x_1.powi(2) * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512)
444                    + n1512 * x_2.powi(2)
445                    + x_1.powi(3) * (n768 * x_2 - n896)
446                    - n1080 * x_2
447                    + n192 * x_1.powi(4)
448                    + n277)
449                * (n31104 * x_2.powi(4) - n6912 * x_2.powi(3)
450                    + x_1 * (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96)
451                    + x_1.powi(2) * (n82944 * x_2.powi(2) - n9216 * x_2 - n192)
452                    - n432 * x_2.powi(2)
453                    + n144 * x_2
454                    + x_1.powi(3) * (n2048 - n36864 * x_2)
455                    + n6144 * x_1.powi(4)
456                    + n11)
457                    .powi(2))
458        + (factor
459            * (n2322432 * x_2.powi(6) - n2985984 * x_2.powi(5)
460                + x_1
461                    * (-n5971968 * x_2.powi(5) + n27463680 * x_2.powi(4)
462                        - n29598720 * x_2.powi(3)
463                        + n4323840 * x_2.powi(2)
464                        + n3113472 * x_2
465                        - n255808)
466                - n5287680 * x_2.powi(4)
467                + x_1.powi(2)
468                    * (-n9400320 * x_2.powi(4)
469                        + n4546560 * x_2.powi(3)
470                        + n23185920 * x_2.powi(2)
471                        - n13570560 * x_2
472                        - n255264)
473                + x_1.powi(3)
474                    * (n7127040 * x_2.powi(3) - n27402240 * x_2.powi(2)
475                        + n10081280 * x_2
476                        + n4231680)
477                + n9454080 * x_2.powi(3)
478                + x_1.powi(4) * (n6266880 * x_2.powi(2) + n2580480 * x_2 - n7175680)
479                - n3896064 * x_2.powi(2)
480                + n327872 * x_2
481                + x_1.powi(5) * (n4079616 - n2654208 * x_2)
482                - n688128 * x_1.powi(6)
483                + n5874))
484            / (n10.ln()
485                * (n192 * x_2.powi(4)
486                    + x_1 * (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080)
487                    - n896 * x_2.powi(3)
488                    + x_1.powi(2) * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512)
489                    + n1512 * x_2.powi(2)
490                    + x_1.powi(3) * (n768 * x_2 - n896)
491                    - n1080 * x_2
492                    + n192 * x_1.powi(4)
493                    + n277)
494                * (n31104 * x_2.powi(4) - n6912 * x_2.powi(3)
495                    + x_1 * (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96)
496                    + x_1.powi(2) * (n82944 * x_2.powi(2) - n9216 * x_2 - n192)
497                    - n432 * x_2.powi(2)
498                    + n144 * x_2
499                    + x_1.powi(3) * (n2048 - n36864 * x_2)
500                    + n6144 * x_1.powi(4)
501                    + n11));
502
503    let b = -(factor
504        * (n768 * x_2.powi(3)
505            + n3 * (n768 * x_1 - n896) * x_2.powi(2)
506            + n2 * (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2
507            + n768 * x_1.powi(3)
508            - n2688 * x_1.powi(2)
509            + n3024 * x_1
510            - n1080)
511        * (n1990656 * x_2.powi(7)
512            + (n2322432 * x_1 - n8515584) * x_2.powi(6)
513            + (-n2985984 * x_1.powi(2) - n2985984 * x_1 + n13284864) * x_2.powi(5)
514            + (-n3133440 * x_1.powi(3) + n13731840 * x_1.powi(2) - n5287680 * x_1 - n9089280)
515                * x_2.powi(4)
516            + (n1781760 * x_1.powi(4) + n1515520 * x_1.powi(3) - n14799360 * x_1.powi(2)
517                + n9454080 * x_1
518                + n2550112)
519                * x_2.powi(3)
520            + (n1253376 * x_1.powi(5) - n6850560 * x_1.powi(4)
521                + n7728640 * x_1.powi(3)
522                + n2161920 * x_1.powi(2)
523                - n3896064 * x_1
524                - n155024)
525                * x_2.powi(2)
526            + (-n442368 * x_1.powi(6) + n516096 * x_1.powi(5) + n2520320 * x_1.powi(4)
527                - n4523520 * x_1.powi(3)
528                + n1556736 * x_1.powi(2)
529                + n327872 * x_1
530                - n21546)
531                * x_2
532            - n98304 * x_1.powi(7)
533            + n679936 * x_1.powi(6)
534            - n1435136 * x_1.powi(5)
535            + n1057920 * x_1.powi(4)
536            - n85088 * x_1.powi(3)
537            - n127904 * x_1.powi(2)
538            + n5874 * x_1
539            + n1167))
540        / (n10.ln()
541            * (n192 * x_2.powi(4)
542                + (n768 * x_1 - n896) * x_2.powi(3)
543                + (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2.powi(2)
544                + (n768 * x_1.powi(3) - n2688 * x_1.powi(2) + n3024 * x_1 - n1080) * x_2
545                + n192 * x_1.powi(4)
546                - n896 * x_1.powi(3)
547                + n1512 * x_1.powi(2)
548                - n1080 * x_1
549                + n277)
550                .powi(2)
551            * (n31104 * x_2.powi(4)
552                + (-n82944 * x_1 - n6912) * x_2.powi(3)
553                + (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2.powi(2)
554                + (-n36864 * x_1.powi(3) - n9216 * x_1.powi(2) + n576 * x_1 + n144) * x_2
555                + n6144 * x_1.powi(4)
556                + n2048 * x_1.powi(3)
557                - n192 * x_1.powi(2)
558                - n96 * x_1
559                + n11))
560        - (factor
561            * (n124416 * x_2.powi(3)
562                + n3 * (-n82944 * x_1 - n6912) * x_2.powi(2)
563                + n2 * (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2
564                - n36864 * x_1.powi(3)
565                - n9216 * x_1.powi(2)
566                + n576 * x_1
567                + n144)
568            * (n1990656 * x_2.powi(7)
569                + (n2322432 * x_1 - n8515584) * x_2.powi(6)
570                + (-n2985984 * x_1.powi(2) - n2985984 * x_1 + n13284864) * x_2.powi(5)
571                + (-n3133440 * x_1.powi(3) + n13731840 * x_1.powi(2) - n5287680 * x_1 - n9089280)
572                    * x_2.powi(4)
573                + (n1781760 * x_1.powi(4) + n1515520 * x_1.powi(3) - n14799360 * x_1.powi(2)
574                    + n9454080 * x_1
575                    + n2550112)
576                    * x_2.powi(3)
577                + (n1253376 * x_1.powi(5) - n6850560 * x_1.powi(4)
578                    + n7728640 * x_1.powi(3)
579                    + n2161920 * x_1.powi(2)
580                    - n3896064 * x_1
581                    - n155024)
582                    * x_2.powi(2)
583                + (-n442368 * x_1.powi(6) + n516096 * x_1.powi(5) + n2520320 * x_1.powi(4)
584                    - n4523520 * x_1.powi(3)
585                    + n1556736 * x_1.powi(2)
586                    + n327872 * x_1
587                    - n21546)
588                    * x_2
589                - n98304 * x_1.powi(7)
590                + n679936 * x_1.powi(6)
591                - n1435136 * x_1.powi(5)
592                + n1057920 * x_1.powi(4)
593                - n85088 * x_1.powi(3)
594                - n127904 * x_1.powi(2)
595                + n5874 * x_1
596                + n1167))
597            / (n10.ln()
598                * (n192 * x_2.powi(4)
599                    + (n768 * x_1 - n896) * x_2.powi(3)
600                    + (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2.powi(2)
601                    + (n768 * x_1.powi(3) - n2688 * x_1.powi(2) + n3024 * x_1 - n1080) * x_2
602                    + n192 * x_1.powi(4)
603                    - n896 * x_1.powi(3)
604                    + n1512 * x_1.powi(2)
605                    - n1080 * x_1
606                    + n277)
607                * (n31104 * x_2.powi(4)
608                    + (-n82944 * x_1 - n6912) * x_2.powi(3)
609                    + (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2.powi(2)
610                    + (-n36864 * x_1.powi(3) - n9216 * x_1.powi(2) + n576 * x_1 + n144) * x_2
611                    + n6144 * x_1.powi(4)
612                    + n2048 * x_1.powi(3)
613                    - n192 * x_1.powi(2)
614                    - n96 * x_1
615                    + n11)
616                    .powi(2))
617        + (factor
618            * (n13934592 * x_2.powi(6)
619                + n6 * (n2322432 * x_1 - n8515584) * x_2.powi(5)
620                + n5 * (-n2985984 * x_1.powi(2) - n2985984 * x_1 + n13284864) * x_2.powi(4)
621                + n4 * (-n3133440 * x_1.powi(3) + n13731840 * x_1.powi(2)
622                    - n5287680 * x_1
623                    - n9089280)
624                    * x_2.powi(3)
625                + n3 * (n1781760 * x_1.powi(4) + n1515520 * x_1.powi(3)
626                    - n14799360 * x_1.powi(2)
627                    + n9454080 * x_1
628                    + n2550112)
629                    * x_2.powi(2)
630                + n2 * (n1253376 * x_1.powi(5) - n6850560 * x_1.powi(4)
631                    + n7728640 * x_1.powi(3)
632                    + n2161920 * x_1.powi(2)
633                    - n3896064 * x_1
634                    - n155024)
635                    * x_2
636                - n442368 * x_1.powi(6)
637                + n516096 * x_1.powi(5)
638                + n2520320 * x_1.powi(4)
639                - n4523520 * x_1.powi(3)
640                + n1556736 * x_1.powi(2)
641                + n327872 * x_1
642                - n21546))
643            / (n10.ln()
644                * (n192 * x_2.powi(4)
645                    + (n768 * x_1 - n896) * x_2.powi(3)
646                    + (n1152 * x_1.powi(2) - n2688 * x_1 + n1512) * x_2.powi(2)
647                    + (n768 * x_1.powi(3) - n2688 * x_1.powi(2) + n3024 * x_1 - n1080) * x_2
648                    + n192 * x_1.powi(4)
649                    - n896 * x_1.powi(3)
650                    + n1512 * x_1.powi(2)
651                    - n1080 * x_1
652                    + n277)
653                * (n31104 * x_2.powi(4)
654                    + (-n82944 * x_1 - n6912) * x_2.powi(3)
655                    + (n82944 * x_1.powi(2) + n13824 * x_1 - n432) * x_2.powi(2)
656                    + (-n36864 * x_1.powi(3) - n9216 * x_1.powi(2) + n576 * x_1 + n144) * x_2
657                    + n6144 * x_1.powi(4)
658                    + n2048 * x_1.powi(3)
659                    - n192 * x_1.powi(2)
660                    - n96 * x_1
661                    + n11));
662
663    let a = -(factor
664        * (n768 * x_1.powi(3)
665            + n3 * (n768 * x_2 - n896) * x_1.powi(2)
666            + n2 * (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1
667            + n768 * x_2.powi(3)
668            - n2688 * x_2.powi(2)
669            + n3024 * x_2
670            - n1080)
671        * (n393216 * x_1.powi(7)
672            + (-n688128 * x_2 - n1490944) * x_1.powi(6)
673            + (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464) * x_1.powi(5)
674            + (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2) - n7175680 * x_2 - n705280)
675                * x_1.powi(4)
676            + (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
677                + n5040640 * x_2.powi(2)
678                + n4231680 * x_2
679                - n118688)
680                * x_1.powi(3)
681            + (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
682                - n6785280 * x_2.powi(2)
683                - n255264 * x_2
684                + n77456)
685                * x_1.powi(2)
686            + (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
687                + n1441280 * x_2.powi(3)
688                + n1556736 * x_2.powi(2)
689                - n255808 * x_2
690                + n5594)
691                * x_1
692            + n331776 * x_2.powi(7)
693            - n497664 * x_2.powi(6)
694            - n1057536 * x_2.powi(5)
695            + n2363520 * x_2.powi(4)
696            - n1298688 * x_2.powi(3)
697            + n163936 * x_2.powi(2)
698            + n5874 * x_2
699            - n1603))
700        / (n10.ln()
701            * (n192 * x_1.powi(4)
702                + (n768 * x_2 - n896) * x_1.powi(3)
703                + (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1.powi(2)
704                + (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080) * x_1
705                + n192 * x_2.powi(4)
706                - n896 * x_2.powi(3)
707                + n1512 * x_2.powi(2)
708                - n1080 * x_2
709                + n277)
710                .powi(2)
711            * (n6144 * x_1.powi(4)
712                + (n2048 - n36864 * x_2) * x_1.powi(3)
713                + (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1.powi(2)
714                + (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96) * x_1
715                + n31104 * x_2.powi(4)
716                - n6912 * x_2.powi(3)
717                - n432 * x_2.powi(2)
718                + n144 * x_2
719                + n11))
720        - (factor
721            * (n24576 * x_1.powi(3)
722                + n3 * (n2048 - n36864 * x_2) * x_1.powi(2)
723                + n2 * (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1
724                - n82944 * x_2.powi(3)
725                + n13824 * x_2.powi(2)
726                + n576 * x_2
727                - n96)
728            * (n393216 * x_1.powi(7)
729                + (-n688128 * x_2 - n1490944) * x_1.powi(6)
730                + (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464) * x_1.powi(5)
731                + (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2) - n7175680 * x_2 - n705280)
732                    * x_1.powi(4)
733                + (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
734                    + n5040640 * x_2.powi(2)
735                    + n4231680 * x_2
736                    - n118688)
737                    * x_1.powi(3)
738                + (-n1880064 * x_2.powi(5) + n1136640 * x_2.powi(4) + n7728640 * x_2.powi(3)
739                    - n6785280 * x_2.powi(2)
740                    - n255264 * x_2
741                    + n77456)
742                    * x_1.powi(2)
743                + (-n995328 * x_2.powi(6) + n5492736 * x_2.powi(5) - n7399680 * x_2.powi(4)
744                    + n1441280 * x_2.powi(3)
745                    + n1556736 * x_2.powi(2)
746                    - n255808 * x_2
747                    + n5594)
748                    * x_1
749                + n331776 * x_2.powi(7)
750                - n497664 * x_2.powi(6)
751                - n1057536 * x_2.powi(5)
752                + n2363520 * x_2.powi(4)
753                - n1298688 * x_2.powi(3)
754                + n163936 * x_2.powi(2)
755                + n5874 * x_2
756                - n1603))
757            / (n10.ln()
758                * (n192 * x_1.powi(4)
759                    + (n768 * x_2 - n896) * x_1.powi(3)
760                    + (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1.powi(2)
761                    + (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080) * x_1
762                    + n192 * x_2.powi(4)
763                    - n896 * x_2.powi(3)
764                    + n1512 * x_2.powi(2)
765                    - n1080 * x_2
766                    + n277)
767                * (n6144 * x_1.powi(4)
768                    + (n2048 - n36864 * x_2) * x_1.powi(3)
769                    + (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1.powi(2)
770                    + (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96) * x_1
771                    + n31104 * x_2.powi(4)
772                    - n6912 * x_2.powi(3)
773                    - n432 * x_2.powi(2)
774                    + n144 * x_2
775                    + n11)
776                    .powi(2))
777        + (factor
778            * (n2752512 * x_1.powi(6)
779                + n6 * (-n688128 * x_2 - n1490944) * x_1.powi(5)
780                + n5 * (-n1327104 * x_2.powi(2) + n4079616 * x_2 + n1854464) * x_1.powi(4)
781                + n4 * (n2088960 * x_2.powi(3) + n1290240 * x_2.powi(2)
782                    - n7175680 * x_2
783                    - n705280)
784                    * x_1.powi(3)
785                + n3 * (n1781760 * x_2.powi(4) - n9134080 * x_2.powi(3)
786                    + n5040640 * x_2.powi(2)
787                    + n4231680 * x_2
788                    - n118688)
789                    * x_1.powi(2)
790                + n2 * (-n1880064 * x_2.powi(5)
791                    + n1136640 * x_2.powi(4)
792                    + n7728640 * x_2.powi(3)
793                    - n6785280 * x_2.powi(2)
794                    - n255264 * x_2
795                    + n77456)
796                    * x_1
797                - n995328 * x_2.powi(6)
798                + n5492736 * x_2.powi(5)
799                - n7399680 * x_2.powi(4)
800                + n1441280 * x_2.powi(3)
801                + n1556736 * x_2.powi(2)
802                - n255808 * x_2
803                + n5594))
804            / (n10.ln()
805                * (n192 * x_1.powi(4)
806                    + (n768 * x_2 - n896) * x_1.powi(3)
807                    + (n1152 * x_2.powi(2) - n2688 * x_2 + n1512) * x_1.powi(2)
808                    + (n768 * x_2.powi(3) - n2688 * x_2.powi(2) + n3024 * x_2 - n1080) * x_1
809                    + n192 * x_2.powi(4)
810                    - n896 * x_2.powi(3)
811                    + n1512 * x_2.powi(2)
812                    - n1080 * x_2
813                    + n277)
814                * (n6144 * x_1.powi(4)
815                    + (n2048 - n36864 * x_2) * x_1.powi(3)
816                    + (n82944 * x_2.powi(2) - n9216 * x_2 - n192) * x_1.powi(2)
817                    + (-n82944 * x_2.powi(3) + n13824 * x_2.powi(2) + n576 * x_2 - n96) * x_1
818                    + n31104 * x_2.powi(4)
819                    - n6912 * x_2.powi(3)
820                    - n432 * x_2.powi(2)
821                    + n144 * x_2
822                    + n11));
823
824    [[a, offdiag], [offdiag, b]]
825}
826
827#[cfg(test)]
828mod tests {
829    use super::*;
830    use approx::assert_relative_eq;
831    use finitediff::FiniteDiff;
832    use proptest::prelude::*;
833    use std::{f32, f64};
834
835    #[test]
836    fn test_picheny_optimum() {
837        assert_relative_eq!(
838            picheny(&[0.5_f32, 0.25_f32]),
839            -3.3851993182,
840            epsilon = f32::EPSILON
841        );
842        assert_relative_eq!(
843            picheny(&[0.5_f64, 0.25_f64]),
844            -3.3851993182036826,
845            epsilon = f64::EPSILON
846        );
847
848        let deriv = picheny_derivative(&[0.5_f64, 0.25_f64]);
849        // println!("1: {deriv:?}");
850        for i in 0..2 {
851            assert_relative_eq!(deriv[i], 0.0, epsilon = f64::EPSILON);
852        }
853    }
854
855    proptest! {
856        #[test]
857        fn test_picheny_derivative(a in 0.0..1.0, b in 0.0..1.0) {
858            let param = [a, b];
859            let derivative = picheny_derivative(&param);
860            let derivative_fd = Vec::from(param).central_diff(&|x| picheny(&[x[0], x[1]]));
861            // println!("1: {derivative:?} at {a}/{b}");
862            // println!("2: {derivative_fd:?} at {a}/{b}");
863            for i in 0..derivative.len() {
864                assert_relative_eq!(
865                    derivative[i],
866                    derivative_fd[i],
867                    epsilon = 1e-5,
868                    max_relative = 1e-2
869                );
870            }
871        }
872    }
873
874    proptest! {
875        #[test]
876        fn test_picheny_hessian_finitediff(a in 0.0..1.0, b in 0.0..1.0) {
877            let param = [a, b];
878            let hessian = picheny_hessian(&param);
879            let hessian_fd =
880                Vec::from(param).central_hessian(&|x| picheny_derivative(&[x[0], x[1]]).to_vec());
881            let n = hessian.len();
882            // println!("1: {hessian:?} at {a}/{b}");
883            // println!("2: {hessian_fd:?} at {a}/{b}");
884            for i in 0..n {
885                assert_eq!(hessian[i].len(), n);
886                for j in 0..n {
887                    if hessian_fd[i][j].is_finite() {
888                        assert_relative_eq!(
889                            hessian[i][j],
890                            hessian_fd[i][j],
891                            epsilon = 1e-5,
892                            max_relative = 1e-2
893                        );
894                    }
895                }
896            }
897        }
898    }
899}
argmin_testfunctions/picheny.rs

argmin_testfunctions/
picheny.rs
