sciforge_lib/maths/optimization/
constrained.rs

1pub fn penalty_method(
2    f: fn(&[f64]) -> f64,
3    constraints: &[fn(&[f64]) -> f64],
4    x: &[f64],
5    penalty: f64,
6) -> f64 {
7    let mut cost = f(x);
8    for c in constraints {
9        let violation = c(x).max(0.0);
10        cost += penalty * violation * violation;
11    }
12    cost
13}
14
15pub fn augmented_lagrangian(
16    f: fn(&[f64]) -> f64,
17    constraints: &[fn(&[f64]) -> f64],
18    x: &[f64],
19    lambdas: &[f64],
20    mu: f64,
21) -> f64 {
22    let mut cost = f(x);
23    for (i, c) in constraints.iter().enumerate() {
24        let ci = c(x);
25        let li = if i < lambdas.len() { lambdas[i] } else { 0.0 };
26        cost += li * ci + mu / 2.0 * ci * ci;
27    }
28    cost
29}
30
31pub fn barrier_method(
32    f: fn(&[f64]) -> f64,
33    inequalities: &[fn(&[f64]) -> f64],
34    x: &[f64],
35    t: f64,
36) -> f64 {
37    let mut cost = t * f(x);
38    for g in inequalities {
39        let gi = g(x);
40        if gi >= 0.0 {
41            return f64::INFINITY;
42        }
43        cost -= (-gi).ln();
44    }
45    cost
46}
47
48pub fn project_box(x: &[f64], lower: &[f64], upper: &[f64]) -> Vec<f64> {
49    x.iter()
50        .zip(lower.iter().zip(upper.iter()))
51        .map(|(&xi, (&lo, &hi))| xi.clamp(lo, hi))
52        .collect()
53}
54
55pub fn project_simplex(x: &[f64]) -> Vec<f64> {
56    let mut sorted: Vec<f64> = x.to_vec();
57    sorted.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));
58    let mut cumsum = 0.0;
59    let mut rho = 0;
60    for (i, &si) in sorted.iter().enumerate() {
61        cumsum += si;
62        if si - (cumsum - 1.0) / (i + 1) as f64 > 0.0 {
63            rho = i + 1;
64        }
65    }
66    let theta = (sorted[..rho].iter().sum::<f64>() - 1.0) / rho as f64;
67    x.iter().map(|&xi| (xi - theta).max(0.0)).collect()
68}
69
70pub fn kkt_violation(
71    grad_f: &[f64],
72    constraints: &[f64],
73    lambdas: &[f64],
74    grad_constraints: &[Vec<f64>],
75) -> f64 {
76    let n = grad_f.len();
77    let mut stationarity = vec![0.0; n];
78    for (i, si) in stationarity.iter_mut().enumerate() {
79        *si = grad_f[i];
80        for (j, gc) in grad_constraints.iter().enumerate() {
81            if j < lambdas.len() {
82                *si += lambdas[j] * gc[i];
83            }
84        }
85    }
86    let stat_norm: f64 = stationarity.iter().map(|x| x * x).sum::<f64>().sqrt();
87    let feas: f64 = constraints
88        .iter()
89        .map(|&c| c.max(0.0).powi(2))
90        .sum::<f64>()
91        .sqrt();
92    let comp: f64 = lambdas
93        .iter()
94        .zip(constraints.iter())
95        .map(|(&l, &c)| (l * c).abs())
96        .sum::<f64>();
97    stat_norm + feas + comp
98}
99
100pub fn lagrangian(f: f64, constraints: &[f64], lambdas: &[f64]) -> f64 {
101    let mut l = f;
102    for (i, &c) in constraints.iter().enumerate() {
103        if i < lambdas.len() {
104            l += lambdas[i] * c;
105        }
106    }
107    l
108}
109
110pub fn projected_gradient_step(
111    grad: &[f64],
112    x: &[f64],
113    lr: f64,
114    lower: &[f64],
115    upper: &[f64],
116) -> Vec<f64> {
117    x.iter()
118        .enumerate()
119        .map(|(i, &xi)| (xi - lr * grad[i]).clamp(lower[i], upper[i]))
120        .collect()
121}
122
123pub fn frank_wolfe_step(grad: &[f64], lower: &[f64], upper: &[f64]) -> Vec<f64> {
124    grad.iter()
125        .enumerate()
126        .map(|(i, &gi)| if gi > 0.0 { lower[i] } else { upper[i] })
127        .collect()
128}
129
130pub fn admm_x_update(a: &[Vec<f64>], b: &[f64], z: &[f64], u: &[f64], rho: f64) -> Vec<f64> {
131    let n = b.len();
132    let mut x = vec![0.0; n];
133    for (i, xi) in x.iter_mut().enumerate() {
134        *xi = (b[i] + rho * (z[i] - u[i])) / (a[i][i] + rho);
135    }
136    x
137}
138
139pub fn dual_ascent_step(lambdas: &[f64], constraints: &[f64], step_size: f64) -> Vec<f64> {
140    lambdas
141        .iter()
142        .zip(constraints.iter())
143        .map(|(&l, &c)| (l + step_size * c).max(0.0))
144        .collect()
145}
146
147pub fn feasibility_check(x: &[f64], constraints: &[fn(&[f64]) -> f64], tol: f64) -> bool {
148    constraints.iter().all(|c| c(x) <= tol)
149}
150
151pub fn quadratic_objective(h: &[Vec<f64>], c: &[f64], x: &[f64]) -> f64 {
152    let mut val = 0.0;
153    for (i, &xi) in x.iter().enumerate() {
154        val += c[i] * xi;
155        for (j, &xj) in x.iter().enumerate() {
156            val += 0.5 * h[i][j] * xi * xj;
157        }
158    }
159    val
160}
161
162pub fn linear_constraint_violation(a: &[Vec<f64>], b: &[f64], x: &[f64]) -> f64 {
163    let n = x.len();
164    let mut max_viol = 0.0f64;
165    for (i, ai) in a.iter().enumerate() {
166        let mut dot = 0.0;
167        for j in 0..n {
168            dot += ai[j] * x[j];
169        }
170        let viol = (dot - b[i]).max(0.0);
171        if viol > max_viol {
172            max_viol = viol;
173        }
174    }
175    max_viol
176}
177
178pub fn l1_penalty(
179    f: fn(&[f64]) -> f64,
180    constraints: &[fn(&[f64]) -> f64],
181    x: &[f64],
182    penalty: f64,
183) -> f64 {
184    let mut cost = f(x);
185    for c in constraints {
186        cost += penalty * c(x).max(0.0);
187    }
188    cost
189}
190
191pub fn equality_penalty(
192    f: fn(&[f64]) -> f64,
193    eq_constraints: &[fn(&[f64]) -> f64],
194    x: &[f64],
195    penalty: f64,
196) -> f64 {
197    let mut cost = f(x);
198    for c in eq_constraints {
199        let v = c(x);
200        cost += penalty * v * v;
201    }
202    cost
203}
204
205pub fn merit_function(f: f64, constraints: &[f64], mu: f64) -> f64 {
206    let violation: f64 = constraints.iter().map(|&c| c.max(0.0)).sum();
207    f + mu * violation
208}
sciforge_lib/maths/optimization/constrained.rs

sciforge_lib/maths/optimization/
constrained.rs
