1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
use super::Constraint;
#[derive(Copy, Clone)]
/// A simplex with level $\alpha$ is a set of the form
/// $\Delta_\alpha^n = \\{x \in \mathbb{R}^n {}:{} x \geq 0, \sum_i x_i = \alpha\\}$,
/// where $\alpha$ is a positive constant.
pub struct Simplex {
/// Simplex level
alpha: f64,
}
impl Simplex {
/// Construct a new simplex with given (positive) $\alpha$. The user does not need
/// to specify the dimension of the simplex.
pub fn new(alpha: f64) -> Self {
assert!(alpha > 0.0, "alpha is nonpositive");
Simplex { alpha }
}
}
impl Constraint for Simplex {
/// Project onto $\Delta_\alpha^n$ using Condat's fast projection algorithm.
///
/// See: Laurent Condat. Fast Projection onto the Simplex and the $\ell_1$ Ball.
/// <em>Mathematical Programming, Series A,</em> Springer, 2016, 158 (1), pp.575-585.
/// ⟨<a href="https://dx.doi.org/10.1007/s10107-015-0946-6">10.1007/s10107-015-0946-6</a>⟩.
fn project(&self, x: &mut [f64]) {
let a = &self.alpha;
// ---- step 1
let mut v = Vec::<f64>::with_capacity(x.len()); // vector containing x[0]
v.push(x[0]);
let mut v_size_old: i64 = -1; // 64 bit signed int
let mut v_tilde: Vec<f64> = Vec::new(); // empty vector of f64
let mut rho: f64 = x[0] - a; // 64 bit float
// ---- step 2
x.iter().skip(1).for_each(|x_n| {
if *x_n > rho {
rho += (*x_n - rho) / ((v.len() + 1) as f64);
if rho > *x_n - a {
v.push(*x_n);
} else {
v_tilde.extend(&v);
v = vec![*x_n];
rho = *x_n - a;
}
}
});
// ---- step 3
if !v_tilde.is_empty() {
v_tilde.iter().for_each(|v_t_n| {
if *v_t_n > rho {
v.push(*v_t_n);
rho += (*v_t_n - rho) / (v.len() as f64);
}
});
}
// ---- step 4
let mut keep_running = true;
while keep_running {
let mut hit_list: Vec<usize> = Vec::with_capacity(x.len());
let mut current_len_v = v.len() as i64;
v.iter().enumerate().for_each(|(n, v_n)| {
if *v_n <= rho {
hit_list.push(n);
current_len_v -= 1;
rho += (rho - *v_n) / (current_len_v as f64);
}
});
hit_list.iter().rev().for_each(|target| {
// remove in reverse to keep indexing correct
v.remove(*target);
});
keep_running = current_len_v != v_size_old;
v_size_old = current_len_v;
}
// ---- step 6
let zero: f64 = 0.0;
x.iter_mut().for_each(|x_n| *x_n = zero.max(*x_n - rho));
}
fn is_convex(&self) -> bool {
true
}
}