```  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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
```
```//! A port of the proximal gradient method from `https://github.com/pmelchior/proxmin`.
//!
#![allow(unused)]
#![allow(non_snake_case)]
#![allow(non_upper_case_globals)]

mod examples;
pub(crate) mod misc;
pub mod utils;

use crate::{
misc::*,
utils::NesterovStepper,
};

#[derive(Copy, Clone, Debug)]
pub enum ProximalOptimizerErr {
/// Caused when the length of parameter vectors does not match the number
/// of parameters specified when the optimizer was created.
ParameterLengthMismatch,
/// The objective function's value at the start position could not be compared
/// with itself (perhaps a `NaN` condition?)
StartUnorderable,
/// The candidate solution is no better than the starting position, and
/// may actually be worse.
SolutionNoBetter,
}

/// Proximal Gradient Method (PGM), ported
/// from `https://github.com/pmelchior/proxmin`
///
///
///  Adapted from Combettes 2009, Algorithm 3.4.
///  The accelerated version is Algorithm 3.6 with modifications
///  from Xu & Yin (2015).
///
///  Args:
///  - start_x: starting position
///  - prox_f: proxed function f (the forward-backward step)
///  - step_f: step size, < 1/L with L being the Lipschitz constant of grad f
///  - accelerated: If Nesterov acceleration should be used
///  - relax: (over)relaxation parameter, 0 < relax < 1.5
///  - e_rel: relative error of X
///  - max_iter: maximum iteration, irrespective of residual error
///  - traceback: utils.Traceback to hold variable histories
///
///  Returns: A 3-tuple containing:
///  - The optimized value for X
///  - converged: whether the optimizer has converged within e_rel
///  - error: X^it - X^it-1
fn pgm<F>(start_x: &[f64],
prox_f: F,
step_f: &[f64],
accelerated: bool,
relax: Option<f64>,
e_rel: f64,
max_iter: usize)
-> Result<(Vec<f64>, bool, Vec<f64>), ProximalOptimizerErr>
where F: Fn(&[f64], &[f64]) -> Vec<f64>
{
let mut stepper = NesterovStepper::new(accelerated);

if let Some(relax_val) = relax {
assert!(relax_val > 0.0);
assert!(relax_val < 1.5);
}

let mut X = Vec::from(&start_x[..]);
let mut X_ = vec![0.0; start_x.len()];

let mut it: usize = 0;
let mut converged: bool = false;
while it < max_iter {
let _X;

// use Nesterov acceleration (if omega > 0), automatically incremented
let omega = stepper.omega();
log::trace!("Omega: {}", &omega);
if omega > 0.0 {
// In Python: _X = X + omega*(X - X_)
let tmp1 = vec_sub(&X[..], &X_[..])?;
let tmp2 = vec_mul_scalar(&tmp1[..], omega);
} else {
_X = X.clone();
}

log::trace!("_X: {:?}", &_X);

X_ = X.clone();

X = prox_f(&_X[..], step_f);

log::trace!("X: {:?}", &X);

if let Some(relax_val) = relax {
// In Python: X += (relax-1)*(X - X_)
let tmp1 = relax_val - 1.0;
let tmp2 = vec_sub(&X[..], &X_[..])?;
let tmp3 = vec_mul_scalar(&tmp2[..], tmp1);
}

// test for fixed point convergence
// In Python: converged = utils.l2sq(X - X_) <= e_rel**2*utils.l2sq(X)
let tmp1 = vec_sub(&X[..], &X_[..])?;
let left = utils::l2sq(&tmp1[..]);
let right = utils::l2sq(&X[..]) * e_rel * e_rel;
converged = left <= right;
if converged {
break;
}
it += 1;
}

log::info!("Completed {} iterations", it + 1);
if !converged {
log::warn!("Solution did not converge");
}

let error = vec_sub(&X[..], &X_[..])?;

return Ok((X, converged, error));
}
```