From 7bedc3d62a5421127758a084973b82335138e2a1 Mon Sep 17 00:00:00 2001 From: Danny Wolf Date: Tue, 12 Aug 2025 01:25:22 +0000 Subject: [PATCH] Make core generic over float types - Make a OptFloat trait - Provide PANOC constants implemented for f32 and f64 - Use somewhat arbitrary values for f32 PANOC constants, unclear how good these are. These should be looked at. - Add some tests for f32 (some of which are failing!) --- Cargo.toml | 6 +- src/alm/alm_cache.rs | 73 +++++--- src/alm/alm_factory.rs | 75 ++++---- src/alm/alm_optimizer.rs | 197 ++++++++++---------- src/alm/alm_optimizer_status.rs | 78 ++++---- src/alm/alm_problem.rs | 39 ++-- src/alm/tests.rs | 10 +- src/constraints/affine_space.rs | 40 ++-- src/constraints/ball1.rs | 39 ++-- src/constraints/ball2.rs | 32 ++-- src/constraints/ballinf.rs | 28 ++- src/constraints/cartesian_product.rs | 23 ++- src/constraints/epigraph_squared_norm.rs | 60 +++--- src/constraints/finite.rs | 25 ++- src/constraints/halfspace.rs | 28 ++- src/constraints/hyperplane.rs | 28 ++- src/constraints/mod.rs | 9 +- src/constraints/no_constraints.rs | 9 +- src/constraints/rectangle.rs | 25 ++- src/constraints/simplex.rs | 46 +++-- src/constraints/soc.rs | 30 ++- src/constraints/sphere2.rs | 31 ++-- src/constraints/tests.rs | 13 +- src/constraints/zero.rs | 11 +- src/core/fbs/fbs_cache.rs | 30 +-- src/core/fbs/fbs_engine.rs | 57 +++--- src/core/fbs/fbs_optimizer.rs | 65 +++---- src/core/fbs/tests.rs | 3 +- src/core/mod.rs | 18 +- src/core/opt_float.rs | 113 ++++++++++++ src/core/panoc/mod.rs | 3 + src/core/panoc/panoc_cache.rs | 119 ++++++------ src/core/panoc/panoc_engine.rs | 144 +++++++-------- src/core/panoc/panoc_optimizer.rs | 63 ++++--- src/core/panoc/tests.rs | 22 ++- src/core/panoc/tests_f32.rs | 222 +++++++++++++++++++++++ src/core/problem.rs | 26 ++- src/core/solver_status.rs | 26 ++- src/lib.rs | 6 +- src/lipschitz_estimator.rs | 46 ++--- src/matrix_operations.rs | 7 +- src/mocks.rs | 143 +++++++++------ src/tests.rs | 31 +++- 43 files changed, 1368 insertions(+), 731 deletions(-) create mode 100644 src/core/opt_float.rs create mode 100644 src/core/panoc/tests_f32.rs diff --git a/Cargo.toml b/Cargo.toml index b967529b..56e4367a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -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" } @@ -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" } diff --git a/src/alm/alm_cache.rs b/src/alm/alm_cache.rs index e014c3bc..8d7e6954 100644 --- a/src/alm/alm_cache.rs +++ b/src/alm/alm_cache.rs @@ -1,4 +1,5 @@ use crate::panoc::PANOCCache; +use crate::OptFloat; const DEFAULT_INITIAL_PENALTY: f64 = 10.0; @@ -12,32 +13,35 @@ const DEFAULT_INITIAL_PENALTY: f64 = 10.0; /// of `AlmProblem` /// #[derive(Debug)] -pub struct AlmCache { +pub struct AlmCache +where + T: OptFloat, +{ /// PANOC cache for inner problems - pub(crate) panoc_cache: PANOCCache, + pub(crate) panoc_cache: PANOCCache, /// Lagrange multipliers (next) - pub(crate) y_plus: Option>, + pub(crate) y_plus: Option>, /// Vector $\xi^\nu = (c^\nu, y^\nu)$ - pub(crate) xi: Option>, + pub(crate) xi: Option>, /// 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>, + pub(crate) w_alm_aux: Option>, /// Infeasibility related to PM-type constraints, `w_pm = F2(u)` - pub(crate) w_pm: Option>, + pub(crate) w_pm: Option>, /// (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, @@ -45,7 +49,10 @@ pub struct AlmCache { pub(crate) available_time: Option, } -impl AlmCache { +impl AlmCache +where + T: OptFloat, +{ /// Construct a new instance of `AlmCache` /// /// # Arguments @@ -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, 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, } } @@ -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; } } diff --git a/src/alm/alm_factory.rs b/src/alm/alm_factory.rs index 735c2be5..ed1f94fd 100644 --- a/src/alm/alm_factory.rs +++ b/src/alm/alm_factory.rs @@ -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$, @@ -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: OptFloat, { f: Cost, df: CostGradient, @@ -89,6 +91,7 @@ pub struct AlmFactory< jacobian_mapping_f2_trans: Option, set_c: Option, n2: usize, + _t: std::marker::PhantomData, } impl< @@ -99,6 +102,7 @@ impl< Cost, CostGradient, SetC, + T, > AlmFactory< MappingF1, @@ -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: OptFloat, { /// Construct a new instance of `MockFactory` /// @@ -190,6 +196,7 @@ where jacobian_mapping_f2_trans, set_c, n2, + _t: std::marker::PhantomData, } } @@ -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) @@ -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(()) } @@ -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 @@ -285,16 +292,16 @@ 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) @@ -302,21 +309,21 @@ where 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) @@ -324,7 +331,7 @@ where // 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(()) } @@ -332,7 +339,9 @@ where #[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() { diff --git a/src/alm/alm_optimizer.rs b/src/alm/alm_optimizer.rs index 0fcb1fde..929ea634 100644 --- a/src/alm/alm_optimizer.rs +++ b/src/alm/alm_optimizer.rs @@ -1,9 +1,7 @@ -use crate::{ - alm::*, - constraints, - core::{panoc::PANOCOptimizer, ExitStatus, Optimizer, Problem, SolverStatus}, - matrix_operations, FunctionCallResult, SolverError, -}; +use crate::alm::*; +use crate::core::panoc::PANOCOptimizer; +use crate::core::{ExitStatus, Optimizer, Problem, SolverStatus}; +use crate::{constraints, matrix_operations, FunctionCallResult, OptFloat, SolverError}; const DEFAULT_MAX_OUTER_ITERATIONS: usize = 50; const DEFAULT_MAX_INNER_ITERATIONS: usize = 5000; @@ -13,7 +11,6 @@ const DEFAULT_PENALTY_UPDATE_FACTOR: f64 = 5.0; const DEFAULT_EPSILON_UPDATE_FACTOR: f64 = 0.1; const DEFAULT_INFEAS_SUFFICIENT_DECREASE_FACTOR: f64 = 0.1; const DEFAULT_INITIAL_TOLERANCE: f64 = 0.1; -const SMALL_EPSILON: f64 = std::f64::EPSILON; /// Internal/private structure used by method AlmOptimizer.step /// to return some minimal information about the inner problem @@ -124,17 +121,19 @@ pub struct AlmOptimizer< ConstraintsType, AlmSetC, LagrangeSetY, + T, > where - MappingAlm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - MappingPm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - ParametricGradientType: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, - ParametricCostType: Fn(&[f64], &[f64], &mut f64) -> FunctionCallResult, - ConstraintsType: constraints::Constraint, - AlmSetC: constraints::Constraint, - LagrangeSetY: constraints::Constraint, + MappingAlm: Fn(&[T], &mut [T]) -> FunctionCallResult, + MappingPm: Fn(&[T], &mut [T]) -> FunctionCallResult, + ParametricGradientType: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, + ParametricCostType: Fn(&[T], &[T], &mut T) -> FunctionCallResult, + ConstraintsType: constraints::Constraint, + AlmSetC: constraints::Constraint, + LagrangeSetY: constraints::Constraint, + T: OptFloat, { /// ALM cache (borrowed) - alm_cache: &'life mut AlmCache, + alm_cache: &'life mut AlmCache, /// ALM problem definition (oracle) alm_problem: AlmProblem< MappingAlm, @@ -144,6 +143,7 @@ pub struct AlmOptimizer< ConstraintsType, AlmSetC, LagrangeSetY, + T, >, /// Maximum number of outer iterations max_outer_iterations: usize, @@ -152,19 +152,19 @@ pub struct AlmOptimizer< /// Maximum duration max_duration: Option, /// epsilon for inner AKKT condition - epsilon_tolerance: f64, + epsilon_tolerance: T, /// delta for outer AKKT condition - delta_tolerance: f64, + delta_tolerance: T, /// At every outer iteration, c is multiplied by this scalar - penalty_update_factor: f64, + penalty_update_factor: T, /// The epsilon-tolerance is multiplied by this factor until /// it reaches its target value - epsilon_update_factor: f64, + epsilon_update_factor: T, /// If current_infeasibility <= sufficient_decrease_coeff * previous_infeasibility, /// then the penalty parameter is kept constant - sufficient_decrease_coeff: f64, + sufficient_decrease_coeff: T, // Initial tolerance (for the inner problem) - epsilon_inner_initial: f64, + epsilon_inner_initial: T, } impl< @@ -176,6 +176,7 @@ impl< ConstraintsType, AlmSetC, LagrangeSetY, + T, > AlmOptimizer< 'life, @@ -186,15 +187,17 @@ impl< ConstraintsType, AlmSetC, LagrangeSetY, + T, > where - MappingAlm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - MappingPm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - ParametricGradientType: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, - ParametricCostType: Fn(&[f64], &[f64], &mut f64) -> FunctionCallResult, - ConstraintsType: constraints::Constraint, - AlmSetC: constraints::Constraint, - LagrangeSetY: constraints::Constraint, + MappingAlm: Fn(&[T], &mut [T]) -> FunctionCallResult, + MappingPm: Fn(&[T], &mut [T]) -> FunctionCallResult, + ParametricGradientType: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, + ParametricCostType: Fn(&[T], &[T], &mut T) -> FunctionCallResult, + ConstraintsType: constraints::Constraint, + AlmSetC: constraints::Constraint, + LagrangeSetY: constraints::Constraint, + T: OptFloat + std::fmt::Debug, { /* ---------------------------------------------------------------------------- */ /* CONSTRUCTOR */ @@ -214,7 +217,7 @@ where /// # Example /// /// ```rust - /// use optimization_engine::{alm::*, FunctionCallResult, core::{panoc::*, constraints}}; + /// use optimization_engine::{alm::*, FunctionCallResult, core::{OptFloat, panoc::*, constraints}}; /// /// let tolerance = 1e-8; /// let nx = 10; @@ -250,7 +253,7 @@ where ///``` /// pub fn new( - alm_cache: &'life mut AlmCache, + alm_cache: &'life mut AlmCache, alm_problem: AlmProblem< MappingAlm, MappingPm, @@ -259,6 +262,7 @@ where ConstraintsType, AlmSetC, LagrangeSetY, + T, >, ) -> Self { // set the initial value of the inner tolerance; this step is @@ -266,19 +270,19 @@ where // in #solve (see below) alm_cache .panoc_cache - .set_akkt_tolerance(DEFAULT_INITIAL_TOLERANCE); + .set_akkt_tolerance(T::from(DEFAULT_INITIAL_TOLERANCE).unwrap()); AlmOptimizer { alm_cache, alm_problem, max_outer_iterations: DEFAULT_MAX_OUTER_ITERATIONS, max_inner_iterations: DEFAULT_MAX_INNER_ITERATIONS, max_duration: None, - epsilon_tolerance: DEFAULT_EPSILON_TOLERANCE, - delta_tolerance: DEFAULT_DELTA_TOLERANCE, - penalty_update_factor: DEFAULT_PENALTY_UPDATE_FACTOR, - epsilon_update_factor: DEFAULT_EPSILON_UPDATE_FACTOR, - sufficient_decrease_coeff: DEFAULT_INFEAS_SUFFICIENT_DECREASE_FACTOR, - epsilon_inner_initial: DEFAULT_INITIAL_TOLERANCE, + epsilon_tolerance: T::from(DEFAULT_EPSILON_TOLERANCE).unwrap(), + delta_tolerance: T::from(DEFAULT_DELTA_TOLERANCE).unwrap(), + penalty_update_factor: T::from(DEFAULT_PENALTY_UPDATE_FACTOR).unwrap(), + epsilon_update_factor: T::from(DEFAULT_EPSILON_UPDATE_FACTOR).unwrap(), + sufficient_decrease_coeff: T::from(DEFAULT_INFEAS_SUFFICIENT_DECREASE_FACTOR).unwrap(), + epsilon_inner_initial: T::from(DEFAULT_INITIAL_TOLERANCE).unwrap(), } } @@ -367,8 +371,11 @@ where /// /// The method panics if the specified tolerance is not positive /// - pub fn with_delta_tolerance(mut self, delta_tolerance: f64) -> Self { - assert!(delta_tolerance > 0.0, "delta_tolerance must be positive"); + pub fn with_delta_tolerance(mut self, delta_tolerance: T) -> Self { + assert!( + delta_tolerance > T::zero(), + "delta_tolerance must be positive" + ); self.delta_tolerance = delta_tolerance; self } @@ -387,9 +394,9 @@ where /// /// The method panics if the specified tolerance is not positive /// - pub fn with_epsilon_tolerance(mut self, epsilon_tolerance: f64) -> Self { + pub fn with_epsilon_tolerance(mut self, epsilon_tolerance: T) -> Self { assert!( - epsilon_tolerance > 0.0, + epsilon_tolerance > T::zero(), "epsilon_tolerance must be positive" ); self.epsilon_tolerance = epsilon_tolerance; @@ -415,10 +422,10 @@ where /// The method panics if the update factor is not larger than `1.0 + f64::EPSILON` /// /// - pub fn with_penalty_update_factor(mut self, penalty_update_factor: f64) -> Self { + pub fn with_penalty_update_factor(mut self, penalty_update_factor: T) -> Self { assert!( - penalty_update_factor > 1.0 + SMALL_EPSILON, - "`penalty_update_factor` must be larger than 1.0 + f64::EPSILON" + penalty_update_factor > T::one() + T::epsilon(), + "`penalty_update_factor` must be larger than 1.0 + EPSILON" ); self.penalty_update_factor = penalty_update_factor; self @@ -444,14 +451,11 @@ where /// The method panics if the specified tolerance update factor is not in the /// interval from `f64::EPSILON` to `1.0 - f64::EPSILON`. /// - pub fn with_inner_tolerance_update_factor( - mut self, - inner_tolerance_update_factor: f64, - ) -> Self { + pub fn with_inner_tolerance_update_factor(mut self, inner_tolerance_update_factor: T) -> Self { assert!( - inner_tolerance_update_factor > SMALL_EPSILON - && inner_tolerance_update_factor < 1.0 - SMALL_EPSILON, - "the tolerance update factor needs to be in (f64::EPSILON, 1)" + inner_tolerance_update_factor > T::epsilon() + && inner_tolerance_update_factor < T::one() - T::epsilon(), + "the tolerance update factor needs to be in (EPSILON, 1)" ); self.epsilon_update_factor = inner_tolerance_update_factor; self @@ -480,7 +484,7 @@ where /// `with_inner_tolerance` to do so before invoking `with_initial_inner_tolerance`. /// /// - pub fn with_initial_inner_tolerance(mut self, initial_inner_tolerance: f64) -> Self { + pub fn with_initial_inner_tolerance(mut self, initial_inner_tolerance: T) -> Self { assert!( initial_inner_tolerance >= self.epsilon_tolerance, "the initial tolerance should be no less than the target tolerance" @@ -514,12 +518,12 @@ where /// pub fn with_sufficient_decrease_coefficient( mut self, - sufficient_decrease_coefficient: f64, + sufficient_decrease_coefficient: T, ) -> Self { assert!( - sufficient_decrease_coefficient < 1.0 - SMALL_EPSILON - && sufficient_decrease_coefficient > SMALL_EPSILON, - "sufficient_decrease_coefficient must be in (f64::EPSILON, 1.0 - f64::EPSILON)" + sufficient_decrease_coefficient < T::one() - T::epsilon() + && sufficient_decrease_coefficient > T::epsilon(), + "sufficient_decrease_coefficient must be in (EPSILON, 1.0 - EPSILON)" ); self.sufficient_decrease_coeff = sufficient_decrease_coefficient; self @@ -540,7 +544,7 @@ where /// /// The method will panic if the length of `y_init` is not equal to `n1` /// - pub fn with_initial_lagrange_multipliers(mut self, y_init: &[f64]) -> Self { + pub fn with_initial_lagrange_multipliers(mut self, y_init: &[T]) -> Self { let cache = &mut self.alm_cache; assert!( y_init.len() == self.alm_problem.n1, @@ -570,9 +574,9 @@ where /// The method panics if the specified initial penalty parameter is not /// larger than `f64::EPSILON` /// - pub fn with_initial_penalty(self, c0: f64) -> Self { + pub fn with_initial_penalty(self, c0: T) -> Self { assert!( - c0 > SMALL_EPSILON, + c0 > T::epsilon(), "the initial penalty must be larger than f64::EPSILON" ); if let Some(xi_in_cache) = &mut self.alm_cache.xi { @@ -596,7 +600,7 @@ where } /// Computes PM infeasibility, that is, ||F2(u)|| - fn compute_pm_infeasibility(&mut self, u: &[f64]) -> FunctionCallResult { + fn compute_pm_infeasibility(&mut self, u: &[T]) -> FunctionCallResult { let problem = &self.alm_problem; // ALM problem let cache = &mut self.alm_cache; // ALM cache @@ -613,7 +617,7 @@ where /// /// `y_plus <-- y + c*[F1(u_plus) - Proj_C(F1(u_plus) + y/c)]` /// - fn update_lagrange_multipliers(&mut self, u: &[f64]) -> FunctionCallResult { + fn update_lagrange_multipliers(&mut self, u: &[T]) -> FunctionCallResult { let problem = &self.alm_problem; // ALM problem let cache = &mut self.alm_cache; // ALM cache @@ -647,7 +651,7 @@ where .iter_mut() .zip(y.iter()) .zip(w_alm_aux.iter()) - .for_each(|((y_plus_i, y_i), w_alm_aux_i)| *y_plus_i = w_alm_aux_i + y_i / c); + .for_each(|((y_plus_i, y_i), w_alm_aux_i)| *y_plus_i = *w_alm_aux_i + *y_i / c); // Step #3: y_plus := Proj_C(y_plus) alm_set_c.project(y_plus); @@ -659,7 +663,7 @@ where .zip(w_alm_aux.iter()) .for_each(|((y_plus_i, y_i), w_alm_aux_i)| { // y_plus := y + c * (w_alm_aux - y_plus) - *y_plus_i = y_i + c * (w_alm_aux_i - *y_plus_i) + *y_plus_i = *y_i + c * (*w_alm_aux_i - *y_plus_i) }); } @@ -696,7 +700,7 @@ where /// error in solving the inner problem. /// /// - fn solve_inner_problem(&mut self, u: &mut [f64]) -> Result { + fn solve_inner_problem(&mut self, u: &mut [T]) -> Result, SolverError> { let alm_problem = &self.alm_problem; // Problem let alm_cache = &mut self.alm_cache; // ALM cache @@ -713,11 +717,11 @@ where // Construct psi and psi_grad (as functions of `u` alone); it is // psi(u) = psi(u; xi) and psi_grad(u) = phi_grad(u; xi) // psi: R^nu --> R - let psi = |u: &[f64], psi_val: &mut f64| -> FunctionCallResult { + let psi = |u: &[T], psi_val: &mut T| -> FunctionCallResult { (alm_problem.parametric_cost)(u, xi, psi_val) }; // psi_grad: R^nu --> R^nu - let psi_grad = |u: &[f64], psi_grad: &mut [f64]| -> FunctionCallResult { + let psi_grad = |u: &[T], psi_grad: &mut [T]| -> FunctionCallResult { (alm_problem.parametric_gradient)(u, xi, psi_grad) }; // define the inner problem @@ -750,7 +754,7 @@ where || if let Some(xi) = &cache.xi { let c = xi[0]; cache.iteration > 0 - && cache.delta_y_norm_plus <= c * self.delta_tolerance + SMALL_EPSILON + && cache.delta_y_norm_plus <= c * self.delta_tolerance + T::epsilon() } else { true }; @@ -758,13 +762,13 @@ where // If n2 = 0, there are no PM-type constraints, so this // criterion is automatically satisfied let criterion_2 = - problem.n2 == 0 || cache.f2_norm_plus <= self.delta_tolerance + SMALL_EPSILON; + problem.n2 == 0 || cache.f2_norm_plus <= self.delta_tolerance + T::epsilon(); // Criterion 3: epsilon_nu <= epsilon // This function will panic is there is no akkt_tolerance // This should never happen because we set the AKKT tolerance // in the constructor and can never become `None` again let criterion_3 = - cache.panoc_cache.akkt_tolerance.unwrap() <= self.epsilon_tolerance + SMALL_EPSILON; + cache.panoc_cache.akkt_tolerance.unwrap() <= self.epsilon_tolerance + T::epsilon(); criterion_1 && criterion_2 && criterion_3 } @@ -781,9 +785,9 @@ where let is_alm = problem.n1 > 0; let is_pm = problem.n2 > 0; let criterion_alm = cache.delta_y_norm_plus - <= self.sufficient_decrease_coeff * cache.delta_y_norm + SMALL_EPSILON; + <= self.sufficient_decrease_coeff * cache.delta_y_norm + T::epsilon(); let criterion_pm = - cache.f2_norm_plus <= self.sufficient_decrease_coeff * cache.f2_norm + SMALL_EPSILON; + cache.f2_norm_plus <= self.sufficient_decrease_coeff * cache.f2_norm + T::epsilon(); if is_alm && !is_pm { return criterion_alm; } else if !is_alm && is_pm { @@ -798,14 +802,14 @@ where fn update_penalty_parameter(&mut self) { let cache = &mut self.alm_cache; if let Some(xi) = &mut cache.xi { - xi[0] *= self.penalty_update_factor; + xi[0] = xi[0] * self.penalty_update_factor; } } fn update_inner_akkt_tolerance(&mut self) { let cache = &mut self.alm_cache; // epsilon_{nu+1} := max(epsilon, beta*epsilon_nu) - cache.panoc_cache.set_akkt_tolerance(f64::max( + cache.panoc_cache.set_akkt_tolerance(T::max( cache.panoc_cache.akkt_tolerance.unwrap() * self.epsilon_update_factor, self.epsilon_tolerance, )); @@ -837,7 +841,7 @@ where /// - Shrinks the inner tolerance and /// - Updates the ALM cache /// - fn step(&mut self, u: &mut [f64]) -> Result { + fn step(&mut self, u: &mut [T]) -> Result { // store the exit status of the inner problem in this problem // (we'll need to return it within `InnerProblemStatus`) let mut inner_exit_status: ExitStatus = ExitStatus::Converged; @@ -848,7 +852,7 @@ where // If the inner problem fails miserably, the failure should be propagated // upstream (using `?`). If the inner problem has not converged, that is fine, // we should keep solving. - self.solve_inner_problem(u).map(|status: SolverStatus| { + self.solve_inner_problem(u).map(|status: SolverStatus| { let inner_iters = status.iterations(); self.alm_cache.last_inner_problem_norm_fpr = status.norm_fpr(); self.alm_cache.inner_iteration_count += inner_iters; @@ -885,18 +889,18 @@ where Ok(InnerProblemStatus::new(true, inner_exit_status)) // `true` means do continue the outer iterations } - fn compute_cost_at_solution(&mut self, u: &mut [f64]) -> Result { + fn compute_cost_at_solution(&mut self, u: &mut [T]) -> Result { /* WORK IN PROGRESS */ let alm_problem = &self.alm_problem; // Problem let alm_cache = &mut self.alm_cache; // ALM Cache let mut empty_vec = std::vec::Vec::new(); // Empty vector - let xi: &mut std::vec::Vec = alm_cache.xi.as_mut().unwrap_or(&mut empty_vec); - let mut __c: f64 = 0.0; + let xi: &mut std::vec::Vec = alm_cache.xi.as_mut().unwrap_or(&mut empty_vec); + let mut __c: T = T::zero(); if !xi.is_empty() { __c = xi[0]; - xi[0] = 0.0; + xi[0] = T::zero(); } - let mut cost_value: f64 = 0.0; + let mut cost_value: T = T::zero(); (alm_problem.parametric_cost)(u, xi, &mut cost_value)?; if !xi.is_empty() { xi[0] = __c; @@ -911,7 +915,10 @@ where /// Solve the specified ALM problem /// /// - pub fn solve(&mut self, u: &mut [f64]) -> Result { + pub fn solve(&mut self, u: &mut [T]) -> Result, SolverError> + where + T: OptFloat + std::fmt::Debug, + { let mut num_outer_iterations = 0; // let tic = std::time::Instant::now(); let tic = instant::Instant::now(); @@ -965,7 +972,7 @@ where let c = if let Some(xi) = &self.alm_cache.xi { xi[0] } else { - 0.0 + T::zero() }; let cost = self.compute_cost_at_solution(u)?; @@ -998,13 +1005,12 @@ where #[cfg(test)] mod tests { - use crate::{ - alm::*, - core::{constraints::*, panoc::*, ExitStatus}, - matrix_operations, - mocks::*, - FunctionCallResult, - }; + use crate::alm::*; + use crate::core::constraints::*; + use crate::core::panoc::*; + use crate::core::ExitStatus; + use crate::mocks::*; + use crate::{matrix_operations, FunctionCallResult}; fn make_dummy_alm_problem( n1: usize, @@ -1014,9 +1020,10 @@ mod tests { impl Fn(&[f64], &mut [f64]) -> FunctionCallResult, impl Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, impl Fn(&[f64], &[f64], &mut f64) -> FunctionCallResult, - impl Constraint, - impl Constraint, - impl Constraint, + impl Constraint, + impl Constraint, + impl Constraint, + f64, > { // Main problem data let psi = void_parameteric_cost; @@ -1031,12 +1038,12 @@ mod tests { let set_c = if n1 > 0 { Some(Ball2::new(None, 1.50)) } else { - None:: + None::> }; - let set_y: Option = if n1 > 0 { + let set_y: Option> = if n1 > 0 { Some(Ball2::new(None, 2.0)) } else { - None:: + None::> }; // Penalty-type data let f2: Option = if n2 == 0 { diff --git a/src/alm/alm_optimizer_status.rs b/src/alm/alm_optimizer_status.rs index e6c9bbed..fc2ddeeb 100644 --- a/src/alm/alm_optimizer_status.rs +++ b/src/alm/alm_optimizer_status.rs @@ -1,5 +1,4 @@ -use crate::core::ExitStatus; - +use crate::core::{ExitStatus, OptFloat}; /// Solution statistics for `AlmOptimizer` /// /// This structure has no public fields and no public setter methods. @@ -7,7 +6,10 @@ use crate::core::ExitStatus; /// `AlmOptimizerStatus` instances. /// #[derive(Debug)] -pub struct AlmOptimizerStatus { +pub struct AlmOptimizerStatus +where + T: OptFloat, +{ /// Exit status exit_status: ExitStatus, /// Number of outer iterations @@ -18,23 +20,26 @@ pub struct AlmOptimizerStatus { /// inner solvers num_inner_iterations: usize, /// Norm of the fixed-point residual of the the problem - last_problem_norm_fpr: f64, + last_problem_norm_fpr: T, /// - lagrange_multipliers: Option>, + lagrange_multipliers: Option>, /// Total solve time solve_time: std::time::Duration, /// Last value of penalty parameter - penalty: f64, + penalty: T, /// A measure of infeasibility of constraints F1(u; p) in C - delta_y_norm: f64, + delta_y_norm: T, /// Norm of F2 at the solution, which is a measure of infeasibility /// of constraints F2(u; p) = 0 - f2_norm: f64, + f2_norm: T, /// Value of cost function at optimal solution (optimal cost) - cost: f64, + cost: T, } -impl AlmOptimizerStatus { +impl AlmOptimizerStatus +where + T: OptFloat, +{ /// Constructor for instances of `AlmOptimizerStatus` /// /// This method is only accessibly within this crate. @@ -60,13 +65,13 @@ impl AlmOptimizerStatus { exit_status, num_outer_iterations: 0, num_inner_iterations: 0, - last_problem_norm_fpr: -1.0, + last_problem_norm_fpr: T::from(-1.0).unwrap(), lagrange_multipliers: None, solve_time: std::time::Duration::from_nanos(0), - penalty: 0.0, - delta_y_norm: 0.0, - f2_norm: 0.0, - cost: 0.0, + penalty: T::zero(), + delta_y_norm: T::zero(), + f2_norm: T::zero(), + cost: T::zero(), } } @@ -129,7 +134,7 @@ impl AlmOptimizerStatus { /// Does not panic; it is the responsibility of the caller to provide a vector of /// Lagrange multipliers of correct length /// - pub(crate) fn with_lagrange_multipliers(mut self, lagrange_multipliers: &[f64]) -> Self { + pub(crate) fn with_lagrange_multipliers(mut self, lagrange_multipliers: &[T]) -> Self { self.lagrange_multipliers = Some(vec![]); if let Some(y) = &mut self.lagrange_multipliers { y.extend_from_slice(lagrange_multipliers); @@ -144,9 +149,9 @@ impl AlmOptimizerStatus { /// /// The method panics if the provided penalty parameter is negative /// - pub(crate) fn with_penalty(mut self, penalty: f64) -> Self { + pub(crate) fn with_penalty(mut self, penalty: T) -> Self { assert!( - penalty >= 0.0, + penalty >= T::zero(), "the penalty parameter should not be negative" ); self.penalty = penalty; @@ -161,28 +166,31 @@ impl AlmOptimizerStatus { /// The method panics if the provided norm of the fixed-point residual is /// negative /// - pub(crate) fn with_last_problem_norm_fpr(mut self, last_problem_norm_fpr: f64) -> Self { + pub(crate) fn with_last_problem_norm_fpr(mut self, last_problem_norm_fpr: T) -> Self { assert!( - last_problem_norm_fpr >= 0.0, + last_problem_norm_fpr >= T::zero(), "last_problem_norm_fpr should not be negative" ); self.last_problem_norm_fpr = last_problem_norm_fpr; self } - pub(crate) fn with_delta_y_norm(mut self, delta_y_norm: f64) -> Self { - assert!(delta_y_norm >= 0.0, "delta_y_norm must be nonnegative"); + pub(crate) fn with_delta_y_norm(mut self, delta_y_norm: T) -> Self { + assert!( + delta_y_norm >= T::zero(), + "delta_y_norm must be nonnegative" + ); self.delta_y_norm = delta_y_norm; self } - pub(crate) fn with_f2_norm(mut self, f2_norm: f64) -> Self { - assert!(f2_norm >= 0.0, "f2_norm must be nonnegative"); + pub(crate) fn with_f2_norm(mut self, f2_norm: T) -> Self { + assert!(f2_norm >= T::zero(), "f2_norm must be nonnegative"); self.f2_norm = f2_norm; self } - pub(crate) fn with_cost(mut self, cost: f64) -> Self { + pub(crate) fn with_cost(mut self, cost: T) -> Self { self.cost = cost; self } @@ -192,17 +200,17 @@ impl AlmOptimizerStatus { // ------------------------------------------------- /// Update cost (to be used when the cost needs to be scaled as a result of preconditioning) - pub fn update_cost(&mut self, new_cost: f64) { + pub fn update_cost(&mut self, new_cost: T) { self.cost = new_cost; } /// Update ALM infeasibility - pub fn update_f1_infeasibility(&mut self, new_alm_infeasibility: f64) { + pub fn update_f1_infeasibility(&mut self, new_alm_infeasibility: T) { self.delta_y_norm = new_alm_infeasibility; } /// Update PM infeasibility - pub fn update_f2_norm(&mut self, new_pm_infeasibility: f64) { + pub fn update_f2_norm(&mut self, new_pm_infeasibility: T) { self.f2_norm = new_pm_infeasibility; } @@ -249,7 +257,7 @@ impl AlmOptimizerStatus { /// /// Does not panic /// - pub fn lagrange_multipliers(&self) -> &Option> { + pub fn lagrange_multipliers(&self) -> &Option> { &self.lagrange_multipliers } @@ -259,7 +267,7 @@ impl AlmOptimizerStatus { /// /// Does not panic /// - pub fn last_problem_norm_fpr(&self) -> f64 { + pub fn last_problem_norm_fpr(&self) -> T { self.last_problem_norm_fpr } @@ -278,23 +286,23 @@ impl AlmOptimizerStatus { /// # Panics /// /// Does not panic - pub fn penalty(&self) -> f64 { + pub fn penalty(&self) -> T { self.penalty } /// Norm of Delta y divided by max{c, 1} - measure of infeasibility - pub fn delta_y_norm_over_c(&self) -> f64 { + pub fn delta_y_norm_over_c(&self) -> T { let c = self.penalty(); - self.delta_y_norm / if c < 1.0 { 1.0 } else { c } + self.delta_y_norm / if c < T::one() { T::one() } else { c } } /// Norm of F2(u) - measure of infeasibility of F2(u) = 0 - pub fn f2_norm(&self) -> f64 { + pub fn f2_norm(&self) -> T { self.f2_norm } /// Value of the cost function at the solution - pub fn cost(&self) -> f64 { + pub fn cost(&self) -> T { self.cost } } diff --git a/src/alm/alm_problem.rs b/src/alm/alm_problem.rs index 7827069e..aea3e9ec 100644 --- a/src/alm/alm_problem.rs +++ b/src/alm/alm_problem.rs @@ -1,4 +1,5 @@ -use crate::{constraints::Constraint, FunctionCallResult}; +use crate::constraints::Constraint; +use crate::{FunctionCallResult, OptFloat}; /// Definition of optimization problem to be solved with `AlmOptimizer`. The optimization /// problem has the general form @@ -32,16 +33,18 @@ pub struct AlmProblem< ConstraintsType, AlmSetC, LagrangeSetY, + T, > where // This is function F1: R^xn --> R^n1 (ALM) - MappingAlm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, + MappingAlm: Fn(&[T], &mut [T]) -> FunctionCallResult, // This is function F2: R^xn --> R^n2 (PM) - MappingPm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - ParametricGradientType: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, - ParametricCostType: Fn(&[f64], &[f64], &mut f64) -> FunctionCallResult, - ConstraintsType: Constraint, - AlmSetC: Constraint, - LagrangeSetY: Constraint, + MappingPm: Fn(&[T], &mut [T]) -> FunctionCallResult, + ParametricGradientType: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, + ParametricCostType: Fn(&[T], &[T], &mut T) -> FunctionCallResult, + ConstraintsType: Constraint, + AlmSetC: Constraint, + LagrangeSetY: Constraint, + T: OptFloat, { // // NOTE: the reason why we need to define different set types (ConstraintsType, @@ -67,6 +70,8 @@ pub struct AlmProblem< pub(crate) n1: usize, /// number of PM-type parameters (range dim of F2) pub(crate) n2: usize, + + _t: std::marker::PhantomData, } impl< @@ -77,6 +82,7 @@ impl< ConstraintsType, AlmSetC, LagrangeSetY, + T, > AlmProblem< MappingAlm, @@ -86,15 +92,17 @@ impl< ConstraintsType, AlmSetC, LagrangeSetY, + T, > where - MappingAlm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - MappingPm: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - ParametricGradientType: Fn(&[f64], &[f64], &mut [f64]) -> FunctionCallResult, - ParametricCostType: Fn(&[f64], &[f64], &mut f64) -> FunctionCallResult, - ConstraintsType: Constraint, - AlmSetC: Constraint, - LagrangeSetY: Constraint, + MappingAlm: Fn(&[T], &mut [T]) -> FunctionCallResult, + MappingPm: Fn(&[T], &mut [T]) -> FunctionCallResult, + ParametricGradientType: Fn(&[T], &[T], &mut [T]) -> FunctionCallResult, + ParametricCostType: Fn(&[T], &[T], &mut T) -> FunctionCallResult, + ConstraintsType: Constraint, + AlmSetC: Constraint, + LagrangeSetY: Constraint, + T: OptFloat, { ///Constructs new instance of `AlmProblem` /// @@ -163,6 +171,7 @@ where mapping_f2, n1, n2, + _t: std::marker::PhantomData, } } } diff --git a/src/alm/tests.rs b/src/alm/tests.rs index 0fd1c2e7..b2271b30 100644 --- a/src/alm/tests.rs +++ b/src/alm/tests.rs @@ -1,8 +1,8 @@ -use crate::{ - alm::*, - core::{constraints::*, panoc::*, ExitStatus}, - matrix_operations, mocks, FunctionCallResult, SolverError, -}; +use crate::alm::*; +use crate::core::constraints::*; +use crate::core::panoc::*; +use crate::core::ExitStatus; +use crate::{matrix_operations, mocks, FunctionCallResult, SolverError}; #[test] fn t_create_alm_cache() { diff --git a/src/constraints/affine_space.rs b/src/constraints/affine_space.rs index 0d17066e..a3bf79c9 100644 --- a/src/constraints/affine_space.rs +++ b/src/constraints/affine_space.rs @@ -1,4 +1,5 @@ use super::Constraint; +use crate::core::OptFloat; extern crate modcholesky; extern crate ndarray; @@ -13,16 +14,24 @@ type OpenVec = ArrayBase, Dim<[usize; 1]>>; /// An affine space here is defined as the set of solutions of a linear equation, $Ax = b$, /// that is, $E=\\{x\in\mathbb{R}^n: Ax = b\\}$, which is an affine space. It is assumed that /// the matrix $AA^\intercal$ is full-rank. -pub struct AffineSpace { - a_mat: OpenMat, - b_vec: OpenVec, - l: OpenMat, +pub struct AffineSpace +where + T: OptFloat, +{ + a_mat: OpenMat, + b_vec: OpenVec, + l: OpenMat, p: OpenVec, n_rows: usize, n_cols: usize, } -impl AffineSpace { +impl AffineSpace +where + T: OptFloat + 'static, + ndarray::Array2: + ModCholeskySE99, ndarray::Array1, ndarray::Array1>, +{ /// Construct a new affine space given the matrix $A\in\mathbb{R}^{m\times n}$ and /// the vector $b\in\mathbb{R}^m$ /// @@ -34,7 +43,7 @@ impl AffineSpace { /// ## Returns /// New Affine Space structure /// - pub fn new(a: Vec, b: Vec) -> Self { + pub fn new(a: Vec, b: Vec) -> Self { // Infer dimensions of A and b let n_rows = b.len(); let n_elements_a = a.len(); @@ -67,7 +76,10 @@ impl AffineSpace { } } -impl Constraint for AffineSpace { +impl Constraint for AffineSpace +where + T: OptFloat + 'static, +{ /// Projection onto the set $E = \\{x: Ax = b\\}$, which is computed by /// $$P_E(x) = x - A^\intercal z(x),$$ /// where $z$ is the solution of the linear system @@ -99,7 +111,7 @@ impl Constraint for AffineSpace { /// ``` /// /// The result is stored in `x` and it can be verified that $Ax = b$. - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { let m = self.n_rows; let n = self.n_cols; let chol = &self.l; @@ -114,21 +126,23 @@ impl Constraint for AffineSpace { // Step 1: Solve Ly = b(P) // TODO: Make `y` into an attribute; however, to do this, we need to change // &self to &mut self, which will require a mild refactoring - let mut y = vec![0.; m]; + let mut y = vec![T::zero(); m]; for i in 0..m { y[i] = err[perm[i]]; for j in 0..i { - y[i] -= chol[(i, j)] * y[j]; + let val = y[j]; + y[i] -= chol[(i, j)] * val; } y[i] /= chol[(i, i)]; } // Step 2: Solve L'z(P) = y - let mut z = vec![0.; m]; + let mut z = vec![T::zero(); m]; for i in 1..=m { z[perm[m - i]] = y[m - i]; for j in 1..i { - z[perm[m - i]] -= chol[(m - j, m - i)] * z[perm[m - j]]; + let val = z[perm[m - j]]; + z[perm[m - i]] -= chol[(m - j, m - i)] * val; } z[perm[m - i]] /= chol[(m - i, m - i)]; } @@ -138,7 +152,7 @@ impl Constraint for AffineSpace { let w = self.a_mat.t().dot(&z_arr); // Step 4: x <-- x - A'(AA')\(Ax-b) - x.iter_mut().zip(w.iter()).for_each(|(xi, wi)| *xi -= wi); + x.iter_mut().zip(w.iter()).for_each(|(xi, wi)| *xi -= *wi); } /// Affine sets are convex. diff --git a/src/constraints/ball1.rs b/src/constraints/ball1.rs index 9d30a75e..9f0d2fc1 100644 --- a/src/constraints/ball1.rs +++ b/src/constraints/ball1.rs @@ -1,20 +1,26 @@ -use super::Constraint; -use super::Simplex; +use super::{Constraint, Simplex}; +use crate::core::OptFloat; #[derive(Copy, Clone)] /// A norm-1 ball, that is, a set given by $B_1^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert_1 \leq r\\}$ /// or a ball-1 centered at a point $x_c$, that is, $B_1^{x_c, r} = \\{x \in \mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert_1 \leq r\\}$ -pub struct Ball1<'a> { - center: Option<&'a [f64]>, - radius: f64, - simplex: Simplex, +pub struct Ball1<'a, T> +where + T: OptFloat, +{ + center: Option<&'a [T]>, + radius: T, + simplex: Simplex, } -impl<'a> Ball1<'a> { +impl<'a, T> Ball1<'a, T> +where + T: OptFloat + PartialOrd + Copy, +{ /// Construct a new ball-1 with given center and radius. /// If no `center` is given, then it is assumed to be in the origin - pub fn new(center: Option<&'a [f64]>, radius: f64) -> Self { - assert!(radius > 0.0); + pub fn new(center: Option<&'a [T]>, radius: T) -> Self { + assert!(radius > T::zero()); let simplex = Simplex::new(radius); Ball1 { center, @@ -23,24 +29,27 @@ impl<'a> Ball1<'a> { } } - fn project_on_ball1_centered_at_origin(&self, x: &mut [f64]) { + fn project_on_ball1_centered_at_origin(&self, x: &mut [T]) { if crate::matrix_operations::norm1(x) > self.radius { // u = |x| (copied) - let mut u = vec![0.0; x.len()]; + let mut u = vec![T::zero(); x.len()]; u.iter_mut() .zip(x.iter()) - .for_each(|(ui, &xi)| *ui = f64::abs(xi)); + .for_each(|(ui, &xi)| *ui = xi.abs()); // u = P_simplex(u) self.simplex.project(&mut u); x.iter_mut() .zip(u.iter()) - .for_each(|(xi, &ui)| *xi = f64::signum(*xi) * ui); + .for_each(|(xi, &ui)| *xi = xi.signum() * ui); } } } -impl<'a> Constraint for Ball1<'a> { - fn project(&self, x: &mut [f64]) { +impl<'a, T> Constraint for Ball1<'a, T> +where + T: OptFloat + PartialOrd + Copy, +{ + fn project(&self, x: &mut [T]) { if let Some(center) = &self.center { x.iter_mut() .zip(center.iter()) diff --git a/src/constraints/ball2.rs b/src/constraints/ball2.rs index a26d2944..0990b65d 100644 --- a/src/constraints/ball2.rs +++ b/src/constraints/ball2.rs @@ -1,30 +1,40 @@ use super::Constraint; +use crate::core::OptFloat; #[derive(Copy, Clone)] /// A Euclidean ball, that is, a set given by $B_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert \leq r\\}$ /// or a Euclidean ball centered at a point $x_c$, that is, $B_2^{x_c, r} = \\{x \in \mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert \leq r\\}$ -pub struct Ball2<'a> { - center: Option<&'a [f64]>, - radius: f64, +pub struct Ball2<'a, T> +where + T: OptFloat, +{ + center: Option<&'a [T]>, + radius: T, } -impl<'a> Ball2<'a> { +impl<'a, T> Ball2<'a, T> +where + T: OptFloat, +{ /// Construct a new Euclidean ball with given center and radius /// If no `center` is given, then it is assumed to be in the origin - pub fn new(center: Option<&'a [f64]>, radius: f64) -> Self { - assert!(radius > 0.0); + pub fn new(center: Option<&'a [T]>, radius: T) -> Self { + assert!(radius > T::zero()); Ball2 { center, radius } } } -impl<'a> Constraint for Ball2<'a> { - fn project(&self, x: &mut [f64]) { +impl<'a, T> Constraint for Ball2<'a, T> +where + T: OptFloat, +{ + fn project(&self, x: &mut [T]) { if let Some(center) = &self.center { - let mut norm_difference = 0.0; + let mut norm_difference = T::zero(); x.iter().zip(center.iter()).for_each(|(a, b)| { let diff_ = *a - *b; - norm_difference += diff_ * diff_ + norm_difference += diff_ * diff_; }); norm_difference = norm_difference.sqrt(); @@ -38,7 +48,7 @@ impl<'a> Constraint for Ball2<'a> { let norm_x = crate::matrix_operations::norm2(x); if norm_x > self.radius { let norm_over_radius = norm_x / self.radius; - x.iter_mut().for_each(|x_| *x_ /= norm_over_radius); + x.iter_mut().for_each(|x_| *x_ = *x_ / norm_over_radius); } } } diff --git a/src/constraints/ballinf.rs b/src/constraints/ballinf.rs index fcaeb800..5d7e5afc 100644 --- a/src/constraints/ballinf.rs +++ b/src/constraints/ballinf.rs @@ -1,26 +1,36 @@ use super::Constraint; +use crate::core::OptFloat; #[derive(Copy, Clone)] /// An infinity ball defined as $B_\infty^r = \\{x\in\mathbb{R}^n {}:{} \Vert{}x{}\Vert_{\infty} \leq r\\}$, /// where $\Vert{}\cdot{}\Vert_{\infty}$ is the infinity norm. The infinity ball centered at a point /// $x_c$ is defined as $B_\infty^{x_c,r} = \\{x\in\mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert_{\infty} \leq r\\}$. /// -pub struct BallInf<'a> { - center: Option<&'a [f64]>, - radius: f64, +pub struct BallInf<'a, T> +where + T: OptFloat, +{ + center: Option<&'a [T]>, + radius: T, } -impl<'a> BallInf<'a> { +impl<'a, T> BallInf<'a, T> +where + T: OptFloat, +{ /// Construct a new infinity-norm ball with given center and radius /// If no `center` is given, then it is assumed to be in the origin /// - pub fn new(center: Option<&'a [f64]>, radius: f64) -> Self { - assert!(radius > 0.0); + pub fn new(center: Option<&'a [T]>, radius: T) -> Self { + assert!(radius > T::zero()); BallInf { center, radius } } } -impl<'a> Constraint for BallInf<'a> { +impl<'a, T> Constraint for BallInf<'a, T> +where + T: OptFloat, +{ /// Computes the projection of a given vector `x` on the current infinity ball. /// /// @@ -42,12 +52,12 @@ impl<'a> Constraint for BallInf<'a> { /// /// for all $i=1,\ldots, n$. /// - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { if let Some(center) = &self.center { x.iter_mut() .zip(center.iter()) .filter(|(&mut xi, &ci)| (xi - ci).abs() > self.radius) - .for_each(|(xi, ci)| *xi = ci + (*xi - ci).signum() * self.radius); + .for_each(|(xi, ci)| *xi = *ci + (*xi - *ci).signum() * self.radius); } else { x.iter_mut() .filter(|xi| xi.abs() > self.radius) diff --git a/src/constraints/cartesian_product.rs b/src/constraints/cartesian_product.rs index c298ff06..69dc7e69 100644 --- a/src/constraints/cartesian_product.rs +++ b/src/constraints/cartesian_product.rs @@ -1,5 +1,5 @@ use super::Constraint; - +use crate::core::OptFloat; /// Cartesian product of constraints /// /// Cartesian product of constraints, $C_0, C_1, \ldots, C_{n-1}$, @@ -22,12 +22,18 @@ use super::Constraint; /// for all $i=0,\ldots, n-1$. /// #[derive(Default)] -pub struct CartesianProduct<'a> { +pub struct CartesianProduct<'a, T> +where + T: OptFloat, +{ idx: Vec, - constraints: Vec>, + constraints: Vec + 'a>>, } -impl<'a> CartesianProduct<'a> { +impl<'a, T> CartesianProduct<'a, T> +where + T: OptFloat, +{ /// Construct new instance of Cartesian product of constraints /// /// # Note @@ -116,7 +122,7 @@ impl<'a> CartesianProduct<'a> { /// ``` /// The method will panic if any of the associated projections panics. /// - pub fn add_constraint(mut self, ni: usize, constraint: impl Constraint + 'a) -> Self { + pub fn add_constraint(mut self, ni: usize, constraint: impl Constraint + 'a) -> Self { assert!( self.dimension() < ni, "provided index is smaller than or equal to previous index, or zero" @@ -127,7 +133,10 @@ impl<'a> CartesianProduct<'a> { } } -impl<'a> Constraint for CartesianProduct<'a> { +impl<'a, T> Constraint for CartesianProduct<'a, T> +where + T: OptFloat, +{ /// Project onto Cartesian product of constraints /// /// The given vector `x` is updated with the projection on the set @@ -136,7 +145,7 @@ impl<'a> Constraint for CartesianProduct<'a> { /// /// The method will panic if the dimension of `x` is not equal to the /// dimension of the Cartesian product (see `dimension()`) - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { assert!(x.len() == self.dimension(), "x has wrong size"); let mut j = 0; self.idx diff --git a/src/constraints/epigraph_squared_norm.rs b/src/constraints/epigraph_squared_norm.rs index 560b34bd..48946ec2 100644 --- a/src/constraints/epigraph_squared_norm.rs +++ b/src/constraints/epigraph_squared_norm.rs @@ -1,22 +1,35 @@ -use crate::matrix_operations; - use super::Constraint; +use crate::core::OptFloat; +use crate::matrix_operations; #[derive(Copy, Clone, Default)] /// The epigraph of the squared Eucliden norm is a set of the form /// $X = \\{x = (z, t) \in \mathbb{R}^{n}\times \mathbb{R} {}:{} \\|z\\|^2 \leq t \\}.$ -pub struct EpigraphSquaredNorm {} +pub struct EpigraphSquaredNorm +where + T: OptFloat, +{ + _phantom: std::marker::PhantomData, +} -impl EpigraphSquaredNorm { +impl EpigraphSquaredNorm +where + T: OptFloat, +{ /// Create a new instance of the epigraph of the squared norm. /// /// Note that you do not need to specify the dimension. pub fn new() -> Self { - EpigraphSquaredNorm {} + EpigraphSquaredNorm { + _phantom: std::marker::PhantomData, + } } } -impl Constraint for EpigraphSquaredNorm { +impl Constraint for EpigraphSquaredNorm +where + T: OptFloat + num::FromPrimitive + roots::FloatType, +{ ///Project on the epigraph of the squared Euclidean norm. /// /// The projection is computed as detailed @@ -34,31 +47,32 @@ impl Constraint for EpigraphSquaredNorm { /// let mut x = [1., 2., 3., 4.]; /// epi.project(&mut x); /// ``` - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { let nx = x.len() - 1; assert!(nx > 0, "x must have a length of at least 2"); - let z: &[f64] = &x[..nx]; - let t: f64 = x[nx]; + let z: &[T] = &x[..nx]; + let t: T = x[nx]; let norm_z_sq = matrix_operations::norm2_squared(z); if norm_z_sq <= t { return; } - let theta = 1. - 2. * t; - let a3 = 4.; - let a2 = 4. * theta; + let theta = ::one() - T::from_f64(2.0).unwrap() * t; + let a3 = T::from_f64(4.0).unwrap(); + let a2 = T::from_f64(4.0).unwrap() * theta; let a1 = theta * theta; let a0 = -norm_z_sq; let cubic_poly_roots = roots::find_roots_cubic(a3, a2, a1, a0); - let mut right_root = f64::NAN; - let mut scaling = f64::NAN; + let mut right_root = T::nan(); + let mut scaling = T::nan(); // Find right root cubic_poly_roots.as_ref().iter().for_each(|ri| { - if *ri > 0. { - let denom = 1. + 2. * (*ri - t); - if ((norm_z_sq / (denom * denom)) - *ri).abs() < 1e-6 { + if *ri > ::zero() { + let denom = ::one() + T::from_f64(2.0).unwrap() * (*ri - t); + let tolerance = T::from_f64(1e-6).unwrap(); + if num::Float::abs((norm_z_sq / (denom * denom)) - *ri) < tolerance { right_root = *ri; scaling = denom; } @@ -66,25 +80,27 @@ impl Constraint for EpigraphSquaredNorm { }); // Refinement of root with Newton-Raphson - let mut refinement_error = 1.; + let mut refinement_error = ::one(); let newton_max_iters: usize = 5; - let newton_eps = 1e-14; + let newton_eps = T::from_f64(1e-14).unwrap(); let mut zsol = right_root; let mut iter = 0; while refinement_error > newton_eps && iter < newton_max_iters { let zsol_sq = zsol * zsol; let zsol_cb = zsol_sq * zsol; let p_z = a3 * zsol_cb + a2 * zsol_sq + a1 * zsol + a0; - let dp_z = 3. * a3 * zsol_sq + 2. * a2 * zsol + a1; + let dp_z = T::from_f64(3.0).unwrap() * a3 * zsol_sq + + T::from_f64(2.0).unwrap() * a2 * zsol + + a1; zsol -= p_z / dp_z; - refinement_error = p_z.abs(); + refinement_error = num::Float::abs(p_z); iter += 1; } right_root = zsol; // Projection for xi in x.iter_mut().take(nx) { - *xi /= scaling; + *xi = *xi / scaling; } x[nx] = right_root; } diff --git a/src/constraints/finite.rs b/src/constraints/finite.rs index ca36f431..fc63ae13 100644 --- a/src/constraints/finite.rs +++ b/src/constraints/finite.rs @@ -1,17 +1,23 @@ use super::Constraint; - +use crate::core::OptFloat; /// /// A finite set, $X = \\{x_1, x_2, \ldots, x_n\\}\subseteq\mathbb{R}^n$, given vectors /// $x_i\in\mathbb{R}^n$ /// #[derive(Clone, Copy)] -pub struct FiniteSet<'a> { +pub struct FiniteSet<'a, T> +where + T: OptFloat, +{ /// The data is stored in a Vec-of-Vec datatype, that is, a vector /// of vectors - data: &'a [&'a [f64]], + data: &'a [&'a [T]], } -impl<'a> FiniteSet<'a> { +impl<'a, T> FiniteSet<'a, T> +where + T: OptFloat, +{ /// Construct a finite set, $X = \\{x_1, x_2, \ldots, x_n\\}$, given vectors /// $x_i\in\mathbb{R}^n$ /// @@ -41,7 +47,7 @@ impl<'a> FiniteSet<'a> { /// This method will panic if (i) the given vector of data is empty /// and (ii) if the given vectors have unequal dimensions. /// - pub fn new(data: &'a [&'a [f64]]) -> Self { + pub fn new(data: &'a [&'a [T]]) -> Self { // Do a sanity check... assert!(!data.is_empty(), "empty data not allowed"); let n = data[0].len(); @@ -52,7 +58,10 @@ impl<'a> FiniteSet<'a> { } } -impl<'a> Constraint for FiniteSet<'a> { +impl<'a, T> Constraint for FiniteSet<'a, T> +where + T: OptFloat + std::ops::AddAssign, +{ /// /// Projection on the current finite set /// @@ -84,9 +93,9 @@ impl<'a> Constraint for FiniteSet<'a> { /// /// Does not panic /// - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { let mut idx: usize = 0; - let mut best_distance: f64 = num::Float::infinity(); + let mut best_distance: T = T::infinity(); for (i, v) in self.data.iter().enumerate() { let dist = crate::matrix_operations::norm2_squared_diff(v, x); if dist < best_distance { diff --git a/src/constraints/halfspace.rs b/src/constraints/halfspace.rs index 5d9d4a38..495821d6 100644 --- a/src/constraints/halfspace.rs +++ b/src/constraints/halfspace.rs @@ -1,18 +1,25 @@ use super::Constraint; +use crate::core::OptFloat; use crate::matrix_operations; #[derive(Clone)] /// A halfspace is a set given by $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle \leq b\\}$. -pub struct Halfspace<'a> { +pub struct Halfspace<'a, T> +where + T: OptFloat, +{ /// normal vector - normal_vector: &'a [f64], + normal_vector: &'a [T], /// offset - offset: f64, + offset: T, /// squared Euclidean norm of the normal vector (computed once upon construction) - normal_vector_squared_norm: f64, + normal_vector_squared_norm: T, } -impl<'a> Halfspace<'a> { +impl<'a, T> Halfspace<'a, T> +where + T: OptFloat, +{ /// A halfspace is a set given by $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle \leq b\\}$, /// where $c$ is the normal vector of the halfspace and $b$ is an offset. /// @@ -45,7 +52,7 @@ impl<'a> Halfspace<'a> { /// halfspace.project(&mut x); /// ``` /// - pub fn new(normal_vector: &'a [f64], offset: f64) -> Self { + pub fn new(normal_vector: &'a [T], offset: T) -> Self { let normal_vector_squared_norm = matrix_operations::norm2_squared(normal_vector); Halfspace { normal_vector, @@ -55,7 +62,10 @@ impl<'a> Halfspace<'a> { } } -impl<'a> Constraint for Halfspace<'a> { +impl<'a, T> Constraint for Halfspace<'a, T> +where + T: OptFloat + std::ops::SubAssign, +{ /// Projects on halfspace using the following formula: /// /// $$\begin{aligned} @@ -79,13 +89,13 @@ impl<'a> Constraint for Halfspace<'a> { /// This method panics if the length of `x` is not equal to the dimension /// of the halfspace. /// - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { let inner_product = matrix_operations::inner_product(x, self.normal_vector); if inner_product > self.offset { let factor = (inner_product - self.offset) / self.normal_vector_squared_norm; x.iter_mut() .zip(self.normal_vector.iter()) - .for_each(|(x, normal_vector_i)| *x -= factor * normal_vector_i); + .for_each(|(x, normal_vector_i)| *x -= factor * *normal_vector_i); } } diff --git a/src/constraints/hyperplane.rs b/src/constraints/hyperplane.rs index d3e6dd46..f77bdf21 100644 --- a/src/constraints/hyperplane.rs +++ b/src/constraints/hyperplane.rs @@ -1,18 +1,25 @@ use super::Constraint; +use crate::core::OptFloat; use crate::matrix_operations; #[derive(Clone)] /// A hyperplane is a set given by $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle = b\\}$. -pub struct Hyperplane<'a> { +pub struct Hyperplane<'a, T> +where + T: OptFloat, +{ /// normal vector - normal_vector: &'a [f64], + normal_vector: &'a [T], /// offset - offset: f64, + offset: T, /// squared Euclidean norm of the normal vector (computed once upon construction) - normal_vector_squared_norm: f64, + normal_vector_squared_norm: T, } -impl<'a> Hyperplane<'a> { +impl<'a, T> Hyperplane<'a, T> +where + T: OptFloat, +{ /// A hyperplane is a set given by $H = \\{x \in \mathbb{R}^n {}:{} \langle c, x\rangle = b\\}$, /// where $c$ is the normal vector of the hyperplane and $b$ is an offset. /// @@ -45,7 +52,7 @@ impl<'a> Hyperplane<'a> { /// hyperplane.project(&mut x); /// ``` /// - pub fn new(normal_vector: &'a [f64], offset: f64) -> Self { + pub fn new(normal_vector: &'a [T], offset: T) -> Self { let normal_vector_squared_norm = matrix_operations::norm2_squared(normal_vector); Hyperplane { normal_vector, @@ -55,7 +62,10 @@ impl<'a> Hyperplane<'a> { } } -impl<'a> Constraint for Hyperplane<'a> { +impl<'a, T> Constraint for Hyperplane<'a, T> +where + T: OptFloat + std::ops::SubAssign, +{ /// Projects on the hyperplane using the formula: /// /// $$\begin{aligned} @@ -76,12 +86,12 @@ impl<'a> Constraint for Hyperplane<'a> { /// This method panics if the length of `x` is not equal to the dimension /// of the hyperplane. /// - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { let inner_product = matrix_operations::inner_product(x, self.normal_vector); let factor = (inner_product - self.offset) / self.normal_vector_squared_norm; x.iter_mut() .zip(self.normal_vector.iter()) - .for_each(|(x, nrm_vct)| *x -= factor * nrm_vct); + .for_each(|(x, nrm_vct)| *x -= factor * *nrm_vct); } /// Hyperplanes are convex sets diff --git a/src/constraints/mod.rs b/src/constraints/mod.rs index dd0b7038..54270bad 100644 --- a/src/constraints/mod.rs +++ b/src/constraints/mod.rs @@ -40,11 +40,16 @@ pub use soc::SecondOrderCone; pub use sphere2::Sphere2; pub use zero::Zero; +use crate::core::OptFloat; + /// A set which can be used as a constraint /// /// This trait defines an abstract function that allows to compute projections /// on sets; this is implemented by a series of structures (see below for details) -pub trait Constraint { +pub trait Constraint +where + T: OptFloat, +{ /// Projection onto the set, that is, /// /// $$ @@ -55,7 +60,7 @@ pub trait Constraint { /// /// - `x`: The given vector $x$ is updated with the projection on the set /// - fn project(&self, x: &mut [f64]); + fn project(&self, x: &mut [T]); /// Returns true if and only if the set is convex fn is_convex(&self) -> bool; diff --git a/src/constraints/no_constraints.rs b/src/constraints/no_constraints.rs index 88df8845..2b160a4a 100644 --- a/src/constraints/no_constraints.rs +++ b/src/constraints/no_constraints.rs @@ -1,5 +1,5 @@ use super::Constraint; - +use crate::core::OptFloat; /// The whole space, no constraints #[derive(Default, Clone, Copy)] pub struct NoConstraints {} @@ -12,8 +12,11 @@ impl NoConstraints { } } -impl Constraint for NoConstraints { - fn project(&self, _x: &mut [f64]) {} +impl Constraint for NoConstraints +where + T: OptFloat, +{ + fn project(&self, _x: &mut [T]) {} fn is_convex(&self) -> bool { true diff --git a/src/constraints/rectangle.rs b/src/constraints/rectangle.rs index f8c8cf9f..d1d37a87 100644 --- a/src/constraints/rectangle.rs +++ b/src/constraints/rectangle.rs @@ -1,5 +1,5 @@ use super::Constraint; - +use crate::core::OptFloat; #[derive(Clone, Copy)] /// /// A rectangle, $R = \\{x \in \mathbb{R}^n {}:{} x_{\min} {}\leq{} x {}\leq{} x_{\max}\\}$ @@ -7,12 +7,18 @@ use super::Constraint; /// A set of the form $\\{x \in \mathbb{R}^n {}:{} x_{\min} {}\leq{} x {}\leq{} x_{\max}\\}$, /// where $\leq$ is meant in the element-wise sense and either of $x_{\min}$ and $x_{\max}$ can /// be equal to infinity. -pub struct Rectangle<'a> { - xmin: Option<&'a [f64]>, - xmax: Option<&'a [f64]>, +pub struct Rectangle<'a, T> +where + T: OptFloat, +{ + xmin: Option<&'a [T]>, + xmax: Option<&'a [T]>, } -impl<'a> Rectangle<'a> { +impl<'a, T> Rectangle<'a, T> +where + T: OptFloat, +{ /// Construct a new rectangle with given $x_{\min}$ and $x_{\max}$ /// /// # Arguments @@ -34,7 +40,7 @@ impl<'a> Rectangle<'a> { /// - Both `xmin` and `xmax` have been provided, but they have incompatible /// dimensions /// - pub fn new(xmin: Option<&'a [f64]>, xmax: Option<&'a [f64]>) -> Self { + pub fn new(xmin: Option<&'a [T]>, xmax: Option<&'a [T]>) -> Self { assert!(xmin.is_some() || xmax.is_some()); // xmin or xmax must be Some assert!( xmin.is_none() || xmax.is_none() || xmin.unwrap().len() == xmax.unwrap().len(), @@ -44,8 +50,11 @@ impl<'a> Rectangle<'a> { } } -impl<'a> Constraint for Rectangle<'a> { - fn project(&self, x: &mut [f64]) { +impl<'a, T> Constraint for Rectangle<'a, T> +where + T: OptFloat, +{ + fn project(&self, x: &mut [T]) { if let Some(xmin) = &self.xmin { x.iter_mut().zip(xmin.iter()).for_each(|(x_, xmin_)| { if *x_ < *xmin_ { diff --git a/src/constraints/simplex.rs b/src/constraints/simplex.rs index 65503f1a..5535721e 100644 --- a/src/constraints/simplex.rs +++ b/src/constraints/simplex.rs @@ -1,49 +1,59 @@ use super::Constraint; - +use crate::core::OptFloat; #[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 { +pub struct Simplex +where + T: OptFloat, +{ /// Simplex level - alpha: f64, + alpha: T, } -impl Simplex { +impl Simplex +where + T: OptFloat, +{ /// 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"); + pub fn new(alpha: T) -> Self { + assert!(alpha > T::zero(), "alpha is nonpositive"); Simplex { alpha } } } -impl Constraint for Simplex { +impl Constraint for Simplex +where + T: OptFloat + PartialOrd + Copy + num::FromPrimitive, +{ /// Project onto $\Delta_\alpha^n$ using Condat's fast projection algorithm. /// /// See: Laurent Condat. Fast Projection onto the Simplex and the $\ell_1$ Ball. /// Mathematical Programming, Series A, Springer, 2016, 158 (1), pp.575-585. /// ⟨10.1007/s10107-015-0946-6⟩. - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { let a = &self.alpha; // ---- step 1 - let mut v = Vec::::with_capacity(x.len()); // vector containing x[0] + let mut v = Vec::::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 = Vec::new(); // empty vector of f64 - let mut rho: f64 = x[0] - a; // 64 bit float + let mut v_tilde: Vec = Vec::new(); // empty vector of T + let mut rho: T = x[0] - *a; // T 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 { + let len_plus_one = T::from(v.len() + 1).unwrap(); + rho = rho + (*x_n - rho) / len_plus_one; + if rho > *x_n - *a { v.push(*x_n); } else { v_tilde.extend(&v); v = vec![*x_n]; - rho = *x_n - a; + rho = *x_n - *a; } } }); @@ -53,7 +63,8 @@ impl Constraint for Simplex { 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); + let len_t = T::from(v.len()).unwrap(); + rho = rho + (*v_t_n - rho) / len_t; } }); } @@ -67,7 +78,8 @@ impl Constraint for Simplex { if *v_n <= rho { hit_list.push(n); current_len_v -= 1; - rho += (rho - *v_n) / (current_len_v as f64); + let current_len_t = T::from(current_len_v).unwrap(); + rho = rho + (rho - *v_n) / current_len_t; } }); hit_list.iter().rev().for_each(|target| { @@ -79,7 +91,7 @@ impl Constraint for Simplex { } // ---- step 6 - let zero: f64 = 0.0; + let zero: T = T::zero(); x.iter_mut().for_each(|x_n| *x_n = zero.max(*x_n - rho)); } diff --git a/src/constraints/soc.rs b/src/constraints/soc.rs index 1d6c969a..c03cee78 100644 --- a/src/constraints/soc.rs +++ b/src/constraints/soc.rs @@ -1,4 +1,5 @@ use super::Constraint; +use crate::core::OptFloat; use crate::matrix_operations; #[derive(Clone, Copy)] @@ -17,11 +18,17 @@ use crate::matrix_operations; /// 1996 doctoral dissertation: Projection Algorithms and Monotone Operators /// (p. 40, Theorem 3.3.6). /// -pub struct SecondOrderCone { - alpha: f64, +pub struct SecondOrderCone +where + T: OptFloat, +{ + alpha: T, } -impl SecondOrderCone { +impl SecondOrderCone +where + T: OptFloat, +{ /// Construct a new instance of SecondOrderCone with parameter `alpha` /// /// A second-order cone with parameter alpha is the set @@ -39,13 +46,16 @@ impl SecondOrderCone { /// # Panics /// /// The method panics if the given parameter `alpha` is nonpositive. - pub fn new(alpha: f64) -> SecondOrderCone { - assert!(alpha > 0.0); // alpha must be positive + pub fn new(alpha: T) -> SecondOrderCone { + assert!(alpha > T::zero()); // alpha must be positive SecondOrderCone { alpha } } } -impl Constraint for SecondOrderCone { +impl Constraint for SecondOrderCone +where + T: OptFloat, +{ /// Project on the second-order cone (updates the given vector/slice) /// /// # Arguments @@ -57,7 +67,7 @@ impl Constraint for SecondOrderCone { /// /// The methods panics is the length of `x` is less than 2. /// - fn project(&self, x: &mut [f64]) { + fn project(&self, x: &mut [T]) { // x = (z, r) let n = x.len(); assert!(n >= 2, "x must be of dimension at least 2"); @@ -65,12 +75,12 @@ impl Constraint for SecondOrderCone { let r = x[n - 1]; let norm_z = matrix_operations::norm2(z); if self.alpha * norm_z <= -r { - x.iter_mut().for_each(|v| *v = 0.0); + x.iter_mut().for_each(|v| *v = T::zero()); } else if norm_z > self.alpha * r { - let beta = (self.alpha * norm_z + r) / (self.alpha.powi(2) + 1.0); + let beta = (self.alpha * norm_z + r) / (self.alpha.powi(2) + T::one()); x[..n - 1] .iter_mut() - .for_each(|v| *v *= self.alpha * beta / norm_z); + .for_each(|v| *v = *v * self.alpha * beta / norm_z); x[n - 1] = beta; } } diff --git a/src/constraints/sphere2.rs b/src/constraints/sphere2.rs index 8ea4d344..31d5ffc6 100644 --- a/src/constraints/sphere2.rs +++ b/src/constraints/sphere2.rs @@ -1,23 +1,32 @@ use super::Constraint; - +use crate::core::OptFloat; #[derive(Copy, Clone)] /// A Euclidean sphere, that is, a set given by $S_2^r = \\{x \in \mathbb{R}^n {}:{} \Vert{}x{}\Vert = r\\}$ /// or a Euclidean sphere centered at a point $x_c$, that is, $S_2^{x_c, r} = \\{x \in \mathbb{R}^n {}:{} \Vert{}x-x_c{}\Vert = r\\}$ -pub struct Sphere2<'a> { - center: Option<&'a [f64]>, - radius: f64, +pub struct Sphere2<'a, T> +where + T: OptFloat, +{ + center: Option<&'a [T]>, + radius: T, } -impl<'a> Sphere2<'a> { +impl<'a, T> Sphere2<'a, T> +where + T: OptFloat, +{ /// Construct a new Euclidean sphere with given center and radius /// If no `center` is given, then it is assumed to be in the origin - pub fn new(center: Option<&'a [f64]>, radius: f64) -> Self { - assert!(radius > 0.0); + pub fn new(center: Option<&'a [T]>, radius: T) -> Self { + assert!(radius > T::zero()); Sphere2 { center, radius } } } -impl<'a> Constraint for Sphere2<'a> { +impl<'a, T> Constraint for Sphere2<'a, T> +where + T: OptFloat, +{ /// Projection onto the sphere, $S_{r, c}$ with radius $r$ and center $c$. /// If $x\neq c$, the projection is uniquely defined by /// @@ -33,8 +42,8 @@ impl<'a> Constraint for Sphere2<'a> { /// /// - `x`: The given vector $x$ is updated with the projection on the set /// - fn project(&self, x: &mut [f64]) { - let epsilon = 1e-12; + fn project(&self, x: &mut [T]) { + let epsilon = T::from(1e-12).unwrap(); if let Some(center) = &self.center { let norm_difference = crate::matrix_operations::norm2_squared_diff(x, center).sqrt(); if norm_difference <= epsilon { @@ -52,7 +61,7 @@ impl<'a> Constraint for Sphere2<'a> { return; } let norm_over_radius = self.radius / norm_x; - x.iter_mut().for_each(|x_| *x_ *= norm_over_radius); + x.iter_mut().for_each(|x_| *x_ = *x_ * norm_over_radius); } } diff --git a/src/constraints/tests.rs b/src/constraints/tests.rs index 26343688..ca86d7f8 100644 --- a/src/constraints/tests.rs +++ b/src/constraints/tests.rs @@ -1,7 +1,6 @@ -use crate::matrix_operations; +use rand; use super::*; -use rand; #[test] fn t_zero_set() { @@ -268,7 +267,7 @@ fn t_ball2_at_origin_different_radius_inside() { #[test] fn t_ball2_at_center_different_radius_outside() { - let radius = 1.2; + let radius = 1.2_f64; let mut x = [1.0, 1.0]; let center = [-0.8, -1.1]; let ball = Ball2::new(Some(¢er), radius); @@ -430,7 +429,7 @@ fn t_second_order_cone_case_ii() { #[test] fn t_second_order_cone_case_iii() { - let alpha = 1.5; + let alpha = 1.5_f64; let soc = SecondOrderCone::new(alpha); let mut x = vec![1.0, 1.0, 0.1]; soc.project(&mut x); @@ -615,7 +614,7 @@ fn t_is_convex_soc() { #[test] fn t_is_convex_zero() { let zero = Zero::new(); - assert!(zero.is_convex()); + assert!(>::is_convex(&zero)); } #[test] @@ -897,7 +896,7 @@ fn t_epigraph_squared_norm() { let t = 0.01 * i as f64; let mut x = [1., 2., 3., t]; epi.project(&mut x); - let err = (matrix_operations::norm2_squared(&x[..3]) - x[3]).abs(); + let err = (crate::matrix_operations::norm2_squared(&x[..3]) - x[3]).abs(); assert!(err < 1e-10, "wrong projection on epigraph of squared norm"); } } @@ -987,5 +986,5 @@ fn t_affine_space_single_row() { fn t_affine_space_wrong_dimensions() { let a = vec![0.5, 0.1, 0.2, -0.3, -0.6, 0.3, 0., 0.5, 1.0, 0.1, -1.0]; let b = vec![1., 2., -0.5]; - let _ = AffineSpace::new(a, b); + let _: AffineSpace = AffineSpace::new(a, b); } diff --git a/src/constraints/zero.rs b/src/constraints/zero.rs index 81064d33..5261e2fa 100644 --- a/src/constraints/zero.rs +++ b/src/constraints/zero.rs @@ -1,5 +1,5 @@ use super::Constraint; - +use crate::core::OptFloat; #[derive(Clone, Copy, Default)] /// Set Zero, $\\{0\\}$ pub struct Zero {} @@ -11,11 +11,14 @@ impl Zero { } } -impl Constraint for Zero { +impl Constraint for Zero +where + T: OptFloat, +{ /// Computes the projection on $\\{0\\}$, that is, $\Pi_{\\{0\\}}(x) = 0$ /// for all $x$ - fn project(&self, x: &mut [f64]) { - x.iter_mut().for_each(|xi| *xi = 0.0); + fn project(&self, x: &mut [T]) { + x.iter_mut().for_each(|xi| *xi = T::zero()); } fn is_convex(&self) -> bool { diff --git a/src/core/fbs/fbs_cache.rs b/src/core/fbs/fbs_cache.rs index de134799..eed4ceee 100644 --- a/src/core/fbs/fbs_cache.rs +++ b/src/core/fbs/fbs_cache.rs @@ -2,18 +2,26 @@ //! use std::num::NonZeroUsize; +use crate::core::OptFloat; + /// Cache for the forward-backward splitting (FBS), or projected gradient, algorithm /// /// This struct allocates memory needed for the FBS algorithm -pub struct FBSCache { - pub(crate) work_gradient_u: Vec, - pub(crate) work_u_previous: Vec, - pub(crate) gamma: f64, - pub(crate) tolerance: f64, - pub(crate) norm_fpr: f64, +pub struct FBSCache +where + T: OptFloat, +{ + pub(crate) work_gradient_u: Vec, + pub(crate) work_u_previous: Vec, + pub(crate) gamma: T, + pub(crate) tolerance: T, + pub(crate) norm_fpr: T, } -impl FBSCache { +impl FBSCache +where + T: OptFloat, +{ /// Construct a new instance of `FBSCache` /// /// ## Arguments @@ -37,13 +45,13 @@ impl FBSCache { /// This method will panic if there is no available memory for the required allocation /// (capacity overflow) /// - pub fn new(n: NonZeroUsize, gamma: f64, tolerance: f64) -> FBSCache { + pub fn new(n: NonZeroUsize, gamma: T, tolerance: T) -> FBSCache { FBSCache { - work_gradient_u: vec![0.0; n.get()], - work_u_previous: vec![0.0; n.get()], + work_gradient_u: vec![T::zero(); n.get()], + work_u_previous: vec![T::zero(); n.get()], gamma, tolerance, - norm_fpr: std::f64::INFINITY, + norm_fpr: T::infinity(), } } } diff --git a/src/core/fbs/fbs_engine.rs b/src/core/fbs/fbs_engine.rs index 717e70c6..6d47098f 100644 --- a/src/core/fbs/fbs_engine.rs +++ b/src/core/fbs/fbs_engine.rs @@ -1,29 +1,29 @@ //! FBS Engine //! -use crate::{ - constraints, - core::{fbs::FBSCache, AlgorithmEngine, Problem}, - matrix_operations, FunctionCallResult, SolverError, -}; +use crate::core::fbs::FBSCache; +use crate::core::{AlgorithmEngine, OptFloat, Problem}; +use crate::{constraints, matrix_operations, FunctionCallResult, SolverError}; /// The FBE Engine defines the steps of the FBE algorithm and the termination criterion /// -pub struct FBSEngine<'a, GradientType, ConstraintType, CostType> +pub struct FBSEngine<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { - pub(crate) problem: Problem<'a, GradientType, ConstraintType, CostType>, - pub(crate) cache: &'a mut FBSCache, + pub(crate) problem: Problem<'a, GradientType, ConstraintType, CostType, T>, + pub(crate) cache: &'a mut FBSCache, } -impl<'a, GradientType, ConstraintType, CostType> - FBSEngine<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> + FBSEngine<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { /// Constructor for instances of `FBSEngine` /// @@ -36,13 +36,13 @@ where /// /// An new instance of `FBSEngine` pub fn new( - problem: Problem<'a, GradientType, ConstraintType, CostType>, - cache: &'a mut FBSCache, - ) -> FBSEngine<'a, GradientType, ConstraintType, CostType> { + problem: Problem<'a, GradientType, ConstraintType, CostType, T>, + cache: &'a mut FBSCache, + ) -> FBSEngine<'a, GradientType, ConstraintType, CostType, T> { FBSEngine { problem, cache } } - fn gradient_step(&mut self, u_current: &mut [f64]) { + fn gradient_step(&mut self, u_current: &mut [T]) { assert_eq!( Ok(()), (self.problem.gradf)(u_current, &mut self.cache.work_gradient_u), @@ -56,17 +56,18 @@ where .for_each(|(u, w)| *u -= self.cache.gamma * *w); } - fn projection_step(&mut self, u_current: &mut [f64]) { + fn projection_step(&mut self, u_current: &mut [T]) { self.problem.constraints.project(u_current); } } -impl<'a, GradientType, ConstraintType, CostType> AlgorithmEngine - for FBSEngine<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> AlgorithmEngine + for FBSEngine<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult + 'a, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult + 'a, - ConstraintType: constraints::Constraint + 'a, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult + 'a, + CostType: Fn(&[T], &mut T) -> FunctionCallResult + 'a, + ConstraintType: constraints::Constraint + 'a, + T: OptFloat, { /// Take a forward-backward step and check whether the algorithm should terminate /// @@ -83,7 +84,7 @@ where /// /// The method may panick if the computation of the gradient of the cost function /// or the cost function panics. - fn step(&mut self, u_current: &mut [f64]) -> Result { + fn step(&mut self, u_current: &mut [T]) -> Result { self.cache.work_u_previous.copy_from_slice(u_current); // cache the previous step self.gradient_step(u_current); // compute the gradient self.projection_step(u_current); // project @@ -93,7 +94,7 @@ where Ok(self.cache.norm_fpr > self.cache.tolerance) } - fn init(&mut self, _u_current: &mut [f64]) -> FunctionCallResult { + fn init(&mut self, _u_current: &mut [T]) -> FunctionCallResult { Ok(()) } } diff --git a/src/core/fbs/fbs_optimizer.rs b/src/core/fbs/fbs_optimizer.rs index d714ab67..cb892a9a 100644 --- a/src/core/fbs/fbs_optimizer.rs +++ b/src/core/fbs/fbs_optimizer.rs @@ -1,15 +1,13 @@ //! FBS Algorithm //! -use crate::{ - constraints, - core::{ - fbs::fbs_engine::FBSEngine, fbs::FBSCache, AlgorithmEngine, ExitStatus, Optimizer, Problem, - SolverStatus, - }, - matrix_operations, FunctionCallResult, SolverError, -}; + use std::time; +use crate::core::fbs::fbs_engine::FBSEngine; +use crate::core::fbs::FBSCache; +use crate::core::{AlgorithmEngine, ExitStatus, OptFloat, Optimizer, Problem, SolverStatus}; +use crate::{constraints, matrix_operations, FunctionCallResult, SolverError}; + const MAX_ITER: usize = 100_usize; /// Optimiser using forward-backward splitting iterations (projected gradient) @@ -22,23 +20,25 @@ const MAX_ITER: usize = 100_usize; /// a different optimization problem. /// /// -pub struct FBSOptimizer<'a, GradientType, ConstraintType, CostType> +pub struct FBSOptimizer<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { - fbs_engine: FBSEngine<'a, GradientType, ConstraintType, CostType>, + fbs_engine: FBSEngine<'a, GradientType, ConstraintType, CostType, T>, max_iter: usize, max_duration: Option, } -impl<'a, GradientType, ConstraintType, CostType> - FBSOptimizer<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> + FBSOptimizer<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { /// Constructs a new instance of `FBSOptimizer` /// @@ -47,8 +47,8 @@ where /// - `problem`: problem definition /// - `cache`: instance of `FBSCache` pub fn new( - problem: Problem<'a, GradientType, ConstraintType, CostType>, - cache: &'a mut FBSCache, + problem: Problem<'a, GradientType, ConstraintType, CostType, T>, + cache: &'a mut FBSCache, ) -> Self { FBSOptimizer { fbs_engine: FBSEngine::new(problem, cache), @@ -64,9 +64,9 @@ where /// The method panics if the specified tolerance is not positive pub fn with_tolerance( self, - tolerance: f64, - ) -> FBSOptimizer<'a, GradientType, ConstraintType, CostType> { - assert!(tolerance > 0.0); + tolerance: T, + ) -> FBSOptimizer<'a, GradientType, ConstraintType, CostType, T> { + assert!(tolerance > T::zero()); self.fbs_engine.cache.tolerance = tolerance; self @@ -76,7 +76,7 @@ where pub fn with_max_iter( mut self, max_iter: usize, - ) -> FBSOptimizer<'a, GradientType, ConstraintType, CostType> { + ) -> FBSOptimizer<'a, GradientType, ConstraintType, CostType, T> { self.max_iter = max_iter; self } @@ -85,20 +85,21 @@ where pub fn with_max_duration( mut self, max_duration: time::Duration, - ) -> FBSOptimizer<'a, GradientType, ConstraintType, CostType> { + ) -> FBSOptimizer<'a, GradientType, ConstraintType, CostType, T> { self.max_duration = Some(max_duration); self } } -impl<'life, GradientType, ConstraintType, CostType> Optimizer - for FBSOptimizer<'life, GradientType, ConstraintType, CostType> +impl<'life, GradientType, ConstraintType, CostType, T> Optimizer + for FBSOptimizer<'life, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult + 'life, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult + 'life, - ConstraintType: constraints::Constraint + 'life, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult + 'life, + CostType: Fn(&[T], &mut T) -> FunctionCallResult + 'life, + ConstraintType: constraints::Constraint + 'life, + T: OptFloat, { - fn solve(&mut self, u: &mut [f64]) -> Result { + fn solve(&mut self, u: &mut [T]) -> Result, SolverError> { let now = instant::Instant::now(); // Initialize - propagate error upstream, if any @@ -120,7 +121,7 @@ where } // cost at the solution [propagate error upstream] - let mut cost_value: f64 = 0.0; + let mut cost_value: T = T::zero(); (self.fbs_engine.problem.cost)(u, &mut cost_value)?; if !matrix_operations::is_finite(u) || !cost_value.is_finite() { diff --git a/src/core/fbs/tests.rs b/src/core/fbs/tests.rs index b94425b1..554b02cc 100644 --- a/src/core/fbs/tests.rs +++ b/src/core/fbs/tests.rs @@ -1,8 +1,9 @@ +use std::num::NonZeroUsize; + use super::super::*; use super::*; use crate::constraints; use crate::core::fbs::fbs_engine::FBSEngine; -use std::num::NonZeroUsize; const N_DIM: usize = 2; diff --git a/src/core/mod.rs b/src/core/mod.rs index 46099e1a..a2e47bd2 100644 --- a/src/core/mod.rs +++ b/src/core/mod.rs @@ -4,11 +4,13 @@ //! pub mod fbs; +pub mod opt_float; pub mod panoc; pub mod problem; pub mod solver_status; pub use crate::{constraints, FunctionCallResult, SolverError}; +pub use opt_float::OptFloat; pub use problem::Problem; pub use solver_status::SolverStatus; @@ -29,12 +31,15 @@ pub enum ExitStatus { } /// A general optimizer -pub trait Optimizer { +pub trait Optimizer +where + T: OptFloat + std::fmt::Debug, +{ /// solves a given problem and updates the initial estimate `u` with the solution /// /// Returns the solver status /// - fn solve(&mut self, u: &mut [f64]) -> Result; + fn solve(&mut self, u: &mut [T]) -> Result, SolverError>; } /// Engine supporting an algorithm @@ -46,10 +51,13 @@ pub trait Optimizer { /// It defines what the algorithm does at every step (see `step`) and whether /// the specified termination criterion is satisfied /// -pub trait AlgorithmEngine { +pub trait AlgorithmEngine +where + T: OptFloat + std::fmt::Debug, +{ /// Take a step of the algorithm and return `Ok(true)` only if the iterations should continue - fn step(&mut self, u: &mut [f64]) -> Result; + fn step(&mut self, u: &mut [T]) -> Result; /// Initializes the algorithm - fn init(&mut self, u: &mut [f64]) -> FunctionCallResult; + fn init(&mut self, u: &mut [T]) -> FunctionCallResult; } diff --git a/src/core/opt_float.rs b/src/core/opt_float.rs new file mode 100644 index 00000000..ac9a92ce --- /dev/null +++ b/src/core/opt_float.rs @@ -0,0 +1,113 @@ +//! OptFloat trait that combines Float functionality with optimization constants (and a few traits) +use num::Float; + +/// Trait that combines Float functionality with optimization-specific constants +/// This allows different float types to have different optimization parameters +pub trait OptFloat: + Float + + std::iter::Sum + + num::FromPrimitive + + num::ToPrimitive + + std::fmt::Debug + + std::ops::AddAssign + + std::ops::SubAssign + + std::ops::MulAssign + + std::ops::DivAssign +{ + /// Minimum estimated Lipschitz constant (initial estimate) + fn min_l_estimate() -> Self; + + /// gamma = GAMMA_L_COEFF/L + fn gamma_l_coeff() -> Self; + + /// Delta in the estimation of the initial Lipschitz constant + fn delta_lipschitz() -> Self; + + /// Epsilon in the estimation of the initial Lipschitz constant + fn epsilon_lipschitz() -> Self; + + /// Safety parameter used to check a strict inequality in the update of the Lipschitz constant + fn lipschitz_update_epsilon() -> Self; + + /// Maximum possible Lipschitz constant + fn max_lipschitz_constant() -> Self; +} + +/// Default implementation for f64 with original constants +impl OptFloat for f64 { + fn min_l_estimate() -> Self { + 1e-10 + } + + fn gamma_l_coeff() -> Self { + 0.95 + } + + fn delta_lipschitz() -> Self { + 1e-12 + } + + fn epsilon_lipschitz() -> Self { + 1e-6 + } + + fn lipschitz_update_epsilon() -> Self { + 1e-6 + } + + fn max_lipschitz_constant() -> Self { + 1e9 + } +} + +/// Default implementation for f32 with scaled constants +impl OptFloat for f32 { + fn min_l_estimate() -> Self { + 8.74e-6 + } + + fn gamma_l_coeff() -> Self { + 0.95 + } + + fn delta_lipschitz() -> Self { + 1.32e-6 + } + + fn epsilon_lipschitz() -> Self { + 7.32e-4 + } + + fn lipschitz_update_epsilon() -> Self { + 2.62e-4 + } + + fn max_lipschitz_constant() -> Self { + 1e9 + } +} + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_f64_constants() { + assert_eq!(f64::min_l_estimate(), 1e-10); + assert_eq!(f64::gamma_l_coeff(), 0.95); + assert_eq!(f64::delta_lipschitz(), 1e-12); + assert_eq!(f64::epsilon_lipschitz(), 1e-6); + assert_eq!(f64::lipschitz_update_epsilon(), 1e-6); + assert_eq!(f64::max_lipschitz_constant(), 1e9); + } + + #[test] + fn test_f32_constants() { + assert_eq!(f32::min_l_estimate(), 8.74e-6); + assert_eq!(f32::gamma_l_coeff(), 0.95); + assert_eq!(f32::delta_lipschitz(), 1.32e-6); + assert_eq!(f32::epsilon_lipschitz(), 7.32e-4); + assert_eq!(f32::lipschitz_update_epsilon(), 2.62e-4); + assert_eq!(f32::max_lipschitz_constant(), 1e9); + } +} diff --git a/src/core/panoc/mod.rs b/src/core/panoc/mod.rs index 36773a64..3c747bbb 100644 --- a/src/core/panoc/mod.rs +++ b/src/core/panoc/mod.rs @@ -11,3 +11,6 @@ pub use panoc_optimizer::PANOCOptimizer; #[cfg(test)] mod tests; + +#[cfg(test)] +mod tests_f32; diff --git a/src/core/panoc/panoc_cache.rs b/src/core/panoc/panoc_cache.rs index 3f3c0b13..ddabbddb 100644 --- a/src/core/panoc/panoc_cache.rs +++ b/src/core/panoc/panoc_cache.rs @@ -1,3 +1,5 @@ +use crate::core::OptFloat; + const DEFAULT_SY_EPSILON: f64 = 1e-10; const DEFAULT_CBFGS_EPSILON: f64 = 1e-8; const DEFAULT_CBFGS_ALPHA: f64 = 1.0; @@ -12,32 +14,38 @@ const DEFAULT_CBFGS_ALPHA: f64 = 1.0; /// Subsequently, a `PANOCEngine` is used to construct an instance of `PANOCAlgorithm` /// #[derive(Debug)] -pub struct PANOCCache { - pub(crate) lbfgs: lbfgs::Lbfgs, - pub(crate) gradient_u: Vec, +pub struct PANOCCache +where + T: OptFloat, +{ + pub(crate) lbfgs: lbfgs::Lbfgs, + pub(crate) gradient_u: Vec, /// Stores the gradient of the cost at the previous iteration. This is /// an optional field because it is used (and needs to be allocated) /// only if we need to check the AKKT-specific termination conditions - pub(crate) gradient_u_previous: Option>, - pub(crate) u_half_step: Vec, - pub(crate) gradient_step: Vec, - pub(crate) direction_lbfgs: Vec, - pub(crate) u_plus: Vec, - pub(crate) rhs_ls: f64, - pub(crate) lhs_ls: f64, - pub(crate) gamma_fpr: Vec, - pub(crate) gamma: f64, - pub(crate) tolerance: f64, - pub(crate) norm_gamma_fpr: f64, - pub(crate) tau: f64, - pub(crate) lipschitz_constant: f64, - pub(crate) sigma: f64, - pub(crate) cost_value: f64, + pub(crate) gradient_u_previous: Option>, + pub(crate) u_half_step: Vec, + pub(crate) gradient_step: Vec, + pub(crate) direction_lbfgs: Vec, + pub(crate) u_plus: Vec, + pub(crate) rhs_ls: T, + pub(crate) lhs_ls: T, + pub(crate) gamma_fpr: Vec, + pub(crate) gamma: T, + pub(crate) tolerance: T, + pub(crate) norm_gamma_fpr: T, + pub(crate) tau: T, + pub(crate) lipschitz_constant: T, + pub(crate) sigma: T, + pub(crate) cost_value: T, pub(crate) iteration: usize, - pub(crate) akkt_tolerance: Option, + pub(crate) akkt_tolerance: Option, } -impl PANOCCache { +impl PANOCCache +where + T: OptFloat, +{ /// Construct a new instance of `PANOCCache` /// /// ## Arguments @@ -59,30 +67,30 @@ impl PANOCCache { /// /// It allocates a total of `8*problem_size + 2*lbfgs_memory_size*problem_size + 2*lbfgs_memory_size + 11` floats (`f64`) /// - pub fn new(problem_size: usize, tolerance: f64, lbfgs_memory_size: usize) -> PANOCCache { - assert!(tolerance > 0., "tolerance must be positive"); + pub fn new(problem_size: usize, tolerance: T, lbfgs_memory_size: usize) -> PANOCCache { + assert!(tolerance > T::zero(), "tolerance must be positive"); PANOCCache { - gradient_u: vec![0.0; problem_size], + gradient_u: vec![T::zero(); problem_size], gradient_u_previous: None, - u_half_step: vec![0.0; problem_size], - gamma_fpr: vec![0.0; problem_size], - direction_lbfgs: vec![0.0; problem_size], - gradient_step: vec![0.0; problem_size], - u_plus: vec![0.0; problem_size], - gamma: 0.0, + u_half_step: vec![T::zero(); problem_size], + gamma_fpr: vec![T::zero(); problem_size], + direction_lbfgs: vec![T::zero(); problem_size], + gradient_step: vec![T::zero(); problem_size], + u_plus: vec![T::zero(); problem_size], + gamma: T::zero(), tolerance, - norm_gamma_fpr: std::f64::INFINITY, - lbfgs: lbfgs::Lbfgs::new(problem_size, lbfgs_memory_size) - .with_cbfgs_alpha(DEFAULT_CBFGS_ALPHA) - .with_cbfgs_epsilon(DEFAULT_CBFGS_EPSILON) - .with_sy_epsilon(DEFAULT_SY_EPSILON), - lhs_ls: 0.0, - rhs_ls: 0.0, - tau: 1.0, - lipschitz_constant: 0.0, - sigma: 0.0, - cost_value: 0.0, + norm_gamma_fpr: T::infinity(), + lbfgs: lbfgs::Lbfgs::::new(problem_size, lbfgs_memory_size) + .with_cbfgs_alpha(T::from(DEFAULT_CBFGS_ALPHA).unwrap()) + .with_cbfgs_epsilon(T::from(DEFAULT_CBFGS_EPSILON).unwrap()) + .with_sy_epsilon(T::from(DEFAULT_SY_EPSILON).unwrap()), + lhs_ls: T::zero(), + rhs_ls: T::zero(), + tau: T::one(), + lipschitz_constant: T::zero(), + sigma: T::zero(), + cost_value: T::zero(), iteration: 0, akkt_tolerance: None, } @@ -99,10 +107,13 @@ impl PANOCCache { /// /// The method panics if `akkt_tolerance` is nonpositive /// - pub fn set_akkt_tolerance(&mut self, akkt_tolerance: f64) { - assert!(akkt_tolerance > 0.0, "akkt_tolerance must be positive"); + pub fn set_akkt_tolerance(&mut self, akkt_tolerance: T) { + assert!( + akkt_tolerance > T::zero(), + "akkt_tolerance must be positive" + ); self.akkt_tolerance = Some(akkt_tolerance); - self.gradient_u_previous = Some(vec![0.0; self.gradient_step.len()]); + self.gradient_u_previous = Some(vec![T::zero(); self.gradient_step.len()]); } /// Copies the value of the current cost gradient to `gradient_u_previous`, @@ -117,8 +128,8 @@ impl PANOCCache { } /// Computes the AKKT residual which is defined as `||gamma*(fpr + df - df_previous)||` - fn akkt_residual(&self) -> f64 { - let mut r = 0.0; + fn akkt_residual(&self) -> T { + let mut r = T::zero(); if let Some(df_previous) = &self.gradient_u_previous { // Notation: gamma_fpr_i is the i-th element of gamma_fpr = gamma * fpr, // df_i is the i-th element of the gradient of the cost function at the @@ -129,7 +140,7 @@ impl PANOCCache { .iter() .zip(self.gradient_u.iter()) .zip(df_previous.iter()) - .fold(0.0, |mut sum, ((&gamma_fpr_i, &df_i), &dfp_i)| { + .fold(T::zero(), |mut sum, ((&gamma_fpr_i, &df_i), &dfp_i)| { sum += (gamma_fpr_i + self.gamma * (df_i - dfp_i)).powi(2); sum }) @@ -175,14 +186,14 @@ impl PANOCCache { /// and `gamma` to 0.0 pub fn reset(&mut self) { self.lbfgs.reset(); - self.lhs_ls = 0.0; - self.rhs_ls = 0.0; - self.tau = 1.0; - self.lipschitz_constant = 0.0; - self.sigma = 0.0; - self.cost_value = 0.0; + self.lhs_ls = T::zero(); + self.rhs_ls = T::zero(); + self.tau = T::one(); + self.lipschitz_constant = T::zero(); + self.sigma = T::zero(); + self.cost_value = T::zero(); self.iteration = 0; - self.gamma = 0.0; + self.gamma = T::zero(); } /// Sets the CBFGS parameters `alpha` and `epsilon` @@ -202,7 +213,7 @@ impl PANOCCache { /// The method panics if alpha or epsilon are nonpositive and if sy_epsilon /// is negative. /// - pub fn with_cbfgs_parameters(mut self, alpha: f64, epsilon: f64, sy_epsilon: f64) -> Self { + pub fn with_cbfgs_parameters(mut self, alpha: T, epsilon: T, sy_epsilon: T) -> Self { self.lbfgs = self .lbfgs .with_cbfgs_alpha(alpha) diff --git a/src/core/panoc/panoc_engine.rs b/src/core/panoc/panoc_engine.rs index bb4ad484..d1c322f7 100644 --- a/src/core/panoc/panoc_engine.rs +++ b/src/core/panoc/panoc_engine.rs @@ -1,52 +1,32 @@ -use crate::{ - constraints, - core::{panoc::PANOCCache, AlgorithmEngine, Problem}, - matrix_operations, FunctionCallResult, SolverError, -}; - -/// Mimum estimated Lipschitz constant (initial estimate) -const MIN_L_ESTIMATE: f64 = 1e-10; - -/// gamma = GAMMA_L_COEFF/L -const GAMMA_L_COEFF: f64 = 0.95; - -//const SIGMA_COEFF: f64 = 0.49; - -/// Delta in the estimation of the initial Lipschitz constant -const DELTA_LIPSCHITZ: f64 = 1e-12; - -/// Epsilon in the estimation of the initial Lipschitz constant -const EPSILON_LIPSCHITZ: f64 = 1e-6; - -/// Safety parameter used to check a strict inequality in the update of the Lipschitz constant -const LIPSCHITZ_UPDATE_EPSILON: f64 = 1e-6; +use crate::core::panoc::PANOCCache; +use crate::core::{AlgorithmEngine, OptFloat, Problem}; +use crate::{constraints, matrix_operations, FunctionCallResult, SolverError}; /// Maximum iterations of updating the Lipschitz constant const MAX_LIPSCHITZ_UPDATE_ITERATIONS: usize = 10; -/// Maximum possible Lipschitz constant -const MAX_LIPSCHITZ_CONSTANT: f64 = 1e9; - /// Maximum number of linesearch iterations const MAX_LINESEARCH_ITERATIONS: u32 = 10; /// Engine for PANOC algorithm -pub struct PANOCEngine<'a, GradientType, ConstraintType, CostType> +pub struct PANOCEngine<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { - problem: Problem<'a, GradientType, ConstraintType, CostType>, - pub(crate) cache: &'a mut PANOCCache, + problem: Problem<'a, GradientType, ConstraintType, CostType, T>, + pub(crate) cache: &'a mut PANOCCache, } -impl<'a, GradientType, ConstraintType, CostType> - PANOCEngine<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> + PANOCEngine<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { /// Construct a new Engine for PANOC /// @@ -60,28 +40,28 @@ where /// /// pub fn new( - problem: Problem<'a, GradientType, ConstraintType, CostType>, - cache: &'a mut PANOCCache, - ) -> PANOCEngine<'a, GradientType, ConstraintType, CostType> { + problem: Problem<'a, GradientType, ConstraintType, CostType, T>, + cache: &'a mut PANOCCache, + ) -> PANOCEngine<'a, GradientType, ConstraintType, CostType, T> { PANOCEngine { problem, cache } } /// Estimate the local Lipschitz constant at `u` - fn estimate_loc_lip(&mut self, u: &mut [f64]) -> FunctionCallResult { + fn estimate_loc_lip(&mut self, u: &mut [T]) -> FunctionCallResult { let mut lipest = crate::lipschitz_estimator::LipschitzEstimator::new( u, &self.problem.gradf, &mut self.cache.gradient_u, ) - .with_delta(DELTA_LIPSCHITZ) - .with_epsilon(EPSILON_LIPSCHITZ); + .with_delta(T::delta_lipschitz()) + .with_epsilon(T::epsilon_lipschitz()); self.cache.lipschitz_constant = lipest.estimate_local_lipschitz()?; Ok(()) } /// Computes the FPR and its norm - fn compute_fpr(&mut self, u_current: &[f64]) { + fn compute_fpr(&mut self, u_current: &[T]) { // compute the FPR: // fpr ← u - u_half_step let cache = &mut self.cache; @@ -90,13 +70,13 @@ where .iter_mut() .zip(u_current.iter()) .zip(cache.u_half_step.iter()) - .for_each(|((fpr, u), uhalf)| *fpr = u - uhalf); + .for_each(|((fpr, u), uhalf)| *fpr = *u - *uhalf); // compute the norm of FPR cache.norm_gamma_fpr = matrix_operations::norm2(&cache.gamma_fpr); } /// Computes a gradient step; does not compute the gradient - fn gradient_step(&mut self, u_current: &[f64]) { + fn gradient_step(&mut self, u_current: &[T]) { // take a gradient step: // gradient_step ← u_current - gamma * gradient let cache = &mut self.cache; @@ -132,7 +112,7 @@ where } /// Computes an LBFGS direction; updates `cache.direction_lbfgs` - fn lbfgs_direction(&mut self, u_current: &[f64]) { + fn lbfgs_direction(&mut self, u_current: &[T]) { let cache = &mut self.cache; // update the LBFGS buffer cache.lbfgs.update_hessian(&cache.gamma_fpr, u_current); @@ -147,7 +127,7 @@ where /// Returns the RHS of the Lipschitz update /// Computes rhs = cost + LIP_EPS * |f| - gamma * + (L/2/gamma) ||gamma * fpr||^2 - fn lipschitz_check_rhs(&mut self) -> f64 { + fn lipschitz_check_rhs(&mut self) -> T { let cache = &mut self.cache; let gamma = cache.gamma; let cost_value = cache.cost_value; @@ -156,13 +136,14 @@ where matrix_operations::inner_product(&cache.gradient_u, &cache.gamma_fpr); // rhs ← cost + LIP_EPS * |f| - + (L/2/gamma) ||gamma_fpr||^2 - cost_value + LIPSCHITZ_UPDATE_EPSILON * cost_value.abs() - inner_prod_grad_fpr - + (GAMMA_L_COEFF / (2.0 * gamma)) * (cache.norm_gamma_fpr.powi(2)) + cost_value + T::lipschitz_update_epsilon() * cost_value.abs() - inner_prod_grad_fpr + + (T::gamma_l_coeff() / (T::from(2.0).unwrap() * gamma)) + * (cache.norm_gamma_fpr.powi(2)) } /// Updates the estimate of the Lipscthiz constant - fn update_lipschitz_constant(&mut self, u_current: &[f64]) -> FunctionCallResult { - let mut cost_u_half_step = 0.0; + fn update_lipschitz_constant(&mut self, u_current: &[T]) -> FunctionCallResult { + let mut cost_u_half_step = T::zero(); // Compute the cost at the half step (self.problem.cost)(&self.cache.u_half_step, &mut cost_u_half_step)?; @@ -174,13 +155,13 @@ where while cost_u_half_step > self.lipschitz_check_rhs() && it_lipschitz_search < MAX_LIPSCHITZ_UPDATE_ITERATIONS - && self.cache.lipschitz_constant < MAX_LIPSCHITZ_CONSTANT + && self.cache.lipschitz_constant < T::max_lipschitz_constant() { self.cache.lbfgs.reset(); // invalidate the L-BFGS buffer // update L, sigma and gamma... - self.cache.lipschitz_constant *= 2.; - self.cache.gamma /= 2.; + self.cache.lipschitz_constant = self.cache.lipschitz_constant * T::from(2.0).unwrap(); + self.cache.gamma = self.cache.gamma / T::from(2.0).unwrap(); // recompute the half step... self.gradient_step(u_current); // updates self.cache.gradient_step @@ -194,17 +175,18 @@ where self.compute_fpr(u_current); it_lipschitz_search += 1; } - self.cache.sigma = (1.0 - GAMMA_L_COEFF) / (4.0 * self.cache.gamma); + self.cache.sigma = + (T::one() - T::gamma_l_coeff()) / (T::from(4.0).unwrap() * self.cache.gamma); Ok(()) } /// Computes u_plus ← u - gamma * (1-tau) * fpr - tau * dir, - fn compute_u_plus(&mut self, u: &[f64]) { + fn compute_u_plus(&mut self, u: &[T]) { let cache = &mut self.cache; let _gamma = cache.gamma; let tau = cache.tau; - let temp_ = 1.0 - tau; + let temp_ = T::from(1.0).unwrap() - tau; cache .u_plus .iter_mut() @@ -228,15 +210,17 @@ where // + 0.5 * dist squared / gamma // - sigma * norm_gamma_fpr^2 let fbe = cache.cost_value - - 0.5 * cache.gamma * matrix_operations::norm2_squared(&cache.gradient_u) - + 0.5 * dist_squared / cache.gamma; + - T::from(0.5).unwrap() + * cache.gamma + * matrix_operations::norm2_squared(&cache.gradient_u) + + T::from(0.5).unwrap() * dist_squared / cache.gamma; let sigma_fpr_sq = cache.sigma * cache.norm_gamma_fpr.powi(2); cache.rhs_ls = fbe - sigma_fpr_sq; } /// Computes the left hand side of the line search condition and compares it with the RHS; /// returns `true` if and only if lhs > rhs (when the line search should continue) - fn line_search_condition(&mut self, u: &[f64]) -> Result { + fn line_search_condition(&mut self, u: &[T]) -> Result { let gamma = self.cache.gamma; // u_plus ← u - (1-tau)*gamma_fpr + tau*direction @@ -259,14 +243,16 @@ where // Update the LHS of the line search condition self.cache.lhs_ls = self.cache.cost_value - - 0.5 * gamma * matrix_operations::norm2_squared(&self.cache.gradient_u) - + 0.5 * dist_squared / self.cache.gamma; + - T::from(0.5).unwrap() + * gamma + * matrix_operations::norm2_squared(&self.cache.gradient_u) + + T::from(0.5).unwrap() * dist_squared / self.cache.gamma; Ok(self.cache.lhs_ls > self.cache.rhs_ls) } /// Update without performing a line search; this is executed at the first iteration - fn update_no_linesearch(&mut self, u_current: &mut [f64]) -> FunctionCallResult { + fn update_no_linesearch(&mut self, u_current: &mut [T]) -> FunctionCallResult { u_current.copy_from_slice(&self.cache.u_half_step); // set u_current ← u_half_step (self.problem.cost)(u_current, &mut self.cache.cost_value)?; // cost value (self.problem.gradf)(u_current, &mut self.cache.gradient_u)?; // compute gradient @@ -277,17 +263,17 @@ where } /// Performs a line search to select tau - fn linesearch(&mut self, u_current: &mut [f64]) -> FunctionCallResult { + fn linesearch(&mut self, u_current: &mut [T]) -> FunctionCallResult { // perform line search self.compute_rhs_ls(); // compute the right hand side of the line search - self.cache.tau = 1.0; // initialise tau ← 1.0 + self.cache.tau = T::from(1.0).unwrap(); // initialise tau ← 1.0 let mut num_ls_iters = 0; while self.line_search_condition(u_current)? && num_ls_iters < MAX_LINESEARCH_ITERATIONS { - self.cache.tau /= 2.0; + self.cache.tau = self.cache.tau / T::from(2.0).unwrap(); num_ls_iters += 1; } if num_ls_iters == MAX_LINESEARCH_ITERATIONS { - self.cache.tau = 0.; + self.cache.tau = T::zero(); u_current.copy_from_slice(&self.cache.u_half_step); } // Sets `u_current` to `u_plus` (u_current ← u_plus) @@ -298,12 +284,13 @@ where } /// Implementation of the `step` and `init` methods of [trait.AlgorithmEngine.html] -impl<'a, GradientType, ConstraintType, CostType> AlgorithmEngine - for PANOCEngine<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> AlgorithmEngine + for PANOCEngine<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat + std::fmt::Debug, { /// PANOC step /// @@ -315,7 +302,7 @@ where /// iterate of PANOC /// /// - fn step(&mut self, u_current: &mut [f64]) -> Result { + fn step(&mut self, u_current: &mut [T]) -> Result { // caches the previous gradient vector (copies df to df_previous) self.cache.cache_previous_gradient(); @@ -349,12 +336,14 @@ where /// gradient of the cost at the initial point, initial estimates for `gamma` and `sigma`, /// a gradient step and a half step (projected gradient step) /// - fn init(&mut self, u_current: &mut [f64]) -> FunctionCallResult { + fn init(&mut self, u_current: &mut [T]) -> FunctionCallResult { self.cache.reset(); (self.problem.cost)(u_current, &mut self.cache.cost_value)?; // cost value self.estimate_loc_lip(u_current)?; // computes the gradient as well! (self.cache.gradient_u) - self.cache.gamma = GAMMA_L_COEFF / f64::max(self.cache.lipschitz_constant, MIN_L_ESTIMATE); - self.cache.sigma = (1.0 - GAMMA_L_COEFF) / (4.0 * self.cache.gamma); + self.cache.gamma = + T::gamma_l_coeff() / T::max(self.cache.lipschitz_constant, T::min_l_estimate()); + self.cache.sigma = + (T::one() - T::gamma_l_coeff()) / (T::from(4.0).unwrap() * self.cache.gamma); self.gradient_step(u_current); // updated self.cache.gradient_step self.half_step(); // updates self.cache.u_half_step @@ -368,11 +357,10 @@ where #[cfg(test)] mod tests { - use crate::constraints; use crate::core::panoc::panoc_engine::PANOCEngine; use crate::core::panoc::*; use crate::core::Problem; - use crate::mocks; + use crate::{constraints, mocks}; #[test] fn t_compute_fpr() { diff --git a/src/core/panoc/panoc_optimizer.rs b/src/core/panoc/panoc_optimizer.rs index 602abd45..921e4fbf 100644 --- a/src/core/panoc/panoc_optimizer.rs +++ b/src/core/panoc/panoc_optimizer.rs @@ -1,37 +1,36 @@ //! PANOC optimizer //! -use crate::{ - constraints, - core::{ - panoc::panoc_engine::PANOCEngine, panoc::PANOCCache, AlgorithmEngine, ExitStatus, - Optimizer, Problem, SolverStatus, - }, - matrix_operations, FunctionCallResult, SolverError, -}; use std::time; +use crate::core::panoc::panoc_engine::PANOCEngine; +use crate::core::panoc::PANOCCache; +use crate::core::{AlgorithmEngine, ExitStatus, OptFloat, Optimizer, Problem, SolverStatus}; +use crate::{constraints, matrix_operations, FunctionCallResult, SolverError}; + const MAX_ITER: usize = 100_usize; /// Optimizer using the PANOC algorithm /// /// -pub struct PANOCOptimizer<'a, GradientType, ConstraintType, CostType> +pub struct PANOCOptimizer<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { - panoc_engine: PANOCEngine<'a, GradientType, ConstraintType, CostType>, + panoc_engine: PANOCEngine<'a, GradientType, ConstraintType, CostType, T>, max_iter: usize, max_duration: Option, } -impl<'a, GradientType, ConstraintType, CostType> - PANOCOptimizer<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> + PANOCOptimizer<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { /// Constructor of `PANOCOptimizer` /// @@ -44,8 +43,8 @@ where /// /// Does not panic pub fn new( - problem: Problem<'a, GradientType, ConstraintType, CostType>, - cache: &'a mut PANOCCache, + problem: Problem<'a, GradientType, ConstraintType, CostType, T>, + cache: &'a mut PANOCCache, ) -> Self { PANOCOptimizer { panoc_engine: PANOCEngine::new(problem, cache), @@ -62,8 +61,8 @@ where /// ## Panics /// /// The method panics if the specified tolerance is not positive - pub fn with_tolerance(self, tolerance: f64) -> Self { - assert!(tolerance > 0.0, "tolerance must be larger than 0"); + pub fn with_tolerance(self, tolerance: T) -> Self { + assert!(tolerance > T::zero(), "tolerance must be larger than 0"); self.panoc_engine.cache.tolerance = tolerance; self @@ -90,8 +89,11 @@ where /// The method panics if the provided value of the AKKT-specific tolerance is /// not positive. /// - pub fn with_akkt_tolerance(self, akkt_tolerance: f64) -> Self { - assert!(akkt_tolerance > 0.0, "akkt_tolerance must be positive"); + pub fn with_akkt_tolerance(self, akkt_tolerance: T) -> Self { + assert!( + akkt_tolerance > T::zero(), + "akkt_tolerance must be positive" + ); self.panoc_engine.cache.set_akkt_tolerance(akkt_tolerance); self } @@ -115,14 +117,15 @@ where } } -impl<'life, GradientType, ConstraintType, CostType> Optimizer - for PANOCOptimizer<'life, GradientType, ConstraintType, CostType> +impl<'life, GradientType, ConstraintType, CostType, T> Optimizer + for PANOCOptimizer<'life, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult + 'life, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint + 'life, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult + 'life, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint + 'life, + T: OptFloat, { - fn solve(&mut self, u: &mut [f64]) -> Result { + fn solve(&mut self, u: &mut [T]) -> Result, SolverError> { let now = instant::Instant::now(); /* diff --git a/src/core/panoc/tests.rs b/src/core/panoc/tests.rs index afd89ef3..f7afbe2d 100644 --- a/src/core/panoc/tests.rs +++ b/src/core/panoc/tests.rs @@ -36,18 +36,22 @@ fn t_panoc_init() { println!("cache = {:#?}", &panoc_cache); } -fn print_panoc_engine( - panoc_engine: &PANOCEngine, +fn print_panoc_engine( + panoc_engine: &PANOCEngine, ) where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat + std::fmt::Debug, { println!("> fpr = {:?}", &panoc_engine.cache.gamma_fpr); - println!("> fpr = {:.2e}", panoc_engine.cache.norm_gamma_fpr); - println!("> L = {:.3}", panoc_engine.cache.lipschitz_constant); - println!("> gamma = {:.10}", panoc_engine.cache.gamma); - println!("> tau = {:.3}", panoc_engine.cache.tau); + println!("> fpr = {:.2?}", panoc_engine.cache.norm_gamma_fpr); + println!( + "> L = {:.3?}", + panoc_engine.cache.lipschitz_constant + ); + println!("> gamma = {:.10?}", panoc_engine.cache.gamma); + println!("> tau = {:.3?}", panoc_engine.cache.tau); println!("> lbfgs dir = {:.11?}", panoc_engine.cache.direction_lbfgs); } diff --git a/src/core/panoc/tests_f32.rs b/src/core/panoc/tests_f32.rs new file mode 100644 index 00000000..5e3df846 --- /dev/null +++ b/src/core/panoc/tests_f32.rs @@ -0,0 +1,222 @@ +use crate::core::panoc::panoc_engine::PANOCEngine; +use crate::core::panoc::*; +use crate::core::*; +use crate::{mocks, FunctionCallResult}; + +const N_DIM: usize = 2; +#[test] +fn t_panoc_init() { + let radius = 0.2_f32; + let ball = constraints::Ball2::new(None, radius); + let problem = Problem::new(&ball, mocks::my_gradient, mocks::my_cost); + let mut panoc_cache = PANOCCache::new(N_DIM, 1e-6, 5); + + { + let mut panoc_engine = PANOCEngine::new(problem, &mut panoc_cache); + let mut u = [0.75, -1.4]; + panoc_engine.init(&mut u).unwrap(); + assert!(2.549_509_967_743_775 > panoc_engine.cache.lipschitz_constant); + assert!(0.372_620_625_931_781 < panoc_engine.cache.gamma, "gamma"); + println!("----------- {} ", panoc_engine.cache.cost_value); + unit_test_utils::assert_nearly_equal( + 6.34125, + panoc_engine.cache.cost_value, + 1e-4, + 1e-10, + "cost value", + ); + unit_test_utils::assert_nearly_equal_array( + &[0.35, -3.05], + &panoc_engine.cache.gradient_u, + 1e-4, + 1e-10, + "gradient at u", + ); + } + println!("cache = {:#?}", &panoc_cache); +} + +fn print_panoc_engine( + panoc_engine: &PANOCEngine, +) where + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat + std::fmt::Debug, +{ + println!("> fpr = {:?}", &panoc_engine.cache.gamma_fpr); + println!("> fpr = {:.2?}", panoc_engine.cache.norm_gamma_fpr); + println!( + "> L = {:.3?}", + panoc_engine.cache.lipschitz_constant + ); + println!("> gamma = {:.10?}", panoc_engine.cache.gamma); + println!("> tau = {:.3?}", panoc_engine.cache.tau); + println!("> lbfgs dir = {:.11?}", panoc_engine.cache.direction_lbfgs); +} + +#[test] +fn t_test_panoc_basic() { + let bounds = constraints::Ball2::new(None, 0.2_f32); + let problem = Problem::new(&bounds, mocks::my_gradient, mocks::my_cost); + let tolerance = 1e-9; + let mut panoc_cache = PANOCCache::new(2, tolerance, 5); + let mut panoc_engine = PANOCEngine::new(problem, &mut panoc_cache); + + let mut u = [0.0, 0.0]; + panoc_engine.init(&mut u).unwrap(); + panoc_engine.step(&mut u).unwrap(); + let fpr0 = panoc_engine.cache.norm_gamma_fpr; + println!("fpr0 = {}", fpr0); + + for i in 1..=100 { + println!("----------------------------------------------------"); + println!("> iter = {}", i); + print_panoc_engine(&panoc_engine); + println!("> u = {:.14?}", u); + if panoc_engine.step(&mut u) != Ok(true) { + break; + } + } + println!("final |fpr| = {}", panoc_engine.cache.norm_gamma_fpr); + assert!(panoc_engine.cache.norm_gamma_fpr <= tolerance); + unit_test_utils::assert_nearly_equal_array(&u, &mocks::SOLUTION_A_F32, 1e-6, 1e-8, ""); +} + +#[test] +fn t_test_panoc_hard() { + let radius: f32 = 0.05_f32; + let bounds = constraints::Ball2::new(None, radius); + let problem = Problem::new( + &bounds, + mocks::hard_quadratic_gradient, + mocks::hard_quadratic_cost, + ); + let n: usize = 3; + let lbfgs_memory: usize = 10; + let tolerance_fpr: f32 = 1e-12; + let mut panoc_cache = PANOCCache::new(n, tolerance_fpr, lbfgs_memory); + let mut panoc_engine = PANOCEngine::new(problem, &mut panoc_cache); + + let mut u = [-20.0, 10., 0.2]; + panoc_engine.init(&mut u).unwrap(); + + println!("L = {}", panoc_engine.cache.lipschitz_constant); + println!("gamma = {}", panoc_engine.cache.gamma); + println!("sigma = {}", panoc_engine.cache.sigma); + + let mut i = 1; + println!("\n*** ITERATION 1"); + while panoc_engine.step(&mut u) == Ok(true) && i < 100 { + i += 1; + println!("+ u_plus = {:?}", u); + println!("\n*** ITERATION {:3}", i); + } + + println!("\nsol = {:?}", u); + assert!(panoc_engine.cache.norm_gamma_fpr <= tolerance_fpr); + unit_test_utils::assert_nearly_equal_array(&u, &mocks::SOLUTION_HARD_F32, 1e-6, 1e-8, ""); +} + +#[test] +fn t_test_panoc_rosenbrock() { + let tolerance = 1e-12_f32; + let a_param = 1.0; + let b_param = 100.0; + let cost_gradient = |u: &[f32], grad: &mut [f32]| -> FunctionCallResult { + mocks::rosenbrock_grad(a_param, b_param, u, grad); + Ok(()) + }; + let cost_function = |u: &[f32], c: &mut f32| -> FunctionCallResult { + *c = mocks::rosenbrock_cost(a_param, b_param, u); + Ok(()) + }; + let bounds = constraints::Ball2::new(None, 1.0); + let problem = Problem::new(&bounds, cost_gradient, cost_function); + let mut panoc_cache = PANOCCache::new(2, tolerance, 2).with_cbfgs_parameters(2.0, 1e-6, 1e-12); + let mut panoc_engine = PANOCEngine::new(problem, &mut panoc_cache); + let mut u_solution = [-1.5, 0.9]; + panoc_engine.init(&mut u_solution).unwrap(); + let mut idx = 1; + while panoc_engine.step(&mut u_solution) == Ok(true) && idx < 50 { + idx += 1; + } + assert!(panoc_engine.cache.norm_gamma_fpr <= tolerance); + println!("u = {:?}", u_solution); +} + +#[test] +fn t_zero_gamma_l() { + let tolerance = 1e-8_f32; + let mut panoc_cache = PANOCCache::new(1, tolerance, 5); + let u = &mut [1e6]; + + // Define the cost function and its gradient. + let df = |u: &[f32], grad: &mut [f32]| -> Result<(), SolverError> { + grad[0] = u[0].signum(); + + Ok(()) + }; + + let f = |u: &[f32], c: &mut f32| -> Result<(), SolverError> { + *c = u[0].abs(); + Ok(()) + }; + + let bounds = constraints::NoConstraints::new(); + + // Problem statement. + let problem = Problem::new(&bounds, df, f); + + let mut panoc_engine = PANOCOptimizer::new(problem, &mut panoc_cache).with_max_iter(100); + + // Invoke the solver. + let _status = panoc_engine.solve(u); + println!("norm_gamma_fpr = {}", panoc_cache.norm_gamma_fpr); + println!("u = {:?}", u); + println!("iters = {}", panoc_cache.iteration); + assert!(panoc_cache.norm_gamma_fpr <= tolerance); +} + +#[test] +fn t_zero_gamma_huber() { + let tolerance = 1e-8_f32; + let mut panoc_cache = PANOCCache::new(1, tolerance, 10); + let u = &mut [1e6]; + let huber_delta = 1e-6; + + // Define the cost function and its gradient. + let df = |u: &[f32], grad: &mut [f32]| -> Result<(), SolverError> { + let u_abs = u[0].abs(); + if u_abs >= huber_delta { + grad[0] = huber_delta * u[0].signum(); + } else { + grad[0] = u[0]; + } + Ok(()) + }; + + let huber_norm = |u: &[f32], y: &mut f32| -> Result<(), SolverError> { + let u_abs = u[0].abs(); + if u_abs <= huber_delta { + *y = 0.5 * u[0].powi(2); + } else { + *y = huber_delta * (u_abs - 0.5 * huber_delta); + } + Ok(()) + }; + + let bounds = constraints::BallInf::new(None, 10000.); + + // Problem statement. + let problem = Problem::new(&bounds, df, huber_norm); + + let mut panoc_engine = PANOCOptimizer::new(problem, &mut panoc_cache).with_max_iter(100); + + // Invoke the solver. + let _status = panoc_engine.solve(u); + println!("norm_gamma_fpr = {}", panoc_cache.norm_gamma_fpr); + println!("u = {:?}", u); + println!("iters = {}", panoc_cache.iteration); + assert!(panoc_cache.norm_gamma_fpr <= tolerance); +} diff --git a/src/core/problem.rs b/src/core/problem.rs index c6e61dae..e227caac 100644 --- a/src/core/problem.rs +++ b/src/core/problem.rs @@ -6,8 +6,8 @@ //! Cost functions are user defined. They can either be defined in Rust or in //! C (and then invoked from Rust via an interface such as icasadi). //! +use crate::core::OptFloat; use crate::{constraints, FunctionCallResult}; - /// Definition of an optimisation problem /// /// The definition of an optimisation problem involves: @@ -15,11 +15,12 @@ use crate::{constraints, FunctionCallResult}; /// - the cost function /// - the set of constraints, which is described by implementations of /// [Constraint](../../panoc_rs/constraints/trait.Constraint.html) -pub struct Problem<'a, GradientType, ConstraintType, CostType> +pub struct Problem<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { /// constraints pub(crate) constraints: &'a ConstraintType, @@ -27,13 +28,17 @@ where pub(crate) gradf: GradientType, /// cost function pub(crate) cost: CostType, + /// phantom data for float type + _phantom: std::marker::PhantomData, } -impl<'a, GradientType, ConstraintType, CostType> Problem<'a, GradientType, ConstraintType, CostType> +impl<'a, GradientType, ConstraintType, CostType, T> + Problem<'a, GradientType, ConstraintType, CostType, T> where - GradientType: Fn(&[f64], &mut [f64]) -> FunctionCallResult, - CostType: Fn(&[f64], &mut f64) -> FunctionCallResult, - ConstraintType: constraints::Constraint, + GradientType: Fn(&[T], &mut [T]) -> FunctionCallResult, + CostType: Fn(&[T], &mut T) -> FunctionCallResult, + ConstraintType: constraints::Constraint, + T: OptFloat, { /// Construct a new instance of an optimisation problem /// @@ -50,11 +55,12 @@ where constraints: &'a ConstraintType, cost_gradient: GradientType, cost: CostType, - ) -> Problem<'a, GradientType, ConstraintType, CostType> { + ) -> Problem<'a, GradientType, ConstraintType, CostType, T> { Problem { constraints, gradf: cost_gradient, cost, + _phantom: std::marker::PhantomData, } } } diff --git a/src/core/solver_status.rs b/src/core/solver_status.rs index 9098d4bf..0ff80852 100644 --- a/src/core/solver_status.rs +++ b/src/core/solver_status.rs @@ -1,16 +1,19 @@ //! Status of the result of a solver (number of iterations, etc) //! //! -use crate::core::ExitStatus; use std::time; +use crate::core::{ExitStatus, OptFloat}; /// Solver status /// /// This structure contais information about the solver status. Instances of /// `SolverStatus` are returned by optimizers. /// #[derive(Debug, PartialEq, Copy, Clone)] -pub struct SolverStatus { +pub struct SolverStatus +where + T: OptFloat, +{ /// exit status of the algorithm exit_status: ExitStatus, /// number of iterations for convergence @@ -18,12 +21,15 @@ pub struct SolverStatus { /// time it took to solve solve_time: time::Duration, /// norm of the fixed-point residual (FPR) - fpr_norm: f64, + fpr_norm: T, /// cost value at the candidate solution - cost_value: f64, + cost_value: T, } -impl SolverStatus { +impl SolverStatus +where + T: OptFloat, +{ /// Constructs a new instance of SolverStatus /// /// ## Arguments @@ -39,9 +45,9 @@ impl SolverStatus { exit_status: ExitStatus, num_iter: usize, solve_time: time::Duration, - fpr_norm: f64, - cost_value: f64, - ) -> SolverStatus { + fpr_norm: T, + cost_value: T, + ) -> SolverStatus { SolverStatus { exit_status, num_iter, @@ -67,12 +73,12 @@ impl SolverStatus { } /// norm of the fixed point residual - pub fn norm_fpr(&self) -> f64 { + pub fn norm_fpr(&self) -> T { self.fpr_norm } /// value of the cost at the solution - pub fn cost_value(&self) -> f64 { + pub fn cost_value(&self) -> T { self.cost_value } diff --git a/src/lib.rs b/src/lib.rs index 2363cfc1..a6ff729b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -58,15 +58,13 @@ pub mod core; pub mod lipschitz_estimator; pub mod matrix_operations; -pub use crate::core::fbs; -pub use crate::core::panoc; -pub use crate::core::{AlgorithmEngine, Optimizer, Problem}; - /* Use Jemalloc if the feature `jem` is activated */ #[cfg(not(target_env = "msvc"))] #[cfg(feature = "jem")] use jemallocator::Jemalloc; +pub use crate::core::{fbs, panoc, AlgorithmEngine, OptFloat, Optimizer, Problem}; + #[cfg(not(target_env = "msvc"))] #[cfg(feature = "jem")] #[global_allocator] diff --git a/src/lipschitz_estimator.rs b/src/lipschitz_estimator.rs index 22d45e43..5fe9fb00 100644 --- a/src/lipschitz_estimator.rs +++ b/src/lipschitz_estimator.rs @@ -37,37 +37,40 @@ //! ``` //! +use crate::core::OptFloat; use crate::{matrix_operations, SolverError}; const DEFAULT_DELTA: f64 = 1e-6; const DEFAULT_EPSILON: f64 = 1e-6; /// Structure for the computation of estimates of the Lipschitz constant of mappings -pub struct LipschitzEstimator<'a, F> +pub struct LipschitzEstimator<'a, F, T> where - F: Fn(&[f64], &mut [f64]) -> Result<(), SolverError>, + F: Fn(&[T], &mut [T]) -> Result<(), SolverError>, + T: OptFloat, { /// `u_decision_var` is the point where the Lipschitz constant is estimated - u_decision_var: &'a mut [f64], + u_decision_var: &'a mut [T], /// internally allocated workspace memory - workspace: Vec, + workspace: Vec, /// `function_value_at_u` a vector which is updated with the /// value of the given function, `F`, at `u`; the provided value /// of `function_value_at_u_p` is not used - function_value_at_u: &'a mut [f64], + function_value_at_u: &'a mut [T], /// /// Function whose Lipschitz constant is to be approximated /// /// For example, in optimization, this is the gradient (Jacobian matrix) /// of the cost function (this is a closure) function: &'a F, - epsilon_lip: f64, - delta_lip: f64, + epsilon_lip: T, + delta_lip: T, } -impl<'a, F> LipschitzEstimator<'a, F> +impl<'a, F, T> LipschitzEstimator<'a, F, T> where - F: Fn(&[f64], &mut [f64]) -> Result<(), SolverError>, + F: Fn(&[T], &mut [T]) -> Result<(), SolverError>, + T: OptFloat, { /// Creates a new instance of this structure /// @@ -88,18 +91,18 @@ where /// /// pub fn new( - u_: &'a mut [f64], + u_: &'a mut [T], f_: &'a F, - function_value_: &'a mut [f64], - ) -> LipschitzEstimator<'a, F> { + function_value_: &'a mut [T], + ) -> LipschitzEstimator<'a, F, T> { let n: usize = u_.len(); LipschitzEstimator { u_decision_var: u_, - workspace: vec![0.0_f64; n], + workspace: vec![T::zero(); n], function_value_at_u: function_value_, function: f_, - epsilon_lip: DEFAULT_EPSILON, - delta_lip: DEFAULT_DELTA, + epsilon_lip: T::from(DEFAULT_EPSILON).unwrap(), + delta_lip: T::from(DEFAULT_DELTA).unwrap(), } } @@ -113,8 +116,8 @@ where /// # Panics /// The method will panic if `delta` is non positive /// - pub fn with_delta(mut self, delta: f64) -> Self { - assert!(delta > 0.0); + pub fn with_delta(mut self, delta: T) -> Self { + assert!(delta > T::zero()); self.delta_lip = delta; self } @@ -129,11 +132,12 @@ where /// # Panics /// The method will panic if `epsilon` is non positive /// - pub fn with_epsilon(mut self, epsilon: f64) -> Self { - assert!(epsilon > 0.0); + pub fn with_epsilon(mut self, epsilon: T) -> Self { + assert!(epsilon > T::zero()); self.epsilon_lip = epsilon; self } + /// /// Getter method for the Jacobian /// @@ -143,7 +147,7 @@ where /// /// If `estimate_local_lipschitz` has not been computed, the result /// will point to a zero vector. - pub fn get_function_value(&self) -> &[f64] { + pub fn get_function_value(&self) -> &[T] { self.function_value_at_u } @@ -180,7 +184,7 @@ where /// No rust-side panics, unless the C function which is called via this interface /// fails. /// - pub fn estimate_local_lipschitz(&mut self) -> Result { + pub fn estimate_local_lipschitz(&mut self) -> Result { // function_value = gradient(u, p) (self.function)(self.u_decision_var, self.function_value_at_u)?; let epsilon_lip = self.epsilon_lip; diff --git a/src/matrix_operations.rs b/src/matrix_operations.rs index f1eafe08..cea79a6c 100644 --- a/src/matrix_operations.rs +++ b/src/matrix_operations.rs @@ -36,10 +36,11 @@ //! ``` //! -use num::{Float, Zero}; use std::iter::Sum; use std::ops::Mul; +use num::{Float, Zero}; + /// Calculate the inner product of two vectors #[inline(always)] pub fn inner_product(a: &[T], b: &[T]) -> T @@ -74,10 +75,10 @@ where #[inline(always)] pub fn norm2_squared_diff(a: &[T], b: &[T]) -> T where - T: Float + Sum + Mul + std::ops::AddAssign, + T: Float + Sum + Mul, { a.iter().zip(b.iter()).fold(T::zero(), |mut sum, (&x, &y)| { - sum += (x - y).powi(2); + sum = sum + (x - y).powi(2); sum }) } diff --git a/src/mocks.rs b/src/mocks.rs index a051da21..640b3086 100644 --- a/src/mocks.rs +++ b/src/mocks.rs @@ -1,81 +1,105 @@ -use crate::{matrix_operations, SolverError}; +use crate::{matrix_operations, OptFloat, SolverError}; pub const SOLUTION_A: [f64; 2] = [-0.148_959_718_255_77, 0.133_457_867_273_39]; +pub const SOLUTION_A_F32: [f32; 2] = [-0.148_959_7, 0.133_457_9]; + pub const SOLUTION_HARD: [f64; 3] = [ -0.041_123_164_672_281, -0.028_440_417_469_206, 0.000_167_276_757_790, ]; +pub const SOLUTION_HARD_F32: [f32; 3] = [-0.041_123_1, -0.028_440_52, 0.000_167_275_72]; -pub fn lipschitz_mock(u: &[f64], g: &mut [f64]) -> Result<(), SolverError> { - g[0] = 3.0 * u[0]; - g[1] = 2.0 * u[1]; - g[2] = 4.5; +pub fn lipschitz_mock(u: &[T], g: &mut [T]) -> Result<(), SolverError> { + g[0] = T::from(3.0).unwrap() * u[0]; + g[1] = T::from(2.0).unwrap() * u[1]; + g[2] = T::from(4.5).unwrap(); Ok(()) } -pub fn void_parameteric_cost(_u: &[f64], _p: &[f64], _cost: &mut f64) -> Result<(), SolverError> { +pub fn void_parameteric_cost( + _u: &[T], + _p: &[T], + _cost: &mut T, +) -> Result<(), SolverError> { Ok(()) } -pub fn void_parameteric_gradient( - _u: &[f64], - _p: &[f64], - _grad: &mut [f64], +pub fn void_parameteric_gradient( + _u: &[T], + _p: &[T], + _grad: &mut [T], ) -> Result<(), SolverError> { Ok(()) } -pub fn void_mapping(_u: &[f64], _result: &mut [f64]) -> Result<(), SolverError> { +pub fn void_mapping(_u: &[T], _result: &mut [T]) -> Result<(), SolverError> { Ok(()) } -pub fn void_cost(_u: &[f64], _cost: &mut f64) -> Result<(), SolverError> { +pub fn void_cost(_u: &[T], _cost: &mut T) -> Result<(), SolverError> { Ok(()) } -pub fn void_gradient(_u: &[f64], _grad: &mut [f64]) -> Result<(), SolverError> { +pub fn void_gradient(_u: &[T], _grad: &mut [T]) -> Result<(), SolverError> { Ok(()) } -pub fn my_cost(u: &[f64], cost: &mut f64) -> Result<(), SolverError> { - *cost = 0.5 * (u[0].powi(2) + 2. * u[1].powi(2) + 2.0 * u[0] * u[1]) + u[0] - u[1] + 3.0; +pub fn my_cost(u: &[T], cost: &mut T) -> Result<(), SolverError> { + *cost = T::from(0.5).unwrap() + * (u[0].powi(2) + + T::from(2.0).unwrap() * u[1].powi(2) + + T::from(2.0).unwrap() * u[0] * u[1]) + + u[0] + - u[1] + + T::from(3.0).unwrap(); Ok(()) } -pub fn my_gradient(u: &[f64], grad: &mut [f64]) -> Result<(), SolverError> { - grad[0] = u[0] + u[1] + 1.0; - grad[1] = u[0] + 2. * u[1] - 1.0; +pub fn my_gradient(u: &[T], grad: &mut [T]) -> Result<(), SolverError> { + grad[0] = u[0] + u[1] + T::one(); + grad[1] = u[0] + T::from(2.0).unwrap() * u[1] - T::one(); Ok(()) } -pub fn rosenbrock_cost(a: f64, b: f64, u: &[f64]) -> f64 { +pub fn rosenbrock_cost(a: T, b: T, u: &[T]) -> T { (a - u[0]).powi(2) + b * (u[1] - u[0].powi(2)).powi(2) } -pub fn rosenbrock_grad(a: f64, b: f64, u: &[f64], grad: &mut [f64]) { - grad[0] = 2.0 * u[0] - 2.0 * a - 4.0 * b * u[0] * (-u[0].powi(2) + u[1]); - grad[1] = b * (-2.0 * u[0].powi(2) + 2.0 * u[1]); +pub fn rosenbrock_grad(a: T, b: T, u: &[T], grad: &mut [T]) { + grad[0] = T::from(2.0).unwrap() * u[0] + - T::from(2.0).unwrap() * a + - T::from(4.0).unwrap() * b * u[0] * (-u[0].powi(2) + u[1]); + grad[1] = b * (T::from(-2.0).unwrap() * u[0].powi(2) + T::from(2.0).unwrap() * u[1]); } -pub fn hard_quadratic_cost(u: &[f64], cost: &mut f64) -> Result<(), SolverError> { - *cost = (4. * u[0].powi(2)) / 2. - + 5.5 * u[1].powi(2) - + 500.5 * u[2].powi(2) - + 5. * u[0] * u[1] - + 25. * u[0] * u[2] - + 5. * u[1] * u[2] +pub fn hard_quadratic_cost(u: &[T], cost: &mut T) -> Result<(), SolverError> { + *cost = T::from(4.0).unwrap() * u[0].powi(2) / T::from(2.0).unwrap() + + T::from(5.5).unwrap() * u[1].powi(2) + + T::from(500.5).unwrap() * u[2].powi(2) + + T::from(5.0).unwrap() * u[0] * u[1] + + T::from(25.0).unwrap() * u[0] * u[2] + + T::from(5.0).unwrap() * u[1] * u[2] + u[0] + u[1] + u[2]; Ok(()) } -pub fn hard_quadratic_gradient(u: &[f64], grad: &mut [f64]) -> Result<(), SolverError> { +pub fn hard_quadratic_gradient(u: &[T], grad: &mut [T]) -> Result<(), SolverError> { // norm(Hessian) = 1000.653 (Lipschitz gradient) - grad[0] = 4. * u[0] + 5. * u[1] + 25. * u[2] + 1.; - grad[1] = 5. * u[0] + 11. * u[1] + 5. * u[2] + 1.; - grad[2] = 25. * u[0] + 5. * u[1] + 1001. * u[2] + 1.; + grad[0] = T::from(4.0).unwrap() * u[0] + + T::from(5.0).unwrap() * u[1] + + T::from(25.0).unwrap() * u[2] + + T::one(); + grad[1] = T::from(5.0).unwrap() * u[0] + + T::from(11.0).unwrap() * u[1] + + T::from(5.0).unwrap() * u[2] + + T::one(); + grad[2] = T::from(25.0).unwrap() * u[0] + + T::from(5.0).unwrap() * u[1] + + T::from(1001.0).unwrap() * u[2] + + T::one(); Ok(()) } @@ -85,28 +109,33 @@ pub fn hard_quadratic_gradient(u: &[f64], grad: &mut [f64]) -> Result<(), Solver /// /// where `m` is the length of `xi`. It is assumed that the length of /// `u` is larger than the length of `xi` -pub fn psi_cost_dummy(u: &[f64], xi: &[f64], cost: &mut f64) -> Result<(), SolverError> { +pub fn psi_cost_dummy(u: &[T], xi: &[T], cost: &mut T) -> Result<(), SolverError> { let u_len = u.len(); let xi_len = xi.len(); assert!(u_len > xi_len); - let sum_u = u.iter().fold(0.0, |mut sum, ui| { - sum += ui; + let sum_u = u.iter().fold(T::zero(), |mut sum, ui| { + sum = sum + *ui; sum }); // psi_cost = 0.5*SUM(ui^2) + xi[0] * sum_u - *cost = - 0.5 * u.iter().fold(0.0, |mut sum_of_squares, ui| { - sum_of_squares += ui.powi(2); + *cost = T::from(0.5).unwrap() + * u.iter().fold(T::zero(), |mut sum_of_squares, ui| { + sum_of_squares = sum_of_squares + ui.powi(2); sum_of_squares - }) + xi[0] * sum_u; + }) + + xi[0] * sum_u; // psi_cost += xi[1..m]'*u[0..m-1] let m = std::cmp::min(u_len, xi_len - 1); - *cost += matrix_operations::inner_product(&u[..m], &xi[1..=m]); + *cost = *cost + matrix_operations::inner_product(&u[..m], &xi[1..=m]); Ok(()) } /// Gradient of `psi_cost` -pub fn psi_gradient_dummy(u: &[f64], xi: &[f64], grad: &mut [f64]) -> Result<(), SolverError> { +pub fn psi_gradient_dummy( + u: &[T], + xi: &[T], + grad: &mut [T], +) -> Result<(), SolverError> { let u_len = u.len(); let xi_len = xi.len(); assert!( @@ -115,11 +144,11 @@ pub fn psi_gradient_dummy(u: &[f64], xi: &[f64], grad: &mut [f64]) -> Result<(), ); assert!(u_len == grad.len(), "u and grad must have equal lengths"); grad.copy_from_slice(u); - grad.iter_mut().for_each(|grad_i| *grad_i += xi[0]); + grad.iter_mut().for_each(|grad_i| *grad_i = *grad_i + xi[0]); xi[1..] .iter() .zip(grad.iter_mut()) - .for_each(|(xi_i, grad_i)| *grad_i += xi_i); + .for_each(|(xi_i, grad_i)| *grad_i = *grad_i + *xi_i); Ok(()) } @@ -132,11 +161,11 @@ pub fn psi_gradient_dummy(u: &[f64], xi: &[f64], grad: &mut [f64]) -> Result<(), /// /// It is `F1: R^3 --> R^2` /// -pub fn mapping_f1_affine(u: &[f64], f1u: &mut [f64]) -> Result<(), SolverError> { +pub fn mapping_f1_affine(u: &[T], f1u: &mut [T]) -> Result<(), SolverError> { assert!(u.len() == 3, "the length of u must be equal to 3"); assert!(f1u.len() == 2, "the length of F1(u) must be equal to 2"); - f1u[0] = 2.0 * u[0] + u[2] - 1.0; - f1u[1] = u[0] + 3.0 * u[1]; + f1u[0] = T::from(2.0).unwrap() * u[0] + u[2] - T::one(); + f1u[1] = u[0] + T::from(3.0).unwrap() * u[1]; Ok(()) } @@ -151,15 +180,15 @@ pub fn mapping_f1_affine(u: &[f64], f1u: &mut [f64]) -> Result<(), SolverError> /// d1 ] /// ``` /// -pub fn mapping_f1_affine_jacobian_product( - _u: &[f64], - d: &[f64], - res: &mut [f64], +pub fn mapping_f1_affine_jacobian_product( + _u: &[T], + d: &[T], + res: &mut [T], ) -> Result<(), SolverError> { assert!(d.len() == 2, "the length of d must be equal to 3"); assert!(res.len() == 3, "the length of res must be equal to 3"); - res[0] = 2.0 * d[0] + d[1]; - res[1] = 3.0 * d[1]; + res[0] = T::from(2.0).unwrap() * d[0] + d[1]; + res[1] = T::from(3.0).unwrap() * d[1]; res[2] = d[0]; Ok(()) } @@ -169,15 +198,15 @@ pub fn mapping_f1_affine_jacobian_product( /// ``` /// f0(u) = 0.5*u'*u + 1'*u /// ``` -pub fn f0(u: &[f64], cost: &mut f64) -> Result<(), SolverError> { - *cost = 0.5 * matrix_operations::norm2_squared(u) + matrix_operations::sum(u); +pub fn f0(u: &[T], cost: &mut T) -> Result<(), SolverError> { + *cost = T::from(0.5).unwrap() * matrix_operations::norm2_squared(u) + matrix_operations::sum(u); Ok(()) } -pub fn d_f0(u: &[f64], grad: &mut [f64]) -> Result<(), SolverError> { +pub fn d_f0(u: &[T], grad: &mut [T]) -> Result<(), SolverError> { grad.iter_mut() .zip(u.iter()) - .for_each(|(grad_i, u_i)| *grad_i = u_i + 1.0); + .for_each(|(grad_i, u_i)| *grad_i = *u_i + T::one()); Ok(()) } diff --git a/src/tests.rs b/src/tests.rs index b8c5e1ca..59d93160 100644 --- a/src/tests.rs +++ b/src/tests.rs @@ -1,9 +1,12 @@ -use crate::{constraints::*, core::fbs::*, core::*}; use std::num::NonZeroUsize; +use crate::constraints::*; +use crate::core::fbs::*; +use crate::core::*; + #[test] fn t_access() { - let radius = 0.2; + let radius = 0.2_f64; let box_constraints = Ball2::new(None, radius); let problem = Problem::new( &box_constraints, @@ -24,3 +27,27 @@ fn t_access() { assert!((-0.14896 - u[0]).abs() < 1e-4); assert!((0.13346 - u[1]).abs() < 1e-4); } + +#[test] +fn t_access_f32() { + let radius = 0.2f32; + let box_constraints = Ball2::new(None, radius); + let problem = Problem::new( + &box_constraints, + super::mocks::my_gradient, + super::mocks::my_cost, + ); + let gamma = 0.1f32; + let tolerance = 1e-6f32; + + let mut fbs_cache = FBSCache::new(NonZeroUsize::new(2).unwrap(), gamma, tolerance); + let mut u = [0.0f32; 2]; + let mut optimizer = FBSOptimizer::new(problem, &mut fbs_cache); + + let status = optimizer.solve(&mut u).unwrap(); + + assert!(status.has_converged()); + assert!(status.norm_fpr() < tolerance); + assert!((-0.14896f32 - u[0]).abs() < 1e-4f32); + assert!((0.13346f32 - u[1]).abs() < 1e-4f32); +}