Skip to content

Commit 128ec4c

Browse files
stefan-kstefan-k
committed
Added derivative and Hessian for Levy test function
1 parent 54273b0 commit 128ec4c

File tree

2 files changed

+272
-2
lines changed

2 files changed

+272
-2
lines changed

‎crates/argmin-testfunctions/proptest-regressions/levy.txt

+1
Original file line numberDiff line numberDiff line change
@@ -6,3 +6,4 @@
66
# everyone who runs the test benefits from these saved cases.
77
cc 0ee9ddf2f929b971ef215e55319945f24fe713dc7f41242b8e80074112a4647d # shrinks to a = 0.0, b = -340.45222947026997
88
cc a9994e5eb0e4c5161426f9fdec86fb46ee637f17a62a251c5e07c58123f960e9 # shrinks to a = 3.219661442452062, b = -2.2175779827682494
9+
cc 0dee901c5d94a9a9b32352c229a6c9d196aa5a27c27f0636bfbd5b66b7e95fc5 # shrinks to a = 0.0, b = 6.072194039937718, c = 0.0

‎crates/argmin-testfunctions/src/levy.rs

+271-2
Original file line numberDiff line numberDiff line change
@@ -54,17 +54,188 @@ where
5454
let n10 = T::from_f64(10.0).unwrap();
5555
let pi = T::from_f64(PI).unwrap();
5656

57-
let w = |x: T| n1 - (x - n1) / n4;
57+
let w = |x: T| n1 + (x - n1) / n4;
5858

5959
(pi * w(param[0])).sin().powi(2)
6060
+ param[1..(plen - 1)]
6161
.iter()
6262
.map(|x| w(*x))
63-
.map(|wi: T| (wi - n1).powi(2) * (n1 + n10 * (pi * wi + n1)))
63+
.map(|wi: T| (wi - n1).powi(2) * (n1 + n10 * (pi * wi + n1).sin().powi(2)))
6464
.sum()
6565
+ (w(param[plen - 1]) - n1).powi(2) * (n1 + (n2 * pi * w(param[plen - 1])).sin().powi(2))
6666
}
6767

68+
/// Derivative of Levy test function
69+
pub fn levy_derivative<T>(param: &[T]) -> Vec<T>
70+
where
71+
T: Float + FromPrimitive + Sum,
72+
{
73+
let d = param.len();
74+
assert!(d >= 2);
75+
76+
let n1 = T::from_f64(1.0).unwrap();
77+
let n2 = T::from_f64(2.0).unwrap();
78+
let n4 = T::from_f64(4.0).unwrap();
79+
let n5 = T::from_f64(5.0).unwrap();
80+
let n8 = T::from_f64(8.0).unwrap();
81+
let n10 = T::from_f64(10.0).unwrap();
82+
let n16 = T::from_f64(16.0).unwrap();
83+
let pi = T::from_f64(PI).unwrap();
84+
85+
param
86+
.iter()
87+
.enumerate()
88+
.map(|(i, x)| (i, x, pi * ((*x - n1) / n4 + n1)))
89+
.map(|(i, &x, wp)| {
90+
if i == 0 {
91+
pi / n2 * wp.cos() * wp.sin()
92+
} else if i == d - 1 {
93+
((n2 * wp).sin().powi(2) + n1) * (x - n1) / n8
94+
+ pi / n16 * (n2 * wp).cos() * (n2 * wp).sin() * (x - n1).powi(2)
95+
} else {
96+
(n10 * (wp + n1).sin().powi(2) + n1) * (x - n1) / n8
97+
+ n5 / n16 * pi * (wp + n1).cos() * (wp + n1).sin() * (x - n1).powi(2)
98+
}
99+
})
100+
.collect()
101+
}
102+
103+
/// Derivative of Levy test function
104+
///
105+
/// This is the const generics version, which requires the number of parameters to be known
106+
/// at compile time.
107+
pub fn levy_derivative_const<const N: usize, T>(param: &[T; N]) -> [T; N]
108+
where
109+
T: Float + FromPrimitive + Sum,
110+
{
111+
assert!(N >= 2);
112+
113+
let n1 = T::from_f64(1.0).unwrap();
114+
let n0 = T::from_f64(0.0).unwrap();
115+
let n2 = T::from_f64(2.0).unwrap();
116+
let n4 = T::from_f64(4.0).unwrap();
117+
let n5 = T::from_f64(5.0).unwrap();
118+
let n8 = T::from_f64(8.0).unwrap();
119+
let n10 = T::from_f64(10.0).unwrap();
120+
let n16 = T::from_f64(16.0).unwrap();
121+
let pi = T::from_f64(PI).unwrap();
122+
123+
let mut out = [n0; N];
124+
125+
param
126+
.iter()
127+
.zip(out.iter_mut())
128+
.enumerate()
129+
.map(|(i, (x, o))| (i, x, pi * ((*x - n1) / n4 + n1), o))
130+
.map(|(i, &x, wp, o)| {
131+
*o = if i == 0 {
132+
pi / n2 * wp.cos() * wp.sin()
133+
} else if i == N - 1 {
134+
((n2 * wp).sin().powi(2) + n1) * (x - n1) / n8
135+
+ pi / n16 * (n2 * wp).cos() * (n2 * wp).sin() * (x - n1).powi(2)
136+
} else {
137+
(n10 * (wp + n1).sin().powi(2) + n1) * (x - n1) / n8
138+
+ n5 / n16 * pi * (wp + n1).cos() * (wp + n1).sin() * (x - n1).powi(2)
139+
}
140+
})
141+
.count();
142+
out
143+
}
144+
145+
/// Hessian of Levy test function
146+
pub fn levy_hessian<T>(param: &[T]) -> Vec<Vec<T>>
147+
where
148+
T: Float + FromPrimitive + Sum,
149+
{
150+
let d = param.len();
151+
assert!(d >= 2);
152+
153+
let x = param;
154+
155+
let n0 = T::from_f64(0.0).unwrap();
156+
let n1 = T::from_f64(1.0).unwrap();
157+
let n2 = T::from_f64(2.0).unwrap();
158+
let n4 = T::from_f64(4.0).unwrap();
159+
let n5 = T::from_f64(5.0).unwrap();
160+
let n6 = T::from_f64(6.0).unwrap();
161+
let n8 = T::from_f64(8.0).unwrap();
162+
let n10 = T::from_f64(10.0).unwrap();
163+
let n32 = T::from_f64(32.0).unwrap();
164+
let n64 = T::from_f64(64.0).unwrap();
165+
let pi = T::from_f64(PI).unwrap();
166+
let pi2 = T::from_f64(PI.powi(2)).unwrap();
167+
168+
let mut out = vec![vec![n0; d]; d];
169+
170+
for i in 0..d {
171+
let xin1 = x[i] - n1;
172+
let wp = pi * (xin1 / n4 + n1);
173+
out[i][i] = if i == 0 {
174+
pi2 / n8 * (wp.cos().powi(2) - wp.sin().powi(2))
175+
} else if i == d - 1 {
176+
-(n4 * pi * xin1 * (pi * x[i]).sin() + (pi2 * xin1.powi(2) - n2) * (pi * x[i]).cos()
177+
- n6)
178+
/ n32
179+
} else {
180+
let wp1cos = (wp + n1).cos();
181+
let wp1sin = (wp + n1).sin();
182+
n5 / n4 * pi * wp1cos * wp1sin * xin1
183+
+ n5 / n64 * pi2 * xin1.powi(2) * (wp1cos.powi(2) - wp1sin.powi(2))
184+
+ (n10 * wp1sin.powi(2) + n1) / n8
185+
}
186+
}
187+
188+
out
189+
}
190+
191+
/// Hessian of Levy test function
192+
///
193+
/// This is the const generics version, which requires the number of parameters to be known
194+
/// at compile time.
195+
pub fn levy_hessian_const<const N: usize, T>(param: &[T; N]) -> [[T; N]; N]
196+
where
197+
T: Float + FromPrimitive + Sum,
198+
{
199+
assert!(N >= 2);
200+
201+
let x = param;
202+
203+
let n0 = T::from_f64(0.0).unwrap();
204+
let n1 = T::from_f64(1.0).unwrap();
205+
let n2 = T::from_f64(2.0).unwrap();
206+
let n4 = T::from_f64(4.0).unwrap();
207+
let n5 = T::from_f64(5.0).unwrap();
208+
let n6 = T::from_f64(6.0).unwrap();
209+
let n8 = T::from_f64(8.0).unwrap();
210+
let n10 = T::from_f64(10.0).unwrap();
211+
let n32 = T::from_f64(32.0).unwrap();
212+
let n64 = T::from_f64(64.0).unwrap();
213+
let pi = T::from_f64(PI).unwrap();
214+
let pi2 = T::from_f64(PI.powi(2)).unwrap();
215+
216+
let mut out = [[n0; N]; N];
217+
218+
for i in 0..N {
219+
let xin1 = x[i] - n1;
220+
let wp = pi * (xin1 / n4 + n1);
221+
out[i][i] = if i == 0 {
222+
pi2 / n8 * (wp.cos().powi(2) - wp.sin().powi(2))
223+
} else if i == N - 1 {
224+
-(n4 * pi * xin1 * (pi * x[i]).sin() + (pi2 * xin1.powi(2) - n2) * (pi * x[i]).cos()
225+
- n6)
226+
/ n32
227+
} else {
228+
let wp1cos = (wp + n1).cos();
229+
let wp1sin = (wp + n1).sin();
230+
n5 / n4 * pi * wp1cos * wp1sin * xin1
231+
+ n5 / n64 * pi2 * xin1.powi(2) * (wp1cos.powi(2) - wp1sin.powi(2))
232+
+ (n10 * wp1sin.powi(2) + n1) / n8
233+
}
234+
}
235+
236+
out
237+
}
238+
68239
/// Levy test function No. 13
69240
///
70241
/// Defined as
@@ -180,6 +351,16 @@ mod tests {
180351
fn test_levy_optimum() {
181352
assert_relative_eq!(levy(&[1_f32, 1_f32, 1_f32]), 0.0, epsilon = f32::EPSILON);
182353
assert_relative_eq!(levy(&[1_f64, 1_f64, 1_f64]), 0.0, epsilon = f64::EPSILON);
354+
355+
let deriv = levy_derivative(&[1_f64, 1_f64, 1_f64]);
356+
for i in 0..2 {
357+
assert_relative_eq!(deriv[i], 0.0, epsilon = 1e-12, max_relative = 1e-12);
358+
}
359+
360+
let deriv = levy_derivative_const(&[1_f64, 1_f64, 1_f64]);
361+
for i in 0..2 {
362+
assert_relative_eq!(deriv[i], 0.0, epsilon = 1e-12, max_relative = 1e-12);
363+
}
183364
}
184365

185366
#[test]
@@ -243,4 +424,92 @@ mod tests {
243424
}
244425
}
245426
}
427+
428+
proptest! {
429+
#[test]
430+
fn test_levy_derivative_finitediff(a in -10.0..10.0, b in -10.0..10.0, c in -10.0..10.0) {
431+
let param = [a, b, c];
432+
let derivative = levy_derivative(&param);
433+
let derivative_fd = Vec::from(param).central_diff(&|x| levy(&[x[0], x[1], x[2]]));
434+
// println!("1: {derivative:?} at {param:?}");
435+
// println!("2: {derivative_fd:?} at {param:?}");
436+
for i in 0..derivative.len() {
437+
assert_relative_eq!(
438+
derivative[i],
439+
derivative_fd[i],
440+
epsilon = 1e-6,
441+
max_relative = 1e-4,
442+
);
443+
}
444+
}
445+
}
446+
447+
proptest! {
448+
#[test]
449+
fn test_levy_derivative_const_finitediff(a in -10.0..10.0, b in -10.0..10.0, c in -10.0..10.0) {
450+
let param = [a, b, c];
451+
let derivative = levy_derivative_const(&param);
452+
let derivative_fd = Vec::from(param).central_diff(&|x| levy(&[x[0], x[1], x[2]]));
453+
// println!("1: {derivative:?} at {param:?}");
454+
// println!("2: {derivative_fd:?} at {param:?}");
455+
for i in 0..derivative.len() {
456+
assert_relative_eq!(
457+
derivative[i],
458+
derivative_fd[i],
459+
epsilon = 1e-6,
460+
max_relative = 1e-4,
461+
);
462+
}
463+
}
464+
}
465+
466+
proptest! {
467+
#[test]
468+
fn test_levy_hessian_finitediff(a in -10.0..10.0, b in -10.0..10.0, c in -10.0..10.0) {
469+
let param = [a, b, c];
470+
let hessian = levy_hessian(&param);
471+
let hessian_fd = Vec::from(param).central_hessian(&|x| levy_derivative(&x).to_vec());
472+
let n = hessian.len();
473+
// println!("1: {hessian:?} at {param:?}");
474+
// println!("2: {hessian_fd:?} at {param:?}");
475+
for i in 0..n {
476+
assert_eq!(hessian[i].len(), n);
477+
for j in 0..n {
478+
if hessian[i][j].is_finite() {
479+
assert_relative_eq!(
480+
hessian[i][j],
481+
hessian_fd[i][j],
482+
epsilon = f64::EPSILON,
483+
max_relative = 1e-3
484+
);
485+
}
486+
}
487+
}
488+
}
489+
}
490+
491+
proptest! {
492+
#[test]
493+
fn test_levy_hessian_const_finitediff(a in -10.0..10.0, b in -10.0..10.0, c in -10.0..10.0) {
494+
let param = [a, b, c];
495+
let hessian = levy_hessian_const(&param);
496+
let hessian_fd = Vec::from(param).central_hessian(&|x| levy_derivative(&x).to_vec());
497+
let n = hessian.len();
498+
// println!("1: {hessian:?} at {param:?}");
499+
// println!("2: {hessian_fd:?} at {param:?}");
500+
for i in 0..n {
501+
assert_eq!(hessian[i].len(), n);
502+
for j in 0..n {
503+
if hessian[i][j].is_finite() {
504+
assert_relative_eq!(
505+
hessian[i][j],
506+
hessian_fd[i][j],
507+
epsilon = f64::EPSILON,
508+
max_relative = 1e-3
509+
);
510+
}
511+
}
512+
}
513+
}
514+
}
246515
}

0 commit comments

Comments
 (0)