Skip to content

Make core generic over float types #373

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from
Draft

Make core generic over float types #373

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 5 additions & 1 deletion Cargo.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ rustdoc-args = ["--html-in-header", "katex-header.html"]
num = "0.4"

# Our own stuff - L-BFGS: limited-memory BFGS directions
lbfgs = "0.2"
lbfgs = { version = "0.2" }

# Instant is a generic timer that works on Wasm (with wasm-bindgen)
instant = { version = "0.1" }
Expand Down Expand Up @@ -147,3 +147,7 @@ travis-ci = { repository = "alphaville/optimization-engine", branch = "master" }

# Actively maintained badge
maintenance = { status = "actively-developed" }

[patch.crates-io]
# Add support for generic floating point types
lbfgs = { git = "https://github.com/AutoPallet/lbfgs-rs.git", rev = "9f450872cb2da6bc7d3ebb79f32e29550689ecb0" }
73 changes: 46 additions & 27 deletions src/alm/alm_cache.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
use crate::panoc::PANOCCache;
use crate::OptFloat;

const DEFAULT_INITIAL_PENALTY: f64 = 10.0;

Expand All @@ -12,40 +13,46 @@ const DEFAULT_INITIAL_PENALTY: f64 = 10.0;
/// of `AlmProblem`
///
#[derive(Debug)]
pub struct AlmCache {
pub struct AlmCache<T>
where
T: OptFloat,
{
/// PANOC cache for inner problems
pub(crate) panoc_cache: PANOCCache,
pub(crate) panoc_cache: PANOCCache<T>,
/// Lagrange multipliers (next)
pub(crate) y_plus: Option<Vec<f64>>,
pub(crate) y_plus: Option<Vec<T>>,
/// Vector $\xi^\nu = (c^\nu, y^\nu)$
pub(crate) xi: Option<Vec<f64>>,
pub(crate) xi: Option<Vec<T>>,
/// Infeasibility related to ALM-type constraints
pub(crate) delta_y_norm: f64,
pub(crate) delta_y_norm: T,
/// Delta y at iteration `nu+1`
pub(crate) delta_y_norm_plus: f64,
pub(crate) delta_y_norm_plus: T,
/// Value $\Vert F_2(u^\nu) \Vert$
pub(crate) f2_norm: f64,
pub(crate) f2_norm: T,
/// Value $\Vert F_2(u^{\nu+1}) \Vert$
pub(crate) f2_norm_plus: f64,
pub(crate) f2_norm_plus: T,
/// Auxiliary variable `w`
pub(crate) w_alm_aux: Option<Vec<f64>>,
pub(crate) w_alm_aux: Option<Vec<T>>,
/// Infeasibility related to PM-type constraints, `w_pm = F2(u)`
pub(crate) w_pm: Option<Vec<f64>>,
pub(crate) w_pm: Option<Vec<T>>,
/// (Outer) iteration count
pub(crate) iteration: usize,
/// Counter for inner iterations
pub(crate) inner_iteration_count: usize,
/// Value of the norm of the fixed-point residual for the last
/// solved inner problem
pub(crate) last_inner_problem_norm_fpr: f64,
pub(crate) last_inner_problem_norm_fpr: T,
/// Available time left for ALM/PM computations (the value `None`
/// corresponds to an unspecified available time, i.e., there are
/// no bounds on the maximum time). The maximum time is specified,
/// if at all, in `AlmOptimizer`
pub(crate) available_time: Option<std::time::Duration>,
}

impl AlmCache {
impl<T> AlmCache<T>
where
T: OptFloat,
{
/// Construct a new instance of `AlmCache`
///
/// # Arguments
Expand All @@ -58,30 +65,42 @@ impl AlmCache {
///
/// Does not panic
///
pub fn new(panoc_cache: PANOCCache, n1: usize, n2: usize) -> Self {
pub fn new(panoc_cache: PANOCCache<T>, n1: usize, n2: usize) -> Self {
AlmCache {
panoc_cache,
y_plus: if n1 > 0 { Some(vec![0.0; n1]) } else { None },
y_plus: if n1 > 0 {
Some(vec![T::zero(); n1])
} else {
None
},
// Allocate memory for xi = (c, y) if either n1 or n2 is nonzero,
// otherwise, xi is None
xi: if n1 + n2 > 0 {
let mut xi_init = vec![DEFAULT_INITIAL_PENALTY; 1];
xi_init.append(&mut vec![0.0; n1]);
let mut xi_init = vec![T::from(DEFAULT_INITIAL_PENALTY).unwrap(); 1];
xi_init.append(&mut vec![T::zero(); n1]);
Some(xi_init)
} else {
None
},
// w_alm_aux should be allocated only if n1 > 0
w_alm_aux: if n1 > 0 { Some(vec![0.0; n1]) } else { None },
w_alm_aux: if n1 > 0 {
Some(vec![T::zero(); n1])
} else {
None
},
// w_pm is needed only if n2 > 0
w_pm: if n2 > 0 { Some(vec![0.0; n2]) } else { None },
w_pm: if n2 > 0 {
Some(vec![T::zero(); n2])
} else {
None
},
iteration: 0,
delta_y_norm: 0.0,
delta_y_norm_plus: std::f64::INFINITY,
f2_norm: 0.0,
f2_norm_plus: std::f64::INFINITY,
delta_y_norm: T::zero(),
delta_y_norm_plus: T::infinity(),
f2_norm: T::zero(),
f2_norm_plus: T::infinity(),
inner_iteration_count: 0,
last_inner_problem_norm_fpr: -1.0,
last_inner_problem_norm_fpr: T::from(-1.0).unwrap(),
available_time: None,
}
}
Expand All @@ -92,10 +111,10 @@ impl AlmCache {
pub fn reset(&mut self) {
self.panoc_cache.reset();
self.iteration = 0;
self.f2_norm = 0.0;
self.f2_norm_plus = 0.0;
self.delta_y_norm = 0.0;
self.delta_y_norm_plus = 0.0;
self.f2_norm = T::zero();
self.f2_norm_plus = T::zero();
self.delta_y_norm = T::zero();
self.delta_y_norm_plus = T::zero();
self.inner_iteration_count = 0;
}
}
75 changes: 42 additions & 33 deletions src/alm/alm_factory.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,10 +1,10 @@
use crate::OptFloat;
/* ---------------------------------------------------------------------------- */
/* ALM FACTORY */
/* */
/* About: The user provides f, df, F1, JF1'*d, F2 and C and MockAlmFactory */
/* prepares psi and d_psi, which can be used to define an AlmOptimizer */
/* ---------------------------------------------------------------------------- */

use crate::{constraints::Constraint, matrix_operations, FunctionCallResult};

/// Prepares function $\psi$ and its gradient given the problem data: $f$, $\nabla{}f$,
Expand Down Expand Up @@ -72,14 +72,16 @@ pub struct AlmFactory<
Cost,
CostGradient,
SetC,
T,
> where
Cost: Fn(&[f64], &mut f64) -> FunctionCallResult, // f(u, result)
CostGradient: Fn(&[f64], &mut [f64]) -> FunctionCallResult, // df(u, result)
MappingF1: Fn(&[f64], &mut [f64]) -> FunctionCallResult, // f1(u, result)
JacobianMappingF1Trans: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, // jf1(u, d, result)
MappingF2: Fn(&[f64], &mut [f64]) -> FunctionCallResult, // f2(u, result)
JacobianMappingF2Trans: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, // jf2(u, d, result)
SetC: Constraint,
Cost: Fn(&[T], &mut T) -> FunctionCallResult, // f(u, result)
CostGradient: Fn(&[T], &mut [T]) -> FunctionCallResult, // df(u, result)
MappingF1: Fn(&[T], &mut [T]) -> FunctionCallResult, // f1(u, result)
JacobianMappingF1Trans: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, // jf1(u, d, result)
MappingF2: Fn(&[T], &mut [T]) -> FunctionCallResult, // f2(u, result)
JacobianMappingF2Trans: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, // jf2(u, d, result)
SetC: Constraint<T>,
T: OptFloat,
{
f: Cost,
df: CostGradient,
Expand All @@ -89,6 +91,7 @@ pub struct AlmFactory<
jacobian_mapping_f2_trans: Option<JacobianMappingF2Trans>,
set_c: Option<SetC>,
n2: usize,
_t: std::marker::PhantomData<T>,
}

impl<
Expand All @@ -99,6 +102,7 @@ impl<
Cost,
CostGradient,
SetC,
T,
>
AlmFactory<
MappingF1,
Expand All @@ -108,15 +112,17 @@ impl<
Cost,
CostGradient,
SetC,
T,
>
where
Cost: Fn(&[f64], &mut f64) -> FunctionCallResult, // f(u, result)
CostGradient: Fn(&[f64], &mut [f64]) -> FunctionCallResult, // df(u, result)
MappingF1: Fn(&[f64], &mut [f64]) -> FunctionCallResult, // f1(u, result)
JacobianMappingF1Trans: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, // jf1(u, d, result)
MappingF2: Fn(&[f64], &mut [f64]) -> FunctionCallResult, // f2(u, result)
JacobianMappingF2Trans: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, // jf2(u, d, result)
SetC: Constraint,
Cost: Fn(&[T], &mut T) -> FunctionCallResult, // f(u, result)
CostGradient: Fn(&[T], &mut [T]) -> FunctionCallResult, // df(u, result)
MappingF1: Fn(&[T], &mut [T]) -> FunctionCallResult, // f1(u, result)
JacobianMappingF1Trans: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, // jf1(u, d, result)
MappingF2: Fn(&[T], &mut [T]) -> FunctionCallResult, // f2(u, result)
JacobianMappingF2Trans: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, // jf2(u, d, result)
SetC: Constraint<T>,
T: OptFloat,
{
/// Construct a new instance of `MockFactory`
///
Expand Down Expand Up @@ -190,6 +196,7 @@ where
jacobian_mapping_f2_trans,
set_c,
n2,
_t: std::marker::PhantomData,
}
}

Expand All @@ -215,11 +222,11 @@ where
/// This method returns `Ok(())` if the computation is successful or an appropriate
/// `SolverError` otherwise.
///
pub fn psi(&self, u: &[f64], xi: &[f64], cost: &mut f64) -> FunctionCallResult {
pub fn psi(&self, u: &[T], xi: &[T], cost: &mut T) -> FunctionCallResult {
(self.f)(u, cost)?;
let ny = if !xi.is_empty() { xi.len() - 1 } else { 0 };
let mut f1_u_plus_y_over_c = vec![0.0; ny];
let mut s = vec![0.0; ny];
let mut f1_u_plus_y_over_c = vec![T::zero(); ny];
let mut s = vec![T::zero(); ny];
if let (Some(set_c), Some(mapping_f1)) = (&self.set_c, &self.mapping_f1) {
let penalty_parameter = xi[0];
mapping_f1(u, &mut f1_u_plus_y_over_c)?; // f1_u = F1(u)
Expand All @@ -231,18 +238,18 @@ where
f1_u_plus_y_over_c
.iter_mut()
.zip(y_lagrange_mult.iter())
.for_each(|(ti, yi)| *ti += yi / f64::max(penalty_parameter, 1.0));
.for_each(|(ti, yi)| *ti += *yi / T::max(penalty_parameter, T::one()));
s.copy_from_slice(&f1_u_plus_y_over_c);
set_c.project(&mut s);
*cost += 0.5
*cost += T::from(0.5).unwrap()
* penalty_parameter
* matrix_operations::norm2_squared_diff(&f1_u_plus_y_over_c, &s);
}
if let Some(f2) = &self.mapping_f2 {
let c = xi[0];
let mut z = vec![0.0; self.n2];
let mut z = vec![T::zero(); self.n2];
f2(u, &mut z)?;
*cost += 0.5 * c * matrix_operations::norm2_squared(&z);
*cost += T::from(0.5).unwrap() * c * matrix_operations::norm2_squared(&z);
}
Ok(())
}
Expand All @@ -267,7 +274,7 @@ where
/// This method returns `Ok(())` if the computation is successful or an appropriate
/// `SolverError` otherwise.
///
pub fn d_psi(&self, u: &[f64], xi: &[f64], grad: &mut [f64]) -> FunctionCallResult {
pub fn d_psi(&self, u: &[T], xi: &[T], grad: &mut [T]) -> FunctionCallResult {
let nu = u.len();

// The following statement is needed to account for the case where n1=n2=0
Expand All @@ -285,54 +292,56 @@ where
&self.jacobian_mapping_f1_trans,
) {
let c_penalty_parameter = xi[0];
let mut f1_u_plus_y_over_c = vec![0.0; ny];
let mut s_aux_var = vec![0.0; ny]; // auxiliary variable `s`
let mut f1_u_plus_y_over_c = vec![T::zero(); ny];
let mut s_aux_var = vec![T::zero(); ny]; // auxiliary variable `s`
let y_lagrange_mult = &xi[1..];
let mut jac_prod = vec![0.0; nu];
let mut jac_prod = vec![T::zero(); nu];
mapping_f1(u, &mut f1_u_plus_y_over_c)?; // f1_u_plus_y_over_c = F1(u)
// f1_u_plus_y_over_c = F1(u) + y/c
f1_u_plus_y_over_c
.iter_mut()
.zip(y_lagrange_mult.iter())
.for_each(|(ti, yi)| *ti += yi / c_penalty_parameter);
.for_each(|(ti, yi)| *ti += *yi / c_penalty_parameter);
s_aux_var.copy_from_slice(&f1_u_plus_y_over_c); // s = t
set_c.project(&mut s_aux_var); // s = Proj_C(F1(u) + y/c)

// t = F1(u) + y/c - Proj_C(F1(u) + y/c)
f1_u_plus_y_over_c
.iter_mut()
.zip(s_aux_var.iter())
.for_each(|(ti, si)| *ti -= si);
.for_each(|(ti, si)| *ti -= *si);

jf1t(u, &f1_u_plus_y_over_c, &mut jac_prod)?;

// grad += c*t
grad.iter_mut()
.zip(jac_prod.iter())
.for_each(|(gradi, jac_prodi)| *gradi += c_penalty_parameter * jac_prodi);
.for_each(|(gradi, jac_prodi)| *gradi += c_penalty_parameter * *jac_prodi);
}

// Compute second part: JF2(u)'*F2(u)
if let (Some(f2), Some(jf2)) = (&self.mapping_f2, &self.jacobian_mapping_f2_trans) {
let c = xi[0];
let mut f2u_aux = vec![0.0; self.n2];
let mut jf2u_times_f2u_aux = vec![0.0; nu];
let mut f2u_aux = vec![T::zero(); self.n2];
let mut jf2u_times_f2u_aux = vec![T::zero(); nu];
f2(u, &mut f2u_aux)?; // f2u_aux = F2(u)
jf2(u, &f2u_aux, &mut jf2u_times_f2u_aux)?; // jf2u_times_f2u_aux = JF2(u)'*f2u_aux
// = JF2(u)'*F2(u)

// grad += c * jf2u_times_f2u_aux
grad.iter_mut()
.zip(jf2u_times_f2u_aux.iter())
.for_each(|(gradi, jf2u_times_f2u_aux_i)| *gradi += c * jf2u_times_f2u_aux_i);
.for_each(|(gradi, jf2u_times_f2u_aux_i)| *gradi += c * *jf2u_times_f2u_aux_i);
}
Ok(())
}
}

#[cfg(test)]
mod tests {
use crate::{alm::*, constraints::*, mocks, FunctionCallResult, SolverError};
use crate::alm::*;
use crate::constraints::*;
use crate::{mocks, FunctionCallResult, SolverError};

#[test]
fn t_mocking_alm_factory_psi() {
Expand Down
Loading