Index: examples/cflp.rs ================================================================== --- examples/cflp.rs +++ examples/cflp.rs @@ -18,11 +18,10 @@ #![allow(non_upper_case_globals)] //! Example implementation for a capacitated facility location problem. use better_panic; -use crossbeam::channel::unbounded as channel; use log::{info, Level}; use rustop::opts; use std::error::Error; use std::fmt::Write; use std::io::Write as _; @@ -30,16 +29,13 @@ use env_logger::{self, fmt::Color}; use ordered_float::NotNan; use threadpool::ThreadPool; -use bundle::parallel::{ - DefaultSolver as ParallelSolver, EvalResult, FirstOrderProblem as ParallelProblem, - NoBundleSolver as NoParallelSolver, ResultSender, -}; +use bundle::problem::{EvalResult, FirstOrderProblem as ParallelProblem, ResultSender}; +use bundle::solver::sync::{DefaultSolver, NoBundleSolver}; use bundle::{dvec, DVector, Minorant, Real}; -use bundle::{DefaultSolver, FirstOrderProblem, NoBundleSolver, SimpleEvaluation}; const Nfac: usize = 3; const Ncus: usize = 5; const F: [Real; Nfac] = [1000.0, 1000.0, 1000.0]; const CAP: [Real; Nfac] = [500.0, 500.0, 500.0]; @@ -51,18 +47,16 @@ const DEMAND: [Real; 5] = [80.0, 270.0, 250.0, 160.0, 180.0]; #[derive(Debug)] enum EvalError { Customer(usize), - NoObjective(usize), } impl std::fmt::Display for EvalError { fn fmt(&self, fmt: &mut std::fmt::Formatter) -> Result<(), std::fmt::Error> { match self { EvalError::Customer(i) => writeln!(fmt, "Customer subproblem {} failed", i), - EvalError::NoObjective(i) => writeln!(fmt, "No objective value generated for subproblem {}", i), } } } impl Error for EvalError {} @@ -123,56 +117,10 @@ }, primal, )) } -impl FirstOrderProblem for CFLProblem { - type Err = EvalError; - - type Primal = DVector; - - type EvalResult = SimpleEvaluation; - - fn num_variables(&self) -> usize { - Nfac - } - - fn lower_bounds(&self) -> Option> { - Some(vec![0.0; FirstOrderProblem::num_variables(self)]) - } - - fn num_subproblems(&self) -> usize { - Nfac + Ncus - } - - fn evaluate( - &mut self, - index: usize, - lambda: &[Real], - _nullstep_bound: Real, - _relprec: Real, - ) -> Result { - let (tx, rx) = channel(); - ParallelProblem::evaluate(self, index, Arc::new(lambda.iter().cloned().collect()), index, tx)?; - let mut objective = None; - let mut minorants = vec![]; - - for r in rx { - match r { - Ok(EvalResult::ObjectiveValue { value, .. }) => objective = Some(value), - Ok(EvalResult::Minorant { minorant, primal, .. }) => minorants.push((minorant, primal)), - _ => break, - } - } - - Ok(SimpleEvaluation { - objective: objective.ok_or(EvalError::NoObjective(index))?, - minorants, - }) - } -} - impl ParallelProblem for CFLProblem { type Err = EvalError; type Primal = DVector; @@ -302,33 +250,21 @@ synopsis "Solver a simple capacitated facility location problem"; opt minimal:bool, desc:"Use the minimal master model"; } .parse_or_exit(); if !args.minimal { - let mut slv = DefaultSolver::new(CFLProblem::new())?; - slv.params.max_bundle_size = 5; - slv.terminator.termination_precision = 1e-9; - slv.solve()?; - show_primals(|i| slv.aggregated_primals(i))?; - - let mut slv = ParallelSolver::<_>::new(CFLProblem::new()); + let mut slv = DefaultSolver::<_>::new(CFLProblem::new()); slv.terminator.termination_precision = 1e-9; slv.master.max_bundle_size = 5; slv.solve()?; show_primals(|i| slv.aggregated_primal(i).unwrap())?; } else { - let mut slv = NoBundleSolver::new(CFLProblem::new())?; - slv.params.max_bundle_size = 2; - slv.terminator.termination_precision = 1e-5; - slv.solve()?; - show_primals(|i| slv.aggregated_primals(i))?; - - let mut slv = NoParallelSolver::<_>::new(CFLProblem::new()); + let mut slv = NoBundleSolver::<_>::new(CFLProblem::new()); slv.terminator.termination_precision = 1e-5; slv.solve()?; show_primals(|i| slv.aggregated_primal(i).unwrap())?; } Ok(()) } Index: examples/mmcf.rs ================================================================== --- examples/mmcf.rs +++ examples/mmcf.rs @@ -19,45 +19,24 @@ use env_logger::fmt::Color; use log::{info, Level}; use rustop::opts; use std::io::Write; -use bundle::master::{Builder as MasterBuilder, MasterProblem}; +use bundle::master::{Builder, FullMasterBuilder, MasterProblem, MinimalMasterBuilder}; use bundle::mcf::MMCFProblem; -use bundle::parallel; +use bundle::solver::sync::Solver; use bundle::{terminator::StandardTerminator, weighter::HKWeighter}; -use bundle::{DefaultSolver, FirstOrderProblem, NoBundleSolver, Solver}; -use bundle::{FullMasterBuilder, MinimalMasterBuilder}; use std::error::Error; use std::result::Result; -fn solve_standard(mut slv: Solver) -> Result<(), Box> -where - M: MasterBuilder + Default, - M::MasterProblem: MasterProblem, -{ - slv.weighter.set_weight_bounds(1e-1, 100.0); - slv.terminator.termination_precision = 1e-6; - slv.solve()?; - - let costs: f64 = (0..slv.problem().num_subproblems()) - .map(|i| { - let aggr_primals = slv.aggregated_primals(i); - slv.problem().get_primal_costs(i, &aggr_primals) - }) - .sum(); - info!("Primal costs: {}", costs); - Ok(()) -} - -fn solve_parallel(master: M, mmcf: MMCFProblem) -> Result<(), Box> -where - M: MasterBuilder, - M::MasterProblem: MasterProblem, -{ - let mut slv = parallel::Solver::<_, StandardTerminator, HKWeighter, M>::with_master(mmcf, master); +fn solve_parallel(master: M, mmcf: MMCFProblem) -> Result<(), Box> +where + M: Builder, + M::MasterProblem: MasterProblem, +{ + let mut slv = Solver::<_, StandardTerminator, HKWeighter, M>::with_master(mmcf, master); slv.weighter.set_weight_bounds(1e-1, 100.0); slv.terminator.termination_precision = 1e-6; slv.solve()?; let costs: f64 = (0..slv.problem().num_subproblems()) @@ -105,66 +84,28 @@ let filename = args.file; info!("Reading instance: {}", filename); if !args.minimal { - { - let mut mmcf = MMCFProblem::read_mnetgen(&filename)?; - if args.aggregated { - mmcf.multimodel = false; - mmcf.set_separate_constraints(false); - } else { - mmcf.multimodel = true; - mmcf.set_separate_constraints(args.separate); - } - - let mut solver = DefaultSolver::new(mmcf)?; - solver.params.max_bundle_size = if args.bundle_size <= 1 { - if args.aggregated { - 50 - } else { - 5 - } - } else { - args.bundle_size - }; - solve_standard(solver)?; - } - - println!("---------------------------------"); - { - let mut mmcf = MMCFProblem::read_mnetgen(&filename)?; - mmcf.set_separate_constraints(args.separate); - mmcf.multimodel = true; - - let mut master = FullMasterBuilder::default(); - if args.aggregated { - master.max_bundle_size(if args.bundle_size <= 1 { 50 } else { args.bundle_size }); - master.use_full_aggregation(); - } else { - master.max_bundle_size(if args.bundle_size <= 1 { 5 } else { args.bundle_size }); - } - solve_parallel(master, mmcf)?; - } - } else { - { - let mut mmcf = MMCFProblem::read_mnetgen(&filename)?; - mmcf.multimodel = false; - - let mut solver = NoBundleSolver::new(mmcf)?; - solver.params.max_bundle_size = 2; - solve_standard(solver)?; - } - - println!("---------------------------------"); - { - let mut mmcf = MMCFProblem::read_mnetgen(&filename)?; - mmcf.set_separate_constraints(args.separate); - mmcf.multimodel = true; - - let master = MinimalMasterBuilder::default(); - solve_parallel(master, mmcf)?; - } + let mut mmcf = MMCFProblem::read_mnetgen(&filename)?; + mmcf.set_separate_constraints(args.separate); + mmcf.multimodel = true; + + let mut master = FullMasterBuilder::default(); + if args.aggregated { + master.max_bundle_size(if args.bundle_size <= 1 { 50 } else { args.bundle_size }); + master.use_full_aggregation(); + } else { + master.max_bundle_size(if args.bundle_size <= 1 { 5 } else { args.bundle_size }); + } + solve_parallel(master, mmcf)?; + } else { + let mut mmcf = MMCFProblem::read_mnetgen(&filename)?; + mmcf.set_separate_constraints(args.separate); + mmcf.multimodel = true; + + let master = MinimalMasterBuilder::default(); + solve_parallel(master, mmcf)?; } Ok(()) } Index: examples/quadratic.rs ================================================================== --- examples/quadratic.rs +++ examples/quadratic.rs @@ -25,15 +25,13 @@ use rustop::opts; use std::io::Write; use std::sync::Arc; use std::thread; -use bundle::parallel::{ - DefaultSolver as ParallelSolver, EvalResult, FirstOrderProblem as ParallelProblem, - NoBundleSolver as NoParallelSolver, ResultSender, -}; -use bundle::{DVector, DefaultSolver, FirstOrderProblem, Minorant, NoBundleSolver, Real, SimpleEvaluation}; +use bundle::problem::{EvalResult, FirstOrderProblem as ParallelProblem, ResultSender}; +use bundle::solver::sync::{DefaultSolver, NoBundleSolver}; +use bundle::{DVector, Minorant, Real}; #[derive(Clone)] struct QuadraticProblem { a: [[Real; 2]; 2], b: [Real; 2], @@ -48,54 +46,10 @@ c: 3.0, } } } -impl FirstOrderProblem for QuadraticProblem { - type Err = Box; - type Primal = (); - type EvalResult = SimpleEvaluation<()>; - - fn num_variables(&self) -> usize { - 2 - } - - #[allow(unused_variables)] - fn evaluate( - &mut self, - fidx: usize, - x: &[Real], - nullstep_bnd: Real, - relprec: Real, - ) -> Result { - assert_eq!(fidx, 0); - let mut objective = self.c; - let mut g = dvec![0.0; 2]; - - for i in 0..2 { - g[i] += (0..2).map(|j| self.a[i][j] * x[j]).sum::(); - objective += x[i] * (g[i] + self.b[i]); - g[i] = 2.0 * g[i] + self.b[i]; - } - - debug!("Evaluation at {:?}", x); - debug!(" objective={}", objective); - debug!(" subgradient={}", g); - - Ok(SimpleEvaluation { - objective: objective, - minorants: vec![( - Minorant { - constant: objective, - linear: g, - }, - (), - )], - }) - } -} - impl ParallelProblem for QuadraticProblem { type Err = Box; type Primal = (); fn num_variables(&self) -> usize { @@ -110,37 +64,49 @@ fn stop(&mut self) {} fn evaluate( &mut self, - i: usize, - y: Arc, + fidx: usize, + x: Arc, index: I, tx: ResultSender, ) -> Result<(), Self::Err> where I: Send + Copy + 'static, { - let y = y.clone(); - let mut p = self.clone(); - thread::spawn(move || match FirstOrderProblem::evaluate(&mut p, i, &y, 0.0, 0.0) { - Ok(res) => { - tx.send(Ok(EvalResult::ObjectiveValue { - index, - value: res.objective, - })) - .unwrap(); - for (minorant, primal) in res.minorants { - tx.send(Ok(EvalResult::Minorant { - index, - minorant, - primal, - })) - .unwrap(); - } - } - Err(err) => tx.send(Err(err)).unwrap(), + let x = x.clone(); + let p = self.clone(); + thread::spawn(move || { + assert_eq!(fidx, 0); + let mut objective = p.c; + let mut g = dvec![0.0; 2]; + + for i in 0..2 { + g[i] += (0..2).map(|j| p.a[i][j] * x[j]).sum::(); + objective += x[i] * (g[i] + p.b[i]); + g[i] = 2.0 * g[i] + p.b[i]; + } + + debug!("Evaluation at {:?}", x); + debug!(" objective={}", objective); + debug!(" subgradient={}", g); + + tx.send(Ok(EvalResult::ObjectiveValue { + index, + value: objective, + })) + .unwrap(); + tx.send(Ok(EvalResult::Minorant { + index, + minorant: Minorant { + constant: objective, + linear: g, + }, + primal: (), + })) + .unwrap(); }); Ok(()) } } @@ -171,34 +137,16 @@ synopsis "Solver a simple quadratic optimization problem"; opt minimal:bool, desc:"Use the minimal master model"; } .parse_or_exit(); - { - let f = QuadraticProblem::new(); - if !args.minimal { - let mut solver = DefaultSolver::new(f).map_err(|e| format!("{}", e))?; - solver.weighter.set_weight_bounds(1.0, 1.0); - solver.solve().map_err(|e| format!("{}", e))?; - } else { - let mut solver = NoBundleSolver::new(f).map_err(|e| format!("{}", e))?; - solver.params.max_bundle_size = 2; - solver.weighter.set_weight_bounds(1.0, 1.0); - solver.solve().map_err(|e| format!("{}", e))?; - } - } - - println!("-------------------------"); - - { - let f = QuadraticProblem::new(); - if !args.minimal { - let mut solver = ParallelSolver::<_>::new(f); - solver.solve().map_err(|e| format!("{}", e))?; - } else { - let mut solver = NoParallelSolver::<_>::new(f); - solver.solve().map_err(|e| format!("{}", e))?; - } + let f = QuadraticProblem::new(); + if !args.minimal { + let mut solver = DefaultSolver::<_>::new(f); + solver.solve().map_err(|e| format!("{}", e))?; + } else { + let mut solver = NoBundleSolver::<_>::new(f); + solver.solve().map_err(|e| format!("{}", e))?; } Ok(()) } ADDED src/data/aggregatable.rs Index: src/data/aggregatable.rs ================================================================== --- /dev/null +++ src/data/aggregatable.rs @@ -0,0 +1,118 @@ +/* + * Copyright (c) 2019 Frank Fischer + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see + */ + +//! Objects that can be combined linearly. + +use super::Real; +use std::borrow::Borrow; + +/// An aggregatable object. +pub trait Aggregatable: Default { + /// Return a scaled version of `other`, i.e. `alpha * other`. + fn new_scaled(alpha: Real, other: A) -> Self + where + A: Borrow; + + /// Add a scaled version of `other` to `self`. + /// + /// This sets `self = self + alpha * other`. + fn add_scaled(&mut self, alpha: Real, other: A) + where + A: Borrow; + + /// Return $\sum\_{i=1}\^n alpha_i m_i$. + /// + /// If `aggregates` is empty return the default value. + fn combine(aggregates: I) -> Self + where + I: IntoIterator, + A: Borrow, + { + let mut it = aggregates.into_iter(); + let mut x; + if let Some((alpha, y)) = it.next() { + x = Self::new_scaled(alpha, y); + } else { + return Self::default(); + } + + for (alpha, y) in it { + x.add_scaled(alpha, y); + } + + x + } +} + +/// Implement for empty tuples. +impl Aggregatable for () { + fn new_scaled(_alpha: Real, _other: A) -> Self + where + A: Borrow, + { + } + + fn add_scaled(&mut self, _alpha: Real, _other: A) + where + A: Borrow, + { + } +} + +/// Implement for scalar values. +impl Aggregatable for Real { + fn new_scaled(alpha: Real, other: A) -> Self + where + A: Borrow, + { + alpha * other.borrow() + } + + fn add_scaled(&mut self, alpha: Real, other: A) + where + A: Borrow, + { + *self += alpha * other.borrow() + } +} + +/// Implement for vectors of aggregatable objects. +impl Aggregatable for Vec +where + T: Aggregatable, +{ + fn new_scaled(alpha: Real, other: A) -> Self + where + A: std::borrow::Borrow, + { + other + .borrow() + .iter() + .map(|y| Aggregatable::new_scaled(alpha, y)) + .collect() + } + + fn add_scaled(&mut self, alpha: Real, other: A) + where + A: std::borrow::Borrow, + { + debug_assert_eq!(self.len(), other.borrow().len(), "Vectors must have the same size"); + for (ref mut x, y) in self.iter_mut().zip(other.borrow()) { + x.add_scaled(alpha, y) + } + } +} ADDED src/data/minorant.rs Index: src/data/minorant.rs ================================================================== --- /dev/null +++ src/data/minorant.rs @@ -0,0 +1,106 @@ +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer +// +// This program is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see +// + +//! A linear minorant. + +use super::{Aggregatable, DVector, Real}; + +use std::borrow::Borrow; +use std::fmt; + +/// A linear minorant of a convex function. +/// +/// A linear minorant of a convex function $f \colon \mathbb{R}\^n \to +/// \mathbb{R}$ is a linear function of the form +/// +/// \\[ l \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto \langle g, x +/// \rangle + c \\] +/// +/// such that $l(x) \le f(x)$ for all $x \in \mathbb{R}\^n$. +#[derive(Clone, Debug)] +pub struct Minorant { + /// The constant term. + pub constant: Real, + + /// The linear term. + pub linear: DVector, +} + +impl fmt::Display for Minorant { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "{} + y * {}", self.constant, self.linear)?; + Ok(()) + } +} + +impl Default for Minorant { + fn default() -> Minorant { + Minorant { + constant: 0.0, + linear: dvec![], + } + } +} + +impl Minorant { + /// Return a new 0 minorant. + pub fn new(constant: Real, linear: Vec) -> Minorant { + Minorant { + constant, + linear: DVector(linear), + } + } + + /** + * Evaluate minorant at some point. + * + * This function computes $c + \langle g, x \rangle$ for this minorant + * \\[\ell \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto c + \langle g, x \rangle\\] + * and the given point $x \in \mathbb{R}\^n$. + */ + pub fn eval(&self, x: &DVector) -> Real { + self.constant + self.linear.dot(x) + } + + /** + * Move the center of the minorant. + */ + pub fn move_center(&mut self, alpha: Real, d: &DVector) { + self.constant += alpha * self.linear.dot(d); + } +} + +impl Aggregatable for Minorant { + fn new_scaled(alpha: Real, other: A) -> Self + where + A: Borrow, + { + let m = other.borrow(); + Minorant { + constant: alpha * m.constant, + linear: DVector::scaled(&m.linear, alpha), + } + } + + fn add_scaled(&mut self, alpha: Real, other: A) + where + A: Borrow, + { + let m = other.borrow(); + self.constant += alpha * m.constant; + self.linear.add_scaled(alpha, &m.linear); + } +} ADDED src/data/mod.rs Index: src/data/mod.rs ================================================================== --- /dev/null +++ src/data/mod.rs @@ -0,0 +1,30 @@ +/* + * Copyright (c) 2019 Frank Fischer + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see + */ + +//! General types and data structures. + +pub mod vector; +pub use vector::{DVector, Vector}; + +pub mod aggregatable; +pub use aggregatable::Aggregatable; + +pub mod minorant; +pub use minorant::Minorant; + +/// Type used for real numbers throughout the library. +pub type Real = f64; ADDED src/data/vector.rs Index: src/data/vector.rs ================================================================== --- /dev/null +++ src/data/vector.rs @@ -0,0 +1,287 @@ +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer +// +// This program is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see +// + +//! Finite-dimensional sparse and dense vectors. + +use crate::{Aggregatable, Real}; +use std::fmt; +use std::ops::{Deref, DerefMut}; +// use std::cmp::min; +use std::borrow::Borrow; +use std::iter::FromIterator; +use std::vec::IntoIter; + +#[cfg(feature = "blas")] +use {openblas_src as _, rs_blas as blas, std::os::raw::c_int}; + +/// Type of dense vectors. +#[derive(Debug, Clone, PartialEq, Default)] +pub struct DVector(pub Vec); + +impl Deref for DVector { + type Target = Vec; + + fn deref(&self) -> &Vec { + &self.0 + } +} + +impl DerefMut for DVector { + fn deref_mut(&mut self) -> &mut Vec { + &mut self.0 + } +} + +impl fmt::Display for DVector { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + write!(f, "(")?; + for (i, x) in self.iter().enumerate() { + if i > 0 { + write!(f, ", ")?; + } + write!(f, "{}", x)? + } + write!(f, ")")?; + Ok(()) + } +} + +impl FromIterator for DVector { + fn from_iter>(iter: I) -> Self { + DVector(Vec::from_iter(iter)) + } +} + +impl IntoIterator for DVector { + type Item = Real; + type IntoIter = IntoIter; + + fn into_iter(self) -> IntoIter { + self.0.into_iter() + } +} + +/// Type of dense or vectors. +#[derive(Debug, Clone)] +pub enum Vector { + /// A vector with dense storage. + Dense(DVector), + + /** + * A vector with sparse storage. + * + * For each non-zero element this vector stores an index and the + * value of the element in addition to the size of the vector. + */ + Sparse { size: usize, elems: Vec<(usize, Real)> }, +} + +impl fmt::Display for Vector { + fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { + match *self { + Vector::Dense(ref v) => write!(f, "{}", v), + Vector::Sparse { size, ref elems } => { + let mut it = elems.iter(); + write!(f, "{}:(", size)?; + if let Some(&(i, x)) = it.next() { + write!(f, "{}:{}", i, x)?; + for &(i, x) in it { + write!(f, ", {}:{}", i, x)?; + } + } + write!(f, ")") + } + } + } +} + +impl DVector { + /// Set all elements to 0. + pub fn init0(&mut self, size: usize) { + self.clear(); + self.extend((0..size).map(|_| 0.0)); + } + + /// Set self = factor * y. + pub fn scal(&mut self, factor: Real, y: &DVector) { + self.clear(); + self.extend(y.iter().map(|y| factor * y)); + } + + /// Return factor * self. + pub fn scaled(&self, factor: Real) -> DVector { + let mut x = DVector::default(); + x.scal(factor, self); + x + } + + /// Return the inner product with another vector. + pub fn dot(&self, other: &DVector) -> Real { + assert_eq!(self.len(), other.len()); + self.dot_begin(other) + } + + /// Return the inner product with another vector. + /// + /// The inner product is computed on the smaller of the two + /// dimensions. All other elements are assumed to be zero. + pub fn dot_begin(&self, other: &DVector) -> Real { + #[cfg(feature = "blas")] + unsafe { + blas::ddot(self.len().min(other.len()) as c_int, &self, 1, &other, 1) + } + #[cfg(not(feature = "blas"))] + { + self.iter().zip(other.iter()).map(|(x, y)| x * y).sum::() + } + } + + /// Add two vectors and store result in this vector. + pub fn add(&mut self, x: &DVector, y: &DVector) { + assert_eq!(x.len(), y.len()); + self.clear(); + self.extend(x.iter().zip(y.iter()).map(|(a, b)| a + b)); + } + + /// Add two vectors and store result in this vector. + pub fn add_scaled(&mut self, alpha: Real, y: &DVector) { + assert_eq!(self.len(), y.len()); + #[cfg(feature = "blas")] + unsafe { + blas::daxpy(self.len() as c_int, alpha, &y, 1, &mut self[..], 1) + } + #[cfg(not(feature = "blas"))] + { + for (x, y) in self.iter_mut().zip(y.iter()) { + *x += alpha * y; + } + } + } + + /// Add two vectors and store result in this vector. + /// + /// In contrast to `add_scaled`, the two vectors might have + /// different sizes. The size of the resulting vector is the + /// larger of the two vector sizes and the remaining entries of + /// the smaller vector are assumed to be 0.0. + pub fn add_scaled_begin(&mut self, alpha: Real, y: &DVector) { + #[cfg(feature = "blas")] + unsafe { + let n = self.len(); + blas::daxpy(n.min(y.len()) as c_int, alpha, &y, 1, &mut self[..], 1); + } + #[cfg(not(feature = "blas"))] + { + for (x, y) in self.iter_mut().zip(y.iter()) { + *x += alpha * y; + } + } + let n = self.len(); + if n < y.len() { + self.extend(y[n..].iter().map(|y| alpha * y)); + } + } + + /// Return the 2-norm of this vector. + pub fn norm2(&self) -> Real { + #[cfg(feature = "blas")] + unsafe { + blas::dnrm2(self.len() as c_int, &self, 1) + } + #[cfg(not(feature = "blas"))] + { + self.iter().map(|x| x * x).sum::().sqrt() + } + } +} + +impl Aggregatable for DVector { + fn new_scaled(alpha: Real, other: A) -> Self + where + A: Borrow, + { + DVector::scaled(&other.borrow(), alpha) + } + + fn add_scaled(&mut self, alpha: Real, other: A) + where + A: Borrow, + { + DVector::add_scaled(self, alpha, &other.borrow()) + } +} + +impl Vector { + /** + * Return a sparse vector with the given non-zeros. + */ + pub fn new_sparse(n: usize, indices: &[usize], values: &[Real]) -> Vector { + assert_eq!(indices.len(), values.len()); + + if indices.is_empty() { + Vector::Sparse { size: n, elems: vec![] } + } else { + let mut ordered: Vec<_> = (0..n).collect(); + ordered.sort_by_key(|&i| indices[i]); + assert!(*indices.last().unwrap() < n); + let mut elems = Vec::with_capacity(indices.len()); + let mut last_idx = n; + for i in ordered { + let val = unsafe { *values.get_unchecked(i) }; + if val != 0.0 { + let idx = unsafe { *indices.get_unchecked(i) }; + if idx != last_idx { + elems.push((idx, val)); + last_idx = idx; + } else { + elems.last_mut().unwrap().1 += val; + if elems.last_mut().unwrap().1 == 0.0 { + elems.pop(); + last_idx = n; + } + } + } + } + Vector::Sparse { size: n, elems } + } + } + + /** + * Convert vector to a dense vector. + * + * This function always returns a copy of the vector. + */ + pub fn to_dense(&self) -> DVector { + match *self { + Vector::Dense(ref x) => x.clone(), + Vector::Sparse { size: n, elems: ref xs } => { + let mut v = vec![0.0; n]; + for &(i, x) in xs { + unsafe { *v.get_unchecked_mut(i) = x }; + } + DVector(v) + } + } + } +} + +#[test] +fn test_add_scaled_begin() { + let mut x = dvec![1.0; 5]; + let y = dvec![2.0; 7]; + x.add_scaled_begin(3.0, &y); + assert_eq!(x, dvec![7.0, 7.0, 7.0, 7.0, 7.0, 6.0, 6.0]); +} DELETED src/firstorderproblem.rs Index: src/firstorderproblem.rs ================================================================== --- src/firstorderproblem.rs +++ /dev/null @@ -1,180 +0,0 @@ -// Copyright (c) 2016, 2017, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -//! Problem description of a first-order convex optimization problem. - -use crate::solver::UpdateState; -use crate::{Aggregatable, Minorant, Real}; - -use std::result::Result; -use std::vec::IntoIter; - -/** - * Trait for results of an evaluation. - * - * An evaluation returns the function value at the point of evaluation - * and one or more subgradients. - * - * The subgradients (linear minorants) can be obtained by iterating over the result. The - * subgradients are centered around the point of evaluation. - */ -pub trait Evaluation

: IntoIterator { - /// Return the function value at the point of evaluation. - fn objective(&self) -> Real; -} - -/** - * Simple standard evaluation result. - * - * This result consists of the function value and a list of one or - * more minorants and associated primal information. - */ -pub struct SimpleEvaluation

{ - pub objective: Real, - pub minorants: Vec<(Minorant, P)>, -} - -impl

IntoIterator for SimpleEvaluation

{ - type Item = (Minorant, P); - type IntoIter = IntoIter<(Minorant, P)>; - - fn into_iter(self) -> Self::IntoIter { - self.minorants.into_iter() - } -} - -impl

Evaluation

for SimpleEvaluation

{ - fn objective(&self) -> Real { - self.objective - } -} - -/// Problem update information. -/// -/// The solver calls the `update` method of the problem regularly. -/// This method can modify the problem by adding (or removing) -/// variables. The possible updates are encoded in this type. -#[derive(Debug, Clone, Copy)] -pub enum Update { - /// Add a variable with bounds. - /// - /// The initial value of the variable will be the feasible value - /// closest to 0. - AddVariable { lower: Real, upper: Real }, - /// Add a variable with bounds and initial value. - AddVariableValue { lower: Real, upper: Real, value: Real }, - /// Change the current value of a variable. The bounds remain - /// unchanged. - MoveVariable { index: usize, value: Real }, -} - -/** - * Trait for implementing a first-order problem description. - * - */ -pub trait FirstOrderProblem { - /// Error raised by this oracle. - type Err; - - /// The primal information associated with a minorant. - type Primal: Aggregatable; - - /// Custom evaluation result value. - type EvalResult: Evaluation; - - /// Return the number of variables. - fn num_variables(&self) -> usize; - - /** - * Return the lower bounds on the variables. - * - * If no lower bounds a specified, $-\infty$ is assumed. - * - * The lower bounds must be less then or equal the upper bounds. - */ - fn lower_bounds(&self) -> Option> { - None - } - - /** - * Return the upper bounds on the variables. - * - * If no lower bounds a specified, $+\infty$ is assumed. - * - * The upper bounds must be greater than or equal the upper bounds. - */ - fn upper_bounds(&self) -> Option> { - None - } - - /// Return the number of subproblems. - fn num_subproblems(&self) -> usize { - 1 - } - - /** - * Evaluate the i^th subproblem at the given point. - * - * The returned evaluation result must contain (an upper bound on) - * the objective value at $y$ as well as at least one subgradient - * centered at $y$. - * - * If the evaluation process reaches a lower bound on the function - * value at $y$ and this bound is larger than $nullstep_bound$, - * the evaluation may stop and return the lower bound and a - * minorant. In this case the function value is guaranteed to be - * large enough so that the new point is rejected as candidate. - * - * The returned objective value should be an upper bound on the - * true function value within $relprec \cdot (\\|f(y)\\| + 1.0)$, - * otherwise the returned objective should be the maximum of all - * linear minorants at $y$. - * - * Note that `nullstep_bound` and `relprec` are usually only - * useful if there is only a `single` subproblem. - */ - fn evaluate( - &mut self, - i: usize, - y: &[Real], - nullstep_bound: Real, - relprec: Real, - ) -> Result; - - /// Return updates of the problem. - /// - /// The default implementation returns no updates. - fn update(&mut self, _state: &UpdateState) -> Result, Self::Err> { - Ok(vec![]) - } - - /// Return new components for a subgradient. - /// - /// The components are typically generated by some primal information. The - /// corresponding primal along with its subproblem index is passed as a - /// parameter. - /// - /// The default implementation fails because it should never be - /// called. - fn extend_subgradient( - &mut self, - _i: usize, - _primal: &Self::Primal, - _vars: &[usize], - ) -> Result, Self::Err> { - unimplemented!() - } -} Index: src/lib.rs ================================================================== --- src/lib.rs +++ src/lib.rs @@ -14,45 +14,26 @@ // along with this program. If not, see // //! Proximal bundle method implementation. -/// Type used for real numbers throughout the library. -pub type Real = f64; - #[macro_export] macro_rules! dvec { ( $ elem : expr ; $ n : expr ) => { DVector(vec![$elem; $n]) }; ( $ ( $ x : expr ) , * ) => { DVector(vec![$($x),*]) }; ( $ ( $ x : expr , ) * ) => { DVector(vec![$($x,)*]) }; } -pub mod vector; -pub use crate::vector::{DVector, Vector}; - -pub mod minorant; -pub use crate::minorant::{Aggregatable, Minorant}; - -pub mod firstorderproblem; -pub use crate::firstorderproblem::{Evaluation, FirstOrderProblem, SimpleEvaluation, Update}; +mod data; +pub use data::{Aggregatable, DVector, Minorant, Real, Vector}; + +pub mod problem; pub mod solver; -pub use crate::solver::{BundleState, FullMasterBuilder, IterationInfo, Solver, SolverParams, Step, UpdateState}; - -pub mod parallel; pub mod weighter; pub mod terminator; pub mod master; pub mod mcf; - -/// The minimal bundle builder. -pub type MinimalMasterBuilder = master::boxed::Builder; - -/// The default bundle solver with general master problem. -pub type DefaultSolver

= Solver; - -/// A bundle solver with a minimal cutting plane model. -pub type NoBundleSolver

= Solver; Index: src/master/boxed.rs ================================================================== --- src/master/boxed.rs +++ src/master/boxed.rs @@ -12,11 +12,15 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -use crate::master::{self, MasterProblem, SubgradientExtension, UnconstrainedMasterProblem}; +pub mod unconstrained; +use self::unconstrained::UnconstrainedMasterProblem; + +use super::MasterProblem; +pub use super::SubgradientExtension; use crate::{DVector, Minorant, Real}; use itertools::multizip; use log::debug; use std::f64::{EPSILON, INFINITY, NEG_INFINITY}; @@ -376,13 +380,13 @@ /// /// `B` is a builder of the underlying `UnconstrainedMasterProblem`. #[derive(Default)] pub struct Builder(B); -impl master::Builder for Builder +impl super::Builder for Builder where - B: master::unconstrained::Builder, + B: unconstrained::Builder, B::MasterProblem: UnconstrainedMasterProblem, { type MasterProblem = BoxedMasterProblem; fn build(&mut self) -> Result::Err> { ADDED src/master/boxed/unconstrained.rs Index: src/master/boxed/unconstrained.rs ================================================================== --- /dev/null +++ src/master/boxed/unconstrained.rs @@ -0,0 +1,142 @@ +// Copyright (c) 2016, 2017, 2019 Frank Fischer +// +// This program is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see +// + +pub mod cpx; +pub mod minimal; + +use crate::{DVector, Minorant, Real}; + +pub use super::SubgradientExtension; + +use std::error::Error; + +/** + * Trait for master problems without box constraints. + * + * Implementors of this trait are supposed to solve quadratic + * optimization problems of the form + * + * \\[ \min \left\\{ \hat{f}(d) + \frac{u}{2} \\| d \\|\^2 \colon + * d \in \mathbb{R}\^n \right\\}. \\] + * + * where $\hat{f}$ is a piecewise linear model, i.e. + * + * \\[ \hat{f}(d) = \max \\{ \ell_i(d) = c_i + \langle g_i, d \rangle \colon + * i=1,\dotsc,k \\} + * = \max \left\\{ \sum_{i=1}\^k \alpha_i \ell_i(d) \colon + * \alpha \in \Delta \right\\}, \\] + * + * where $\Delta := \left\\{ \alpha \in \mathbb{R}\^k \colon \sum_{i=1}\^k + * \alpha_i = 1 \right\\}$. Note, the unconstrained solver is expected + * to compute *dual* optimal solutions, i.e. the solver must compute + * optimal coefficients $\bar{\alpha}$ for the dual problem + * + * \\[ \max_{\alpha \in \Delta} \min_{d \in \mathbb{R}\^n} + * \sum_{i=1}\^k \alpha_i \ell_i(d) + \frac{u}{2} \\| d \\|\^2. \\] + */ +pub trait UnconstrainedMasterProblem: Send + 'static { + /// Unique index for a minorant. + type MinorantIndex: Copy + Eq; + + /// Error type for this master problem. + type Err: Error + Send + Sync; + + /// Return a new instance of the unconstrained master problem. + fn new() -> Result + where + Self: Sized; + + /// Return the number of subproblems. + fn num_subproblems(&self) -> usize; + + /// Set the number of subproblems (different function models.) + fn set_num_subproblems(&mut self, n: usize) -> Result<(), Self::Err>; + + /// Return the current weight. + fn weight(&self) -> Real; + + /// Set the weight of the quadratic term, must be > 0. + fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err>; + + /// Return the number of minorants of subproblem `fidx`. + fn num_minorants(&self, fidx: usize) -> usize; + + /// Compress the bundle. + /// + /// When some minorants are compressed, the callback is called with the + /// coefficients and indices of the compressed minorants and the index of + /// the new minorant. The callback may be called several times. + fn compress(&mut self, f: F) -> Result<(), Self::Err> + where + F: FnMut(Self::MinorantIndex, &mut dyn Iterator); + + /// Add a new minorant to the model. + fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result; + + /// Add or move some variables. + /// + /// The variables in `changed` have been changed, so the subgradient + /// information must be updated. Furthermore, `nnew` new variables + /// are added. + fn add_vars( + &mut self, + nnew: usize, + changed: &[usize], + extend_subgradient: &mut SubgradientExtension, + ) -> Result<(), Self::Err>; + + /// Solve the master problem. + fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<(), Self::Err>; + + /// Return the current dual optimal solution. + fn dualopt(&self) -> &DVector; + + /// Return the current dual optimal solution value. + fn dualopt_cutval(&self) -> Real; + + /// Return the multiplier associated with a minorant. + fn multiplier(&self, min: Self::MinorantIndex) -> Real; + + /// Return the multipliers associated with a subproblem. + fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box + 'a>; + + /// Return the value of the current model at the given point. + fn eval_model(&self, y: &DVector) -> Real; + + /// Aggregate the given minorants according to the current solution. + /// + /// The (indices of the) minorants to be aggregated get invalid + /// after this operation. The function returns the index of the + /// aggregated minorant along with the coefficients of the convex + /// combination. The index of the new aggregated minorant might or + /// might not be one of indices of the original minorants. + /// + /// # Error + /// The indices of the minorants `mins` must belong to subproblem `fidx`. + fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector), Self::Err>; + + /// Move the center of the master problem along $\alpha \cdot d$. + fn move_center(&mut self, alpha: Real, d: &DVector); +} + +/// A builder for creating unconstrained master problem solvers. +pub trait Builder { + /// The master problem to be build. + type MasterProblem: UnconstrainedMasterProblem; + + /// Create a new master problem. + fn build(&mut self) -> Result::Err>; +} ADDED src/master/boxed/unconstrained/cpx.rs Index: src/master/boxed/unconstrained/cpx.rs ================================================================== --- /dev/null +++ src/master/boxed/unconstrained/cpx.rs @@ -0,0 +1,814 @@ +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer +// +// This program is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see +// + +//! Master problem implementation using CPLEX. + +#![allow(unused_unsafe)] + +use super::{SubgradientExtension, UnconstrainedMasterProblem}; +use crate::{Aggregatable, DVector, Minorant, Real}; + +use c_str_macro::c_str; +use cplex_sys as cpx; +use cplex_sys::trycpx; +use log::debug; + +use std; +use std::f64::NEG_INFINITY; +use std::iter::repeat; +use std::ops::{Deref, DerefMut}; +use std::os::raw::{c_char, c_int}; +use std::ptr; +use std::sync::Arc; + +#[derive(Debug)] +pub enum CplexMasterError { + Cplex(cpx::CplexError), + SubgradientExtension(Box), + NoMinorants, +} + +impl std::fmt::Display for CplexMasterError { + fn fmt(&self, fmt: &mut std::fmt::Formatter) -> std::result::Result<(), std::fmt::Error> { + use CplexMasterError::*; + match self { + Cplex(err) => err.fmt(fmt), + SubgradientExtension(err) => write!(fmt, "Subgradient extension failed: {}", err), + NoMinorants => write!(fmt, "Master problem contains no minorants"), + } + } +} + +impl std::error::Error for CplexMasterError { + fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { + use CplexMasterError::*; + match self { + Cplex(err) => Some(err), + SubgradientExtension(err) => Some(err.as_ref()), + NoMinorants => None, + } + } +} + +impl From for CplexMasterError { + fn from(err: cpx::CplexError) -> Self { + CplexMasterError::Cplex(err) + } +} + +pub type Result = std::result::Result; + +/// A minorant and its unique index. +struct MinorantInfo { + minorant: Minorant, + index: usize, +} + +impl Deref for MinorantInfo { + type Target = Minorant; + fn deref(&self) -> &Minorant { + &self.minorant + } +} + +impl DerefMut for MinorantInfo { + fn deref_mut(&mut self) -> &mut Minorant { + &mut self.minorant + } +} + +/// Maps a global index to a minorant. +#[derive(Clone, Copy)] +struct MinorantIdx { + /// The function (subproblem) index. + fidx: usize, + /// The minorant index within the subproblem. + idx: usize, +} + +/// A submodel. +/// +/// A submodel is a list of subproblems being aggregated in that model. +#[derive(Debug, Clone, Default)] +struct SubModel { + /// The list of subproblems. + subproblems: Vec, + + /// The number of minorants in this subproblem. + /// + /// The is the minimal number of minorants in each contained subproblem. + num_mins: usize, + /// The aggregated minorants of this submodel. + /// + /// This is just the sum of the corresponding minorants of the single + /// functions contained in this submodel. + minorants: Vec, +} + +impl Deref for SubModel { + type Target = Vec; + fn deref(&self) -> &Vec { + &self.subproblems + } +} + +impl DerefMut for SubModel { + fn deref_mut(&mut self) -> &mut Vec { + &mut self.subproblems + } +} + +pub struct CplexMaster { + env: *mut cpx::Env, + + lp: *mut cpx::Lp, + + /// True if the QP must be updated. + force_update: bool, + + /// List of free minorant indices. + freeinds: Vec, + + /// List of minorant indices to be updated. + updateinds: Vec, + + /// Mapping index to minorant. + index2min: Vec, + + /// The quadratic term. + qterm: Vec, + + /// The additional diagonal term to ensure positive definiteness. + qdiag: Real, + + /// The weight of the quadratic term. + weight: Real, + + /// The minorants for each subproblem in the model. + minorants: Vec>, + + /// The callback for the submodel for each subproblem. + select_model: Arc usize>, + + /// For each submodel the list of subproblems contained in that model. + submodels: Vec, + /// For each subproblem the submodel it is contained in. + in_submodel: Vec, + + /// Optimal multipliers for each subproblem in the model. + opt_mults: Vec, + /// Optimal aggregated minorant. + opt_minorant: Minorant, + + /// Maximal bundle size. + pub max_bundle_size: usize, +} + +unsafe impl Send for CplexMaster {} + +impl Drop for CplexMaster { + fn drop(&mut self) { + unsafe { cpx::freeprob(self.env, &mut self.lp) }; + unsafe { cpx::closeCPLEX(&mut self.env) }; + } +} + +impl UnconstrainedMasterProblem for CplexMaster { + type MinorantIndex = usize; + + type Err = CplexMasterError; + + fn new() -> Result { + let env; + + trycpx!({ + let mut status = 0; + env = cpx::openCPLEX(&mut status); + status + }); + + Ok(CplexMaster { + env, + lp: ptr::null_mut(), + force_update: true, + freeinds: vec![], + updateinds: vec![], + index2min: vec![], + qterm: vec![], + qdiag: 0.0, + weight: 1.0, + minorants: vec![], + select_model: Arc::new(|i| i), + submodels: vec![], + in_submodel: vec![], + opt_mults: vec![], + opt_minorant: Minorant::default(), + max_bundle_size: 50, + }) + } + + fn num_subproblems(&self) -> usize { + self.minorants.len() + } + + fn set_num_subproblems(&mut self, n: usize) -> Result<()> { + trycpx!(cpx::setintparam( + self.env, + cpx::Param::Qpmethod.to_c(), + cpx::Alg::Barrier.to_c() + )); + trycpx!(cpx::setdblparam(self.env, cpx::Param::Barepcomp.to_c(), 1e-12)); + + self.index2min.clear(); + self.freeinds.clear(); + self.minorants = (0..n).map(|_| vec![]).collect(); + self.opt_mults = vec![dvec![]; n]; + self.update_submodels(); + + Ok(()) + } + + fn weight(&self) -> Real { + self.weight + } + + fn set_weight(&mut self, weight: Real) -> Result<()> { + assert!(weight > 0.0); + self.weight = weight; + Ok(()) + } + + fn num_minorants(&self, fidx: usize) -> usize { + self.minorants[fidx].len() + } + + fn compress(&mut self, f: F) -> Result<()> + where + F: FnMut(Self::MinorantIndex, &mut dyn Iterator), + { + assert!(self.max_bundle_size >= 2, "Maximal bundle size must be >= 2"); + let mut f = f; + for i in 0..self.num_subproblems() { + let n = self.num_minorants(i); + if n >= self.max_bundle_size { + // aggregate minorants with smallest coefficients + let mut inds = (0..n).collect::>(); + inds.sort_by_key(|&j| -((1e6 * self.opt_mults[i][j]) as isize)); + let inds = inds[self.max_bundle_size - 2..] + .iter() + .map(|&j| self.minorants[i][j].index) + .collect::>(); + let (newindex, coeffs) = self.aggregate(i, &inds)?; + f(newindex, &mut inds.into_iter().zip(coeffs)); + } + } + Ok(()) + } + + fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result { + debug!("Add minorant"); + debug!(" fidx={} index={}: {}", fidx, self.minorants[fidx].len(), minorant); + + // find new unique minorant index + let min_idx = self.minorants[fidx].len(); + let index = if let Some(index) = self.freeinds.pop() { + self.index2min[index] = MinorantIdx { fidx, idx: min_idx }; + index + } else { + self.index2min.push(MinorantIdx { fidx, idx: min_idx }); + self.index2min.len() - 1 + }; + self.updateinds.push(index); + + // store minorant + self.minorants[fidx].push(MinorantInfo { minorant, index }); + self.opt_mults[fidx].push(0.0); + + Ok(index) + } + + fn add_vars( + &mut self, + nnew: usize, + changed: &[usize], + extend_subgradient: &mut SubgradientExtension, + ) -> Result<()> { + debug_assert!(!self.minorants[0].is_empty()); + if changed.is_empty() && nnew == 0 { + return Ok(()); + } + let noldvars = self.minorants[0][0].linear.len(); + let nnewvars = noldvars + nnew; + + let mut changedvars = vec![]; + changedvars.extend_from_slice(changed); + changedvars.extend(noldvars..nnewvars); + for (fidx, mins) in self.minorants.iter_mut().enumerate() { + for m in &mut mins[..] { + let new_subg = + extend_subgradient(fidx, m.index, &changedvars).map_err(CplexMasterError::SubgradientExtension)?; + for (&j, &g) in changed.iter().zip(new_subg.iter()) { + m.linear[j] = g; + } + m.linear.extend_from_slice(&new_subg[changed.len()..]); + } + } + + // update qterm because all minorants have changed + self.force_update = true; + + Ok(()) + } + + fn solve(&mut self, eta: &DVector, _fbound: Real, _augbound: Real, _relprec: Real) -> Result<()> { + if self.force_update || !self.updateinds.is_empty() { + self.init_qp()?; + } + + let nvars = unsafe { cpx::getnumcols(self.env, self.lp) as usize }; + debug_assert_eq!( + nvars, + self.submodels + .iter() + .map(|funs| funs.iter().map(|&fidx| self.minorants[fidx].len()).min().unwrap_or(0)) + .sum::() + ); + if nvars == 0 { + return Err(CplexMasterError::NoMinorants); + } + // update linear costs + { + let mut c = Vec::with_capacity(nvars); + let mut inds = Vec::with_capacity(nvars); + for submodel in &self.submodels { + for i in 0..submodel.num_mins { + let m = &submodel.minorants[i]; + let cost = -m.constant * self.weight - m.linear.dot(eta); + inds.push(c.len() as c_int); + c.push(cost); + } + } + debug_assert_eq!(inds.len(), nvars); + trycpx!(cpx::chgobj( + self.env, + self.lp, + nvars as c_int, + inds.as_ptr(), + c.as_ptr() + )); + } + + trycpx!(cpx::qpopt(self.env, self.lp)); + let mut sol = vec![0.0; nvars]; + trycpx!(cpx::getx(self.env, self.lp, sol.as_mut_ptr(), 0, nvars as c_int - 1)); + + let mut idx = 0; + let mut mults = Vec::with_capacity(nvars); + let mut mins = Vec::with_capacity(nvars); + + for submodel in &self.submodels { + for i in 0..submodel.num_mins { + for &fidx in submodel.iter() { + self.opt_mults[fidx][i] = sol[idx]; + mults.push(sol[idx]); + mins.push(&self.minorants[fidx][i].minorant); + } + idx += 1; + } + // set all multipliers for unused minorants to 0 + for &fidx in submodel.iter() { + for mult in &mut self.opt_mults[fidx][submodel.num_mins..] { + *mult = 0.0; + } + } + } + + self.opt_minorant = Aggregatable::combine(mults.into_iter().zip(mins)); + + Ok(()) + } + + fn dualopt(&self) -> &DVector { + &self.opt_minorant.linear + } + + fn dualopt_cutval(&self) -> Real { + self.opt_minorant.constant + } + + fn multiplier(&self, min: usize) -> Real { + let MinorantIdx { fidx, idx } = self.index2min[min]; + self.opt_mults[fidx][idx] + } + + fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box + 'a> { + Box::new( + self.opt_mults[fidx] + .iter() + .enumerate() + .map(move |(i, alpha)| (self.minorants[fidx][i].index, *alpha)), + ) + } + + fn eval_model(&self, y: &DVector) -> Real { + let mut result = 0.0; + for submodel in &self.submodels { + let mut this_val = NEG_INFINITY; + for m in &submodel.minorants[0..submodel.num_mins] { + this_val = this_val.max(m.eval(y)); + } + result += this_val; + } + result + } + + fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector)> { + assert!(!mins.is_empty(), "No minorants specified to be aggregated"); + + if mins.len() == 1 { + return Ok((mins[0], dvec![1.0])); + } + + // scale coefficients + let mut sum_coeffs = 0.0; + for &i in mins { + debug_assert_eq!( + fidx, self.index2min[i].fidx, + "Minorant {} does not belong to subproblem {} (belongs to: {})", + i, fidx, self.index2min[i].fidx + ); + sum_coeffs += self.opt_mults[fidx][self.index2min[i].idx]; + } + let aggr_coeffs = if sum_coeffs != 0.0 { + mins.iter() + .map(|&i| self.opt_mults[fidx][self.index2min[i].idx] / sum_coeffs) + .collect::() + } else { + dvec![0.0; mins.len()] + }; + + // Compute the aggregated minorant. + let aggr = Aggregatable::combine( + aggr_coeffs.iter().cloned().zip( + mins.iter() + .map(|&index| &self.minorants[fidx][self.index2min[index].idx].minorant), + ), + ); + + // Remove the minorants that have been aggregated. + for &i in mins { + let MinorantIdx { + fidx: min_fidx, + idx: min_idx, + } = self.index2min[i]; + debug_assert_eq!( + fidx, min_fidx, + "Minorant {} does not belong to subproblem {} (belongs to: {})", + i, fidx, min_fidx + ); + + let m = self.minorants[fidx].swap_remove(min_idx); + self.opt_mults[fidx].swap_remove(min_idx); + self.freeinds.push(m.index); + debug_assert_eq!(m.index, i); + + // Update index2min table and mark qterm to be updated. + // This is only necessary if the removed minorant was not the last one. + if min_idx < self.minorants[fidx].len() { + self.index2min[self.minorants[fidx][min_idx].index].idx = min_idx; + self.updateinds.push(self.minorants[fidx][min_idx].index); + } + } + + // Finally add the aggregated minorant. + let aggr_idx = self.add_minorant(fidx, aggr)?; + Ok((aggr_idx, aggr_coeffs)) + } + + fn move_center(&mut self, alpha: Real, d: &DVector) { + for mins in &mut self.minorants { + for m in mins.iter_mut() { + m.move_center(alpha, d); + } + } + for submod in &mut self.submodels { + for m in &mut submod.minorants { + m.move_center(alpha, d); + } + } + } +} + +impl CplexMaster { + /// Set a custom submodel selector. + /// + /// For each subproblem index the selector should return a submodel index. + /// All subproblems with the same submodel index are aggregated in a single + /// cutting plane model. + fn set_submodel_selection(&mut self, selector: F) + where + F: Fn(usize) -> usize + 'static, + { + self.select_model = Arc::new(selector); + self.update_submodels(); + //unimplemented!("Arbitrary submodels are not implemented, yet"); + } + + fn update_submodels(&mut self) { + self.submodels.clear(); + self.in_submodel.resize(self.num_subproblems(), 0); + for fidx in 0..self.num_subproblems() { + let model_idx = (self.select_model)(fidx); + if model_idx >= self.submodels.len() { + self.submodels.resize_with(model_idx + 1, SubModel::default); + } + self.submodels[model_idx].push(fidx); + self.in_submodel[fidx] = model_idx; + } + } + + /// Use a fully disaggregated model. + /// + /// A fully disaggregated model has one separate submodel for each subproblem. + /// Hence, calling this method is equivalent to + /// `CplexMaster::set_submodel_selection(|i| i)`. + pub fn use_full_disaggregation(&mut self) { + self.set_submodel_selection(|i| i) + } + + /// Use a fully aggregated model. + /// + /// A fully aggregated model has one submodel for all subproblems. + /// Hence, calling this method is equivalent to + /// `CplexMaster::set_submodel_selection(|_| 0)`. + pub fn use_full_aggregation(&mut self) { + self.set_submodel_selection(|_| 0) + } + + fn init_qp(&mut self) -> Result<()> { + if self.force_update { + self.updateinds.clear(); + for mins in &self.minorants { + self.updateinds.extend(mins.iter().map(|m| m.index)); + } + } + + let minorants = &self.minorants; + + // Compute the number of minorants in each submodel. + for submodel in self.submodels.iter_mut() { + submodel.num_mins = submodel.iter().map(|&fidx| minorants[fidx].len()).min().unwrap_or(0); + submodel.minorants.resize_with(submodel.num_mins, Minorant::default); + } + + // Only minorants belonging to the first subproblem of each submodel + // must be updated. + // + // We filter all indices that + // 1. belong to a subproblem being the first in its model + // 2. is a valid index (< minimal number of minorants within this submodel) + // and map them to (index, model_index, minorant_index) where + // - `index` is the variable index (i.e. minorant index of first subproblem) + // - `model_index` is the submodel index + // - `minorant_index` is the index of the minorant within the model + let updateinds = self + .updateinds + .iter() + .filter_map(|&index| { + let MinorantIdx { fidx, idx } = self.index2min[index]; + let mod_i = self.in_submodel[fidx]; + let submodel = &self.submodels[mod_i]; + if submodel[0] == fidx && idx < submodel.num_mins { + Some((index, mod_i, idx)) + } else { + None + } + }) + .collect::>(); + + // Compute the aggregated minorants. + for &(_, mod_i, i) in &updateinds { + let submodel = &mut self.submodels[mod_i]; + submodel.minorants[i] = + Aggregatable::combine(submodel.iter().map(|&fidx| (1.0, &minorants[fidx][i].minorant))); + } + + let submodels = &self.submodels; + + // Build quadratic term, this is for all pairs of minorants where each g_i + // is an aggregated minorant of a submodel. + // + // For simplicity we always store the terms in the index of the minorant of the first + // subproblem in each submodel. + let ntotalminorants = self.index2min.len(); + if ntotalminorants > self.qterm.len() { + self.qterm.resize(ntotalminorants, dvec![]); + for i in 0..self.qterm.len() { + self.qterm[i].resize(ntotalminorants, 0.0); + } + } + + // - i is the number of the minorant within the submodel submodel_i + // - idx_i is the unique index of that minorant (of the first subproblem) + // - j is the number of the minorant within the submodel submodel_j + // - idx_j is the unique index of that minorant (of the first subproblem) + for submodel_i in submodels.iter() { + // Compute the minorant g_i for each i + for i in 0..submodel_i.num_mins { + // Store the computed values at the index of the first subproblem in this model. + let idx_i = minorants[submodel_i[0]][i].index; + let g_i = &submodel_i.minorants[i].linear; + // Now compute the inner product with each other minorant + // that has to be updated. + for &(idx_j, mod_j, j) in updateinds.iter() { + let x = submodels[mod_j].minorants[j].linear.dot(g_i); + self.qterm[idx_i][idx_j] = x; + self.qterm[idx_j][idx_i] = x; + } + } + } + + // We verify that the qterm is correct + if cfg!(debug_assertions) { + for submod_i in submodels.iter() { + for i in 0..submod_i.num_mins { + let idx_i = self.minorants[submod_i[0]][i].index; + for submod_j in submodels.iter() { + for j in 0..submod_j.num_mins { + let idx_j = self.minorants[submod_j[0]][j].index; + let x = submod_i.minorants[i].linear.dot(&submod_j.minorants[j].linear); + debug_assert!((x - self.qterm[idx_i][idx_j]) < 1e-6); + } + } + } + } + } + + // main diagonal plus small identity to ensure Q being semi-definite + self.qdiag = 0.0; + for submodel in submodels.iter() { + for i in 0..submodel.num_mins { + let idx = minorants[submodel[0]][i].index; + self.qdiag = Real::max(self.qdiag, self.qterm[idx][idx]); + } + } + self.qdiag *= 1e-8; + + // We have updated everything. + self.updateinds.clear(); + self.force_update = false; + + self.init_cpx_qp() + } + + fn init_cpx_qp(&mut self) -> Result<()> { + if !self.lp.is_null() { + trycpx!(cpx::freeprob(self.env, &mut self.lp)); + } + trycpx!({ + let mut status = 0; + self.lp = cpx::createprob(self.env, &mut status, c_str!("mastercp").as_ptr()); + status + }); + + let nsubmodels = self.submodels.len(); + let submodels = &self.submodels; + let minorants = &self.minorants; + + // add convexity constraints + let sense: Vec = vec!['E' as c_char; nsubmodels]; + let rhs = dvec![1.0; nsubmodels]; + let mut rmatbeg = Vec::with_capacity(nsubmodels); + let mut rmatind = Vec::with_capacity(self.index2min.len()); + let mut rmatval = Vec::with_capacity(self.index2min.len()); + + let mut nvars = 0; + for submodel in submodels.iter() { + if submodel.is_empty() { + // this should only happen if the submodel selector leaves + // holes + continue; + } + + rmatbeg.push(nvars as c_int); + rmatind.extend((nvars as c_int..).take(submodel.num_mins)); + rmatval.extend(repeat(1.0).take(submodel.num_mins)); + nvars += submodel.num_mins; + } + + trycpx!(cpx::addrows( + self.env, + self.lp, + nvars as c_int, + rmatbeg.len() as c_int, + rmatind.len() as c_int, + rhs.as_ptr(), + sense.as_ptr(), + rmatbeg.as_ptr(), + rmatind.as_ptr(), + rmatval.as_ptr(), + ptr::null(), + ptr::null() + )); + + // update coefficients + let mut var_i = 0; + for (mod_i, submodel_i) in submodels.iter().enumerate() { + for i in 0..submodel_i.num_mins { + let idx_i = minorants[submodel_i[0]][i].index; + let mut var_j = 0; + for (mod_j, submodel_j) in submodels.iter().enumerate() { + for j in 0..submodel_j.num_mins { + let idx_j = minorants[submodel_j[0]][j].index; + let q = self.qterm[idx_i][idx_j] + if mod_i != mod_j || i != j { 0.0 } else { self.qdiag }; + trycpx!(cpx::chgqpcoef(self.env, self.lp, var_i as c_int, var_j as c_int, q)); + var_j += 1; + } + } + var_i += 1; + } + } + + Ok(()) + } +} + +pub struct Builder { + /// The maximal bundle size used in the master problem. + pub max_bundle_size: usize, + /// The submodel selector. + select_model: Arc usize>, +} + +impl Default for Builder { + fn default() -> Self { + Builder { + max_bundle_size: 50, + select_model: Arc::new(|i| i), + } + } +} + +impl super::Builder for Builder { + type MasterProblem = CplexMaster; + + fn build(&mut self) -> Result { + let mut cpx = CplexMaster::new()?; + cpx.max_bundle_size = self.max_bundle_size; + cpx.select_model = self.select_model.clone(); + cpx.update_submodels(); + Ok(cpx) + } +} + +impl Builder { + pub fn max_bundle_size(&mut self, s: usize) -> &mut Self { + assert!(s >= 2, "The maximal bundle size must be >= 2"); + self.max_bundle_size = s; + self + } + + /// Set a custom submodel selector. + /// + /// For each subproblem index the selector should return a submodel index. + /// All subproblems with the same submodel index are aggregated in a single + /// cutting plane model. + pub fn submodel_selection(&mut self, selector: F) -> &mut Self + where + F: Fn(usize) -> usize + 'static, + { + self.select_model = Arc::new(selector); + self + } + + /// Use a fully disaggregated model. + /// + /// A fully disaggregated model has one separate submodel for each subproblem. + /// Hence, calling this method is equivalent to + /// `CplexMaster::set_submodel_selection(|i| i)`. + pub fn use_full_disaggregation(&mut self) -> &mut Self { + self.submodel_selection(|i| i) + } + + /// Use a fully aggregated model. + /// + /// A fully aggregated model has one submodel for all subproblems. + /// Hence, calling this method is equivalent to + /// `CplexMaster::set_submodel_selection(|_| 0)`. + pub fn use_full_aggregation(&mut self) -> &mut Self { + self.submodel_selection(|_| 0) + } +} ADDED src/master/boxed/unconstrained/minimal.rs Index: src/master/boxed/unconstrained/minimal.rs ================================================================== --- /dev/null +++ src/master/boxed/unconstrained/minimal.rs @@ -0,0 +1,401 @@ +// Copyright (c) 2016, 2017, 2019 Frank Fischer +// +// This program is free software: you can redistribute it and/or +// modify it under the terms of the GNU General Public License as +// published by the Free Software Foundation, either version 3 of the +// License, or (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, but +// WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program. If not, see +// + +use super::{SubgradientExtension, UnconstrainedMasterProblem}; +use crate::{Aggregatable, DVector, Minorant, Real}; + +use log::debug; + +use std::error::Error; +use std::f64::NEG_INFINITY; +use std::fmt; +use std::result; + +/// Minimal master problem error. +#[derive(Debug)] +pub enum MinimalMasterError { + NoMinorants, + MaxMinorants { subproblem: usize }, + SubgradientExtension(Box), +} + +impl fmt::Display for MinimalMasterError { + fn fmt(&self, fmt: &mut fmt::Formatter) -> result::Result<(), fmt::Error> { + use self::MinimalMasterError::*; + match self { + MaxMinorants { subproblem } => write!( + fmt, + "The minimal master problem allows at most two minorants (subproblem: {})", + subproblem + ), + NoMinorants => write!(fmt, "The master problem does not contain a minorant"), + SubgradientExtension(err) => write!(fmt, "Subgradient extension failed: {}", err), + } + } +} + +impl Error for MinimalMasterError { + fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { + use MinimalMasterError::*; + match self { + SubgradientExtension(err) => Some(err.as_ref()), + _ => None, + } + } +} + +/** + * A minimal master problem with only two minorants. + * + * This is the simplest possible master problem for bundle methods. It + * has only two minorants and only one function model. The advantage + * is that this model can be solved explicitely and very quickly, but + * it is only a very loose approximation of the objective function. + * + * Because of its properties, it can only be used if the problem to be + * solved has a maximal number of minorants of two and only one + * subproblem. + */ +pub struct MinimalMaster { + /// The weight of the quadratic term. + weight: Real, + + /// The minorants in the model. + /// + /// There are up to two minorants, each minorant consists of one part for + /// each subproblem. + /// + /// The minorants for the i-th subproblem have the indices `2*i` and + /// `2*i+1`. + minorants: [Vec; 2], + /// The number of minorants. Only the minorants with index less than this + /// number are valid. + num_minorants: usize, + /// The number of minorants for each subproblem. + num_minorants_of: Vec, + /// The number of subproblems. + num_subproblems: usize, + /// The number of subproblems with at least 1 minorant. + num_subproblems_with_1: usize, + /// The number of subproblems with at least 2 minorants. + num_subproblems_with_2: usize, + /// Optimal multipliers. + opt_mult: [Real; 2], + /// Optimal aggregated minorant. + opt_minorant: Minorant, +} + +impl UnconstrainedMasterProblem for MinimalMaster { + type MinorantIndex = usize; + + type Err = MinimalMasterError; + + fn new() -> Result { + Ok(MinimalMaster { + weight: 1.0, + num_minorants: 0, + num_minorants_of: vec![], + num_subproblems: 0, + num_subproblems_with_1: 0, + num_subproblems_with_2: 0, + minorants: [vec![], vec![]], + opt_mult: [0.0, 0.0], + opt_minorant: Minorant::default(), + }) + } + + fn num_subproblems(&self) -> usize { + self.num_subproblems + } + + fn set_num_subproblems(&mut self, n: usize) -> Result<(), Self::Err> { + self.num_subproblems = n; + self.num_minorants = 0; + self.num_minorants_of = vec![0; n]; + self.num_subproblems_with_1 = 0; + self.num_subproblems_with_2 = 0; + self.minorants = [vec![Minorant::default(); n], vec![Minorant::default(); n]]; + Ok(()) + } + + fn compress(&mut self, f: F) -> Result<(), Self::Err> + where + F: FnMut(Self::MinorantIndex, &mut dyn Iterator), + { + if self.num_minorants == 2 { + debug!("Aggregate"); + debug!(" {} * {:?}", self.opt_mult[0], self.minorants[0]); + debug!(" {} * {:?}", self.opt_mult[1], self.minorants[1]); + + let mut f = f; + for fidx in 0..self.num_subproblems { + f( + 2 * fidx, + &mut self + .opt_mult + .iter() + .enumerate() + .map(|(i, alpha)| (2 * fidx + i, *alpha)), + ); + } + + self.minorants[0] = Aggregatable::combine(self.opt_mult.iter().cloned().zip(&self.minorants)); + self.opt_mult[0] = 1.0; + self.num_minorants = 1; + self.num_minorants_of.clear(); + self.num_minorants_of.resize(self.num_subproblems, 1); + self.num_subproblems_with_2 = 0; + + debug!(" {:?}", self.minorants[0]); + } + Ok(()) + } + + fn weight(&self) -> Real { + self.weight + } + + fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err> { + assert!(weight > 0.0); + self.weight = weight; + Ok(()) + } + + fn num_minorants(&self, fidx: usize) -> usize { + self.num_minorants_of[fidx] + } + + fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result { + if self.num_minorants_of[fidx] >= 2 { + return Err(MinimalMasterError::MaxMinorants { subproblem: fidx }); + } + + let minidx = self.num_minorants_of[fidx]; + self.num_minorants_of[fidx] += 1; + self.minorants[minidx][fidx] = minorant; + + match minidx { + 0 => { + self.num_subproblems_with_1 += 1; + if self.num_subproblems_with_1 == self.num_subproblems { + self.num_minorants = 1; + self.opt_mult[0] = 0.0; + } + Ok(2 * fidx) + } + 1 => { + self.num_subproblems_with_2 += 1; + if self.num_subproblems_with_2 == self.num_subproblems { + self.num_minorants = 2; + self.opt_mult[1] = 0.0; + } + Ok(2 * fidx + 1) + } + _ => unreachable!("Invalid number of minorants in subproblem {}", fidx), + } + } + + fn add_vars( + &mut self, + nnew: usize, + changed: &[usize], + extend_subgradient: &mut SubgradientExtension, + ) -> Result<(), Self::Err> { + if self.num_subproblems_with_1 == 0 { + return Ok(()); + } + + let noldvars = self.minorants[0][self.num_minorants_of.iter().position(|&n| n > 0).unwrap()] + .linear + .len(); + let mut changedvars = vec![]; + changedvars.extend_from_slice(changed); + changedvars.extend(noldvars..noldvars + nnew); + + for fidx in 0..self.num_subproblems { + for i in 0..self.num_minorants_of[fidx] { + let new_subg = extend_subgradient(fidx, 2 * fidx + i, &changedvars) + .map_err(MinimalMasterError::SubgradientExtension)?; + let m = &mut self.minorants[i][fidx]; + for (&j, &g) in changed.iter().zip(new_subg.iter()) { + m.linear[j] = g; + } + m.linear.extend_from_slice(&new_subg[changed.len()..]); + } + } + + Ok(()) + } + + #[allow(unused_variables)] + fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<(), Self::Err> { + for fidx in 0..self.num_subproblems { + for i in 0..self.num_minorants_of[fidx] { + debug!(" min(fidx:{}, i:{}) = {}", fidx, i, self.minorants[i][fidx]); + } + } + + if self.num_minorants == 2 { + let min0 = Minorant::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx]))); + let min1 = Minorant::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[1][fidx]))); + let xx = min0.linear.dot(&min0.linear); + let yy = min1.linear.dot(&min1.linear); + let xy = min0.linear.dot(&min1.linear); + let xeta = min0.linear.dot(eta); + let yeta = min1.linear.dot(eta); + let a = yy - 2.0 * xy + xx; + let b = xy - xx - yeta + xeta; + + let mut alpha2 = 0.0; + if a > 0.0 { + alpha2 = ((min1.constant - min0.constant) * self.weight - b) / a; + alpha2 = alpha2.max(0.0).min(1.0); + } + self.opt_mult[0] = 1.0 - alpha2; + self.opt_mult[1] = alpha2; + self.opt_minorant = Aggregatable::combine(self.opt_mult.iter().cloned().zip([min0, min1].iter())); + } else if self.num_minorants == 1 { + let min0 = Aggregatable::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx]))); + self.opt_minorant = min0; + self.opt_mult[0] = 1.0; + } else { + return Err(MinimalMasterError::NoMinorants); + } + + debug!("Unrestricted"); + debug!(" opt_minorant={}", self.opt_minorant); + debug!(" opt_mult={:?}", &self.opt_mult[0..self.num_minorants]); + + Ok(()) + } + + fn dualopt(&self) -> &DVector { + &self.opt_minorant.linear + } + + fn dualopt_cutval(&self) -> Real { + self.opt_minorant.constant + } + + fn multiplier(&self, min: usize) -> Real { + self.opt_mult[min % 2] + } + + fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box + 'a> { + Box::new( + self.opt_mult + .iter() + .take(self.num_minorants_of[fidx]) + .enumerate() + .map(move |(i, alpha)| (2 * fidx + i, *alpha)), + ) + } + + fn eval_model(&self, y: &DVector) -> Real { + let mut result = NEG_INFINITY; + for mins in &self.minorants[0..self.num_minorants] { + result = result.max(mins.iter().map(|m| m.eval(y)).sum()); + } + result + } + + fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector), Self::Err> { + debug!("Aggregate minorants {:?} of subproblem {}", mins, fidx); + if mins.len() == 2 { + debug_assert_ne!(mins[0], mins[1], "Minorants to be aggregated must be different"); + debug_assert_eq!( + mins[0] / 2, + fidx, + "Minorant {} does not belong to subproblem {}", + mins[0], + fidx + ); + debug_assert_eq!( + mins[1] / 2, + fidx, + "Minorant {} does not belong to subproblem {}", + mins[1], + fidx + ); + debug_assert!( + mins[0] % 2 < self.num_minorants_of[fidx], + "Invalid minorant index for subproblem {}: {}", + fidx, + mins[0] + ); + debug_assert!( + mins[1] % 2 < self.num_minorants_of[fidx], + "Invalid minorant index for subproblem {}: {}", + fidx, + mins[1] + ); + + let min0 = mins[0] % 2; + let min1 = mins[1] % 2; + + debug!("Aggregate"); + debug!(" {} * {}", self.opt_mult[min0], self.minorants[min0][fidx]); + debug!(" {} * {}", self.opt_mult[min1], self.minorants[min1][fidx]); + self.minorants[0][fidx] = Aggregatable::combine( + [ + (self.opt_mult[min0], &self.minorants[min0][fidx]), + (self.opt_mult[min1], &self.minorants[min1][fidx]), + ] + .iter() + .cloned(), + ); + + if self.num_subproblems_with_2 == self.num_subproblems { + self.num_minorants -= 1; + } + self.num_subproblems_with_2 -= 1; + self.num_minorants_of[fidx] -= 1; + + let coeffs = dvec![self.opt_mult[min0], self.opt_mult[min1]]; + + debug!(" {}", self.minorants[0][fidx]); + Ok((2 * fidx, coeffs)) + } else if mins.len() == 1 { + Ok((mins[0], dvec![1.0])) + } else { + panic!("No minorants specified to be aggregated"); + } + } + + fn move_center(&mut self, alpha: Real, d: &DVector) { + for fidx in 0..self.num_subproblems { + for i in 0..self.num_minorants_of[fidx] { + self.minorants[i][fidx].move_center(alpha, d); + } + } + } +} + +pub struct Builder; + +impl Default for Builder { + fn default() -> Self { + Builder + } +} + +impl super::Builder for Builder { + type MasterProblem = MinimalMaster; + + fn build(&mut self) -> Result { + MinimalMaster::new() + } +} DELETED src/master/cpx.rs Index: src/master/cpx.rs ================================================================== --- src/master/cpx.rs +++ /dev/null @@ -1,815 +0,0 @@ -// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -//! Master problem implementation using CPLEX. - -#![allow(unused_unsafe)] - -use crate::master::unconstrained::{self, UnconstrainedMasterProblem}; -use crate::master::SubgradientExtension; -use crate::{Aggregatable, DVector, Minorant, Real}; - -use c_str_macro::c_str; -use cplex_sys as cpx; -use cplex_sys::trycpx; -use log::debug; - -use std; -use std::f64::NEG_INFINITY; -use std::iter::repeat; -use std::ops::{Deref, DerefMut}; -use std::os::raw::{c_char, c_int}; -use std::ptr; -use std::sync::Arc; - -#[derive(Debug)] -pub enum CplexMasterError { - Cplex(cpx::CplexError), - SubgradientExtension(Box), - NoMinorants, -} - -impl std::fmt::Display for CplexMasterError { - fn fmt(&self, fmt: &mut std::fmt::Formatter) -> std::result::Result<(), std::fmt::Error> { - use CplexMasterError::*; - match self { - Cplex(err) => err.fmt(fmt), - SubgradientExtension(err) => write!(fmt, "Subgradient extension failed: {}", err), - NoMinorants => write!(fmt, "Master problem contains no minorants"), - } - } -} - -impl std::error::Error for CplexMasterError { - fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { - use CplexMasterError::*; - match self { - Cplex(err) => Some(err), - SubgradientExtension(err) => Some(err.as_ref()), - NoMinorants => None, - } - } -} - -impl From for CplexMasterError { - fn from(err: cpx::CplexError) -> Self { - CplexMasterError::Cplex(err) - } -} - -pub type Result = std::result::Result; - -/// A minorant and its unique index. -struct MinorantInfo { - minorant: Minorant, - index: usize, -} - -impl Deref for MinorantInfo { - type Target = Minorant; - fn deref(&self) -> &Minorant { - &self.minorant - } -} - -impl DerefMut for MinorantInfo { - fn deref_mut(&mut self) -> &mut Minorant { - &mut self.minorant - } -} - -/// Maps a global index to a minorant. -#[derive(Clone, Copy)] -struct MinorantIdx { - /// The function (subproblem) index. - fidx: usize, - /// The minorant index within the subproblem. - idx: usize, -} - -/// A submodel. -/// -/// A submodel is a list of subproblems being aggregated in that model. -#[derive(Debug, Clone, Default)] -struct SubModel { - /// The list of subproblems. - subproblems: Vec, - - /// The number of minorants in this subproblem. - /// - /// The is the minimal number of minorants in each contained subproblem. - num_mins: usize, - /// The aggregated minorants of this submodel. - /// - /// This is just the sum of the corresponding minorants of the single - /// functions contained in this submodel. - minorants: Vec, -} - -impl Deref for SubModel { - type Target = Vec; - fn deref(&self) -> &Vec { - &self.subproblems - } -} - -impl DerefMut for SubModel { - fn deref_mut(&mut self) -> &mut Vec { - &mut self.subproblems - } -} - -pub struct CplexMaster { - env: *mut cpx::Env, - - lp: *mut cpx::Lp, - - /// True if the QP must be updated. - force_update: bool, - - /// List of free minorant indices. - freeinds: Vec, - - /// List of minorant indices to be updated. - updateinds: Vec, - - /// Mapping index to minorant. - index2min: Vec, - - /// The quadratic term. - qterm: Vec, - - /// The additional diagonal term to ensure positive definiteness. - qdiag: Real, - - /// The weight of the quadratic term. - weight: Real, - - /// The minorants for each subproblem in the model. - minorants: Vec>, - - /// The callback for the submodel for each subproblem. - select_model: Arc usize>, - - /// For each submodel the list of subproblems contained in that model. - submodels: Vec, - /// For each subproblem the submodel it is contained in. - in_submodel: Vec, - - /// Optimal multipliers for each subproblem in the model. - opt_mults: Vec, - /// Optimal aggregated minorant. - opt_minorant: Minorant, - - /// Maximal bundle size. - pub max_bundle_size: usize, -} - -unsafe impl Send for CplexMaster {} - -impl Drop for CplexMaster { - fn drop(&mut self) { - unsafe { cpx::freeprob(self.env, &mut self.lp) }; - unsafe { cpx::closeCPLEX(&mut self.env) }; - } -} - -impl UnconstrainedMasterProblem for CplexMaster { - type MinorantIndex = usize; - - type Err = CplexMasterError; - - fn new() -> Result { - let env; - - trycpx!({ - let mut status = 0; - env = cpx::openCPLEX(&mut status); - status - }); - - Ok(CplexMaster { - env, - lp: ptr::null_mut(), - force_update: true, - freeinds: vec![], - updateinds: vec![], - index2min: vec![], - qterm: vec![], - qdiag: 0.0, - weight: 1.0, - minorants: vec![], - select_model: Arc::new(|i| i), - submodels: vec![], - in_submodel: vec![], - opt_mults: vec![], - opt_minorant: Minorant::default(), - max_bundle_size: 50, - }) - } - - fn num_subproblems(&self) -> usize { - self.minorants.len() - } - - fn set_num_subproblems(&mut self, n: usize) -> Result<()> { - trycpx!(cpx::setintparam( - self.env, - cpx::Param::Qpmethod.to_c(), - cpx::Alg::Barrier.to_c() - )); - trycpx!(cpx::setdblparam(self.env, cpx::Param::Barepcomp.to_c(), 1e-12)); - - self.index2min.clear(); - self.freeinds.clear(); - self.minorants = (0..n).map(|_| vec![]).collect(); - self.opt_mults = vec![dvec![]; n]; - self.update_submodels(); - - Ok(()) - } - - fn weight(&self) -> Real { - self.weight - } - - fn set_weight(&mut self, weight: Real) -> Result<()> { - assert!(weight > 0.0); - self.weight = weight; - Ok(()) - } - - fn num_minorants(&self, fidx: usize) -> usize { - self.minorants[fidx].len() - } - - fn compress(&mut self, f: F) -> Result<()> - where - F: FnMut(Self::MinorantIndex, &mut dyn Iterator), - { - assert!(self.max_bundle_size >= 2, "Maximal bundle size must be >= 2"); - let mut f = f; - for i in 0..self.num_subproblems() { - let n = self.num_minorants(i); - if n >= self.max_bundle_size { - // aggregate minorants with smallest coefficients - let mut inds = (0..n).collect::>(); - inds.sort_by_key(|&j| -((1e6 * self.opt_mults[i][j]) as isize)); - let inds = inds[self.max_bundle_size - 2..] - .iter() - .map(|&j| self.minorants[i][j].index) - .collect::>(); - let (newindex, coeffs) = self.aggregate(i, &inds)?; - f(newindex, &mut inds.into_iter().zip(coeffs)); - } - } - Ok(()) - } - - fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result { - debug!("Add minorant"); - debug!(" fidx={} index={}: {}", fidx, self.minorants[fidx].len(), minorant); - - // find new unique minorant index - let min_idx = self.minorants[fidx].len(); - let index = if let Some(index) = self.freeinds.pop() { - self.index2min[index] = MinorantIdx { fidx, idx: min_idx }; - index - } else { - self.index2min.push(MinorantIdx { fidx, idx: min_idx }); - self.index2min.len() - 1 - }; - self.updateinds.push(index); - - // store minorant - self.minorants[fidx].push(MinorantInfo { minorant, index }); - self.opt_mults[fidx].push(0.0); - - Ok(index) - } - - fn add_vars( - &mut self, - nnew: usize, - changed: &[usize], - extend_subgradient: &mut SubgradientExtension, - ) -> Result<()> { - debug_assert!(!self.minorants[0].is_empty()); - if changed.is_empty() && nnew == 0 { - return Ok(()); - } - let noldvars = self.minorants[0][0].linear.len(); - let nnewvars = noldvars + nnew; - - let mut changedvars = vec![]; - changedvars.extend_from_slice(changed); - changedvars.extend(noldvars..nnewvars); - for (fidx, mins) in self.minorants.iter_mut().enumerate() { - for m in &mut mins[..] { - let new_subg = - extend_subgradient(fidx, m.index, &changedvars).map_err(CplexMasterError::SubgradientExtension)?; - for (&j, &g) in changed.iter().zip(new_subg.iter()) { - m.linear[j] = g; - } - m.linear.extend_from_slice(&new_subg[changed.len()..]); - } - } - - // update qterm because all minorants have changed - self.force_update = true; - - Ok(()) - } - - fn solve(&mut self, eta: &DVector, _fbound: Real, _augbound: Real, _relprec: Real) -> Result<()> { - if self.force_update || !self.updateinds.is_empty() { - self.init_qp()?; - } - - let nvars = unsafe { cpx::getnumcols(self.env, self.lp) as usize }; - debug_assert_eq!( - nvars, - self.submodels - .iter() - .map(|funs| funs.iter().map(|&fidx| self.minorants[fidx].len()).min().unwrap_or(0)) - .sum::() - ); - if nvars == 0 { - return Err(CplexMasterError::NoMinorants); - } - // update linear costs - { - let mut c = Vec::with_capacity(nvars); - let mut inds = Vec::with_capacity(nvars); - for submodel in &self.submodels { - for i in 0..submodel.num_mins { - let m = &submodel.minorants[i]; - let cost = -m.constant * self.weight - m.linear.dot(eta); - inds.push(c.len() as c_int); - c.push(cost); - } - } - debug_assert_eq!(inds.len(), nvars); - trycpx!(cpx::chgobj( - self.env, - self.lp, - nvars as c_int, - inds.as_ptr(), - c.as_ptr() - )); - } - - trycpx!(cpx::qpopt(self.env, self.lp)); - let mut sol = vec![0.0; nvars]; - trycpx!(cpx::getx(self.env, self.lp, sol.as_mut_ptr(), 0, nvars as c_int - 1)); - - let mut idx = 0; - let mut mults = Vec::with_capacity(nvars); - let mut mins = Vec::with_capacity(nvars); - - for submodel in &self.submodels { - for i in 0..submodel.num_mins { - for &fidx in submodel.iter() { - self.opt_mults[fidx][i] = sol[idx]; - mults.push(sol[idx]); - mins.push(&self.minorants[fidx][i].minorant); - } - idx += 1; - } - // set all multipliers for unused minorants to 0 - for &fidx in submodel.iter() { - for mult in &mut self.opt_mults[fidx][submodel.num_mins..] { - *mult = 0.0; - } - } - } - - self.opt_minorant = Aggregatable::combine(mults.into_iter().zip(mins)); - - Ok(()) - } - - fn dualopt(&self) -> &DVector { - &self.opt_minorant.linear - } - - fn dualopt_cutval(&self) -> Real { - self.opt_minorant.constant - } - - fn multiplier(&self, min: usize) -> Real { - let MinorantIdx { fidx, idx } = self.index2min[min]; - self.opt_mults[fidx][idx] - } - - fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box + 'a> { - Box::new( - self.opt_mults[fidx] - .iter() - .enumerate() - .map(move |(i, alpha)| (self.minorants[fidx][i].index, *alpha)), - ) - } - - fn eval_model(&self, y: &DVector) -> Real { - let mut result = 0.0; - for submodel in &self.submodels { - let mut this_val = NEG_INFINITY; - for m in &submodel.minorants[0..submodel.num_mins] { - this_val = this_val.max(m.eval(y)); - } - result += this_val; - } - result - } - - fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector)> { - assert!(!mins.is_empty(), "No minorants specified to be aggregated"); - - if mins.len() == 1 { - return Ok((mins[0], dvec![1.0])); - } - - // scale coefficients - let mut sum_coeffs = 0.0; - for &i in mins { - debug_assert_eq!( - fidx, self.index2min[i].fidx, - "Minorant {} does not belong to subproblem {} (belongs to: {})", - i, fidx, self.index2min[i].fidx - ); - sum_coeffs += self.opt_mults[fidx][self.index2min[i].idx]; - } - let aggr_coeffs = if sum_coeffs != 0.0 { - mins.iter() - .map(|&i| self.opt_mults[fidx][self.index2min[i].idx] / sum_coeffs) - .collect::() - } else { - dvec![0.0; mins.len()] - }; - - // Compute the aggregated minorant. - let aggr = Aggregatable::combine( - aggr_coeffs.iter().cloned().zip( - mins.iter() - .map(|&index| &self.minorants[fidx][self.index2min[index].idx].minorant), - ), - ); - - // Remove the minorants that have been aggregated. - for &i in mins { - let MinorantIdx { - fidx: min_fidx, - idx: min_idx, - } = self.index2min[i]; - debug_assert_eq!( - fidx, min_fidx, - "Minorant {} does not belong to subproblem {} (belongs to: {})", - i, fidx, min_fidx - ); - - let m = self.minorants[fidx].swap_remove(min_idx); - self.opt_mults[fidx].swap_remove(min_idx); - self.freeinds.push(m.index); - debug_assert_eq!(m.index, i); - - // Update index2min table and mark qterm to be updated. - // This is only necessary if the removed minorant was not the last one. - if min_idx < self.minorants[fidx].len() { - self.index2min[self.minorants[fidx][min_idx].index].idx = min_idx; - self.updateinds.push(self.minorants[fidx][min_idx].index); - } - } - - // Finally add the aggregated minorant. - let aggr_idx = self.add_minorant(fidx, aggr)?; - Ok((aggr_idx, aggr_coeffs)) - } - - fn move_center(&mut self, alpha: Real, d: &DVector) { - for mins in &mut self.minorants { - for m in mins.iter_mut() { - m.move_center(alpha, d); - } - } - for submod in &mut self.submodels { - for m in &mut submod.minorants { - m.move_center(alpha, d); - } - } - } -} - -impl CplexMaster { - /// Set a custom submodel selector. - /// - /// For each subproblem index the selector should return a submodel index. - /// All subproblems with the same submodel index are aggregated in a single - /// cutting plane model. - fn set_submodel_selection(&mut self, selector: F) - where - F: Fn(usize) -> usize + 'static, - { - self.select_model = Arc::new(selector); - self.update_submodels(); - //unimplemented!("Arbitrary submodels are not implemented, yet"); - } - - fn update_submodels(&mut self) { - self.submodels.clear(); - self.in_submodel.resize(self.num_subproblems(), 0); - for fidx in 0..self.num_subproblems() { - let model_idx = (self.select_model)(fidx); - if model_idx >= self.submodels.len() { - self.submodels.resize_with(model_idx + 1, SubModel::default); - } - self.submodels[model_idx].push(fidx); - self.in_submodel[fidx] = model_idx; - } - } - - /// Use a fully disaggregated model. - /// - /// A fully disaggregated model has one separate submodel for each subproblem. - /// Hence, calling this method is equivalent to - /// `CplexMaster::set_submodel_selection(|i| i)`. - pub fn use_full_disaggregation(&mut self) { - self.set_submodel_selection(|i| i) - } - - /// Use a fully aggregated model. - /// - /// A fully aggregated model has one submodel for all subproblems. - /// Hence, calling this method is equivalent to - /// `CplexMaster::set_submodel_selection(|_| 0)`. - pub fn use_full_aggregation(&mut self) { - self.set_submodel_selection(|_| 0) - } - - fn init_qp(&mut self) -> Result<()> { - if self.force_update { - self.updateinds.clear(); - for mins in &self.minorants { - self.updateinds.extend(mins.iter().map(|m| m.index)); - } - } - - let minorants = &self.minorants; - - // Compute the number of minorants in each submodel. - for submodel in self.submodels.iter_mut() { - submodel.num_mins = submodel.iter().map(|&fidx| minorants[fidx].len()).min().unwrap_or(0); - submodel.minorants.resize_with(submodel.num_mins, Minorant::default); - } - - // Only minorants belonging to the first subproblem of each submodel - // must be updated. - // - // We filter all indices that - // 1. belong to a subproblem being the first in its model - // 2. is a valid index (< minimal number of minorants within this submodel) - // and map them to (index, model_index, minorant_index) where - // - `index` is the variable index (i.e. minorant index of first subproblem) - // - `model_index` is the submodel index - // - `minorant_index` is the index of the minorant within the model - let updateinds = self - .updateinds - .iter() - .filter_map(|&index| { - let MinorantIdx { fidx, idx } = self.index2min[index]; - let mod_i = self.in_submodel[fidx]; - let submodel = &self.submodels[mod_i]; - if submodel[0] == fidx && idx < submodel.num_mins { - Some((index, mod_i, idx)) - } else { - None - } - }) - .collect::>(); - - // Compute the aggregated minorants. - for &(_, mod_i, i) in &updateinds { - let submodel = &mut self.submodels[mod_i]; - submodel.minorants[i] = - Aggregatable::combine(submodel.iter().map(|&fidx| (1.0, &minorants[fidx][i].minorant))); - } - - let submodels = &self.submodels; - - // Build quadratic term, this is for all pairs of minorants where each g_i - // is an aggregated minorant of a submodel. - // - // For simplicity we always store the terms in the index of the minorant of the first - // subproblem in each submodel. - let ntotalminorants = self.index2min.len(); - if ntotalminorants > self.qterm.len() { - self.qterm.resize(ntotalminorants, dvec![]); - for i in 0..self.qterm.len() { - self.qterm[i].resize(ntotalminorants, 0.0); - } - } - - // - i is the number of the minorant within the submodel submodel_i - // - idx_i is the unique index of that minorant (of the first subproblem) - // - j is the number of the minorant within the submodel submodel_j - // - idx_j is the unique index of that minorant (of the first subproblem) - for submodel_i in submodels.iter() { - // Compute the minorant g_i for each i - for i in 0..submodel_i.num_mins { - // Store the computed values at the index of the first subproblem in this model. - let idx_i = minorants[submodel_i[0]][i].index; - let g_i = &submodel_i.minorants[i].linear; - // Now compute the inner product with each other minorant - // that has to be updated. - for &(idx_j, mod_j, j) in updateinds.iter() { - let x = submodels[mod_j].minorants[j].linear.dot(g_i); - self.qterm[idx_i][idx_j] = x; - self.qterm[idx_j][idx_i] = x; - } - } - } - - // We verify that the qterm is correct - if cfg!(debug_assertions) { - for submod_i in submodels.iter() { - for i in 0..submod_i.num_mins { - let idx_i = self.minorants[submod_i[0]][i].index; - for submod_j in submodels.iter() { - for j in 0..submod_j.num_mins { - let idx_j = self.minorants[submod_j[0]][j].index; - let x = submod_i.minorants[i].linear.dot(&submod_j.minorants[j].linear); - debug_assert!((x - self.qterm[idx_i][idx_j]) < 1e-6); - } - } - } - } - } - - // main diagonal plus small identity to ensure Q being semi-definite - self.qdiag = 0.0; - for submodel in submodels.iter() { - for i in 0..submodel.num_mins { - let idx = minorants[submodel[0]][i].index; - self.qdiag = Real::max(self.qdiag, self.qterm[idx][idx]); - } - } - self.qdiag *= 1e-8; - - // We have updated everything. - self.updateinds.clear(); - self.force_update = false; - - self.init_cpx_qp() - } - - fn init_cpx_qp(&mut self) -> Result<()> { - if !self.lp.is_null() { - trycpx!(cpx::freeprob(self.env, &mut self.lp)); - } - trycpx!({ - let mut status = 0; - self.lp = cpx::createprob(self.env, &mut status, c_str!("mastercp").as_ptr()); - status - }); - - let nsubmodels = self.submodels.len(); - let submodels = &self.submodels; - let minorants = &self.minorants; - - // add convexity constraints - let sense: Vec = vec!['E' as c_char; nsubmodels]; - let rhs = dvec![1.0; nsubmodels]; - let mut rmatbeg = Vec::with_capacity(nsubmodels); - let mut rmatind = Vec::with_capacity(self.index2min.len()); - let mut rmatval = Vec::with_capacity(self.index2min.len()); - - let mut nvars = 0; - for submodel in submodels.iter() { - if submodel.is_empty() { - // this should only happen if the submodel selector leaves - // holes - continue; - } - - rmatbeg.push(nvars as c_int); - rmatind.extend((nvars as c_int..).take(submodel.num_mins)); - rmatval.extend(repeat(1.0).take(submodel.num_mins)); - nvars += submodel.num_mins; - } - - trycpx!(cpx::addrows( - self.env, - self.lp, - nvars as c_int, - rmatbeg.len() as c_int, - rmatind.len() as c_int, - rhs.as_ptr(), - sense.as_ptr(), - rmatbeg.as_ptr(), - rmatind.as_ptr(), - rmatval.as_ptr(), - ptr::null(), - ptr::null() - )); - - // update coefficients - let mut var_i = 0; - for (mod_i, submodel_i) in submodels.iter().enumerate() { - for i in 0..submodel_i.num_mins { - let idx_i = minorants[submodel_i[0]][i].index; - let mut var_j = 0; - for (mod_j, submodel_j) in submodels.iter().enumerate() { - for j in 0..submodel_j.num_mins { - let idx_j = minorants[submodel_j[0]][j].index; - let q = self.qterm[idx_i][idx_j] + if mod_i != mod_j || i != j { 0.0 } else { self.qdiag }; - trycpx!(cpx::chgqpcoef(self.env, self.lp, var_i as c_int, var_j as c_int, q)); - var_j += 1; - } - } - var_i += 1; - } - } - - Ok(()) - } -} - -pub struct Builder { - /// The maximal bundle size used in the master problem. - pub max_bundle_size: usize, - /// The submodel selector. - select_model: Arc usize>, -} - -impl Default for Builder { - fn default() -> Self { - Builder { - max_bundle_size: 50, - select_model: Arc::new(|i| i), - } - } -} - -impl unconstrained::Builder for Builder { - type MasterProblem = CplexMaster; - - fn build(&mut self) -> Result { - let mut cpx = CplexMaster::new()?; - cpx.max_bundle_size = self.max_bundle_size; - cpx.select_model = self.select_model.clone(); - cpx.update_submodels(); - Ok(cpx) - } -} - -impl Builder { - pub fn max_bundle_size(&mut self, s: usize) -> &mut Self { - assert!(s >= 2, "The maximal bundle size must be >= 2"); - self.max_bundle_size = s; - self - } - - /// Set a custom submodel selector. - /// - /// For each subproblem index the selector should return a submodel index. - /// All subproblems with the same submodel index are aggregated in a single - /// cutting plane model. - pub fn submodel_selection(&mut self, selector: F) -> &mut Self - where - F: Fn(usize) -> usize + 'static, - { - self.select_model = Arc::new(selector); - self - } - - /// Use a fully disaggregated model. - /// - /// A fully disaggregated model has one separate submodel for each subproblem. - /// Hence, calling this method is equivalent to - /// `CplexMaster::set_submodel_selection(|i| i)`. - pub fn use_full_disaggregation(&mut self) -> &mut Self { - self.submodel_selection(|i| i) - } - - /// Use a fully aggregated model. - /// - /// A fully aggregated model has one submodel for all subproblems. - /// Hence, calling this method is equivalent to - /// `CplexMaster::set_submodel_selection(|_| 0)`. - pub fn use_full_aggregation(&mut self) -> &mut Self { - self.submodel_selection(|_| 0) - } -} DELETED src/master/minimal.rs Index: src/master/minimal.rs ================================================================== --- src/master/minimal.rs +++ /dev/null @@ -1,402 +0,0 @@ -// Copyright (c) 2016, 2017, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -use crate::master::unconstrained; -use crate::master::{SubgradientExtension, UnconstrainedMasterProblem}; -use crate::{Aggregatable, DVector, Minorant, Real}; - -use log::debug; - -use std::error::Error; -use std::f64::NEG_INFINITY; -use std::fmt; -use std::result; - -/// Minimal master problem error. -#[derive(Debug)] -pub enum MinimalMasterError { - NoMinorants, - MaxMinorants { subproblem: usize }, - SubgradientExtension(Box), -} - -impl fmt::Display for MinimalMasterError { - fn fmt(&self, fmt: &mut fmt::Formatter) -> result::Result<(), fmt::Error> { - use self::MinimalMasterError::*; - match self { - MaxMinorants { subproblem } => write!( - fmt, - "The minimal master problem allows at most two minorants (subproblem: {})", - subproblem - ), - NoMinorants => write!(fmt, "The master problem does not contain a minorant"), - SubgradientExtension(err) => write!(fmt, "Subgradient extension failed: {}", err), - } - } -} - -impl Error for MinimalMasterError { - fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { - use MinimalMasterError::*; - match self { - SubgradientExtension(err) => Some(err.as_ref()), - _ => None, - } - } -} - -/** - * A minimal master problem with only two minorants. - * - * This is the simplest possible master problem for bundle methods. It - * has only two minorants and only one function model. The advantage - * is that this model can be solved explicitely and very quickly, but - * it is only a very loose approximation of the objective function. - * - * Because of its properties, it can only be used if the problem to be - * solved has a maximal number of minorants of two and only one - * subproblem. - */ -pub struct MinimalMaster { - /// The weight of the quadratic term. - weight: Real, - - /// The minorants in the model. - /// - /// There are up to two minorants, each minorant consists of one part for - /// each subproblem. - /// - /// The minorants for the i-th subproblem have the indices `2*i` and - /// `2*i+1`. - minorants: [Vec; 2], - /// The number of minorants. Only the minorants with index less than this - /// number are valid. - num_minorants: usize, - /// The number of minorants for each subproblem. - num_minorants_of: Vec, - /// The number of subproblems. - num_subproblems: usize, - /// The number of subproblems with at least 1 minorant. - num_subproblems_with_1: usize, - /// The number of subproblems with at least 2 minorants. - num_subproblems_with_2: usize, - /// Optimal multipliers. - opt_mult: [Real; 2], - /// Optimal aggregated minorant. - opt_minorant: Minorant, -} - -impl UnconstrainedMasterProblem for MinimalMaster { - type MinorantIndex = usize; - - type Err = MinimalMasterError; - - fn new() -> Result { - Ok(MinimalMaster { - weight: 1.0, - num_minorants: 0, - num_minorants_of: vec![], - num_subproblems: 0, - num_subproblems_with_1: 0, - num_subproblems_with_2: 0, - minorants: [vec![], vec![]], - opt_mult: [0.0, 0.0], - opt_minorant: Minorant::default(), - }) - } - - fn num_subproblems(&self) -> usize { - self.num_subproblems - } - - fn set_num_subproblems(&mut self, n: usize) -> Result<(), Self::Err> { - self.num_subproblems = n; - self.num_minorants = 0; - self.num_minorants_of = vec![0; n]; - self.num_subproblems_with_1 = 0; - self.num_subproblems_with_2 = 0; - self.minorants = [vec![Minorant::default(); n], vec![Minorant::default(); n]]; - Ok(()) - } - - fn compress(&mut self, f: F) -> Result<(), Self::Err> - where - F: FnMut(Self::MinorantIndex, &mut dyn Iterator), - { - if self.num_minorants == 2 { - debug!("Aggregate"); - debug!(" {} * {:?}", self.opt_mult[0], self.minorants[0]); - debug!(" {} * {:?}", self.opt_mult[1], self.minorants[1]); - - let mut f = f; - for fidx in 0..self.num_subproblems { - f( - 2 * fidx, - &mut self - .opt_mult - .iter() - .enumerate() - .map(|(i, alpha)| (2 * fidx + i, *alpha)), - ); - } - - self.minorants[0] = Aggregatable::combine(self.opt_mult.iter().cloned().zip(&self.minorants)); - self.opt_mult[0] = 1.0; - self.num_minorants = 1; - self.num_minorants_of.clear(); - self.num_minorants_of.resize(self.num_subproblems, 1); - self.num_subproblems_with_2 = 0; - - debug!(" {:?}", self.minorants[0]); - } - Ok(()) - } - - fn weight(&self) -> Real { - self.weight - } - - fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err> { - assert!(weight > 0.0); - self.weight = weight; - Ok(()) - } - - fn num_minorants(&self, fidx: usize) -> usize { - self.num_minorants_of[fidx] - } - - fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result { - if self.num_minorants_of[fidx] >= 2 { - return Err(MinimalMasterError::MaxMinorants { subproblem: fidx }); - } - - let minidx = self.num_minorants_of[fidx]; - self.num_minorants_of[fidx] += 1; - self.minorants[minidx][fidx] = minorant; - - match minidx { - 0 => { - self.num_subproblems_with_1 += 1; - if self.num_subproblems_with_1 == self.num_subproblems { - self.num_minorants = 1; - self.opt_mult[0] = 0.0; - } - Ok(2 * fidx) - } - 1 => { - self.num_subproblems_with_2 += 1; - if self.num_subproblems_with_2 == self.num_subproblems { - self.num_minorants = 2; - self.opt_mult[1] = 0.0; - } - Ok(2 * fidx + 1) - } - _ => unreachable!("Invalid number of minorants in subproblem {}", fidx), - } - } - - fn add_vars( - &mut self, - nnew: usize, - changed: &[usize], - extend_subgradient: &mut SubgradientExtension, - ) -> Result<(), Self::Err> { - if self.num_subproblems_with_1 == 0 { - return Ok(()); - } - - let noldvars = self.minorants[0][self.num_minorants_of.iter().position(|&n| n > 0).unwrap()] - .linear - .len(); - let mut changedvars = vec![]; - changedvars.extend_from_slice(changed); - changedvars.extend(noldvars..noldvars + nnew); - - for fidx in 0..self.num_subproblems { - for i in 0..self.num_minorants_of[fidx] { - let new_subg = extend_subgradient(fidx, 2 * fidx + i, &changedvars) - .map_err(MinimalMasterError::SubgradientExtension)?; - let m = &mut self.minorants[i][fidx]; - for (&j, &g) in changed.iter().zip(new_subg.iter()) { - m.linear[j] = g; - } - m.linear.extend_from_slice(&new_subg[changed.len()..]); - } - } - - Ok(()) - } - - #[allow(unused_variables)] - fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<(), Self::Err> { - for fidx in 0..self.num_subproblems { - for i in 0..self.num_minorants_of[fidx] { - debug!(" min(fidx:{}, i:{}) = {}", fidx, i, self.minorants[i][fidx]); - } - } - - if self.num_minorants == 2 { - let min0 = Minorant::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx]))); - let min1 = Minorant::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[1][fidx]))); - let xx = min0.linear.dot(&min0.linear); - let yy = min1.linear.dot(&min1.linear); - let xy = min0.linear.dot(&min1.linear); - let xeta = min0.linear.dot(eta); - let yeta = min1.linear.dot(eta); - let a = yy - 2.0 * xy + xx; - let b = xy - xx - yeta + xeta; - - let mut alpha2 = 0.0; - if a > 0.0 { - alpha2 = ((min1.constant - min0.constant) * self.weight - b) / a; - alpha2 = alpha2.max(0.0).min(1.0); - } - self.opt_mult[0] = 1.0 - alpha2; - self.opt_mult[1] = alpha2; - self.opt_minorant = Aggregatable::combine(self.opt_mult.iter().cloned().zip([min0, min1].iter())); - } else if self.num_minorants == 1 { - let min0 = Aggregatable::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx]))); - self.opt_minorant = min0; - self.opt_mult[0] = 1.0; - } else { - return Err(MinimalMasterError::NoMinorants); - } - - debug!("Unrestricted"); - debug!(" opt_minorant={}", self.opt_minorant); - debug!(" opt_mult={:?}", &self.opt_mult[0..self.num_minorants]); - - Ok(()) - } - - fn dualopt(&self) -> &DVector { - &self.opt_minorant.linear - } - - fn dualopt_cutval(&self) -> Real { - self.opt_minorant.constant - } - - fn multiplier(&self, min: usize) -> Real { - self.opt_mult[min % 2] - } - - fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box + 'a> { - Box::new( - self.opt_mult - .iter() - .take(self.num_minorants_of[fidx]) - .enumerate() - .map(move |(i, alpha)| (2 * fidx + i, *alpha)), - ) - } - - fn eval_model(&self, y: &DVector) -> Real { - let mut result = NEG_INFINITY; - for mins in &self.minorants[0..self.num_minorants] { - result = result.max(mins.iter().map(|m| m.eval(y)).sum()); - } - result - } - - fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector), Self::Err> { - debug!("Aggregate minorants {:?} of subproblem {}", mins, fidx); - if mins.len() == 2 { - debug_assert_ne!(mins[0], mins[1], "Minorants to be aggregated must be different"); - debug_assert_eq!( - mins[0] / 2, - fidx, - "Minorant {} does not belong to subproblem {}", - mins[0], - fidx - ); - debug_assert_eq!( - mins[1] / 2, - fidx, - "Minorant {} does not belong to subproblem {}", - mins[1], - fidx - ); - debug_assert!( - mins[0] % 2 < self.num_minorants_of[fidx], - "Invalid minorant index for subproblem {}: {}", - fidx, - mins[0] - ); - debug_assert!( - mins[1] % 2 < self.num_minorants_of[fidx], - "Invalid minorant index for subproblem {}: {}", - fidx, - mins[1] - ); - - let min0 = mins[0] % 2; - let min1 = mins[1] % 2; - - debug!("Aggregate"); - debug!(" {} * {}", self.opt_mult[min0], self.minorants[min0][fidx]); - debug!(" {} * {}", self.opt_mult[min1], self.minorants[min1][fidx]); - self.minorants[0][fidx] = Aggregatable::combine( - [ - (self.opt_mult[min0], &self.minorants[min0][fidx]), - (self.opt_mult[min1], &self.minorants[min1][fidx]), - ] - .iter() - .cloned(), - ); - - if self.num_subproblems_with_2 == self.num_subproblems { - self.num_minorants -= 1; - } - self.num_subproblems_with_2 -= 1; - self.num_minorants_of[fidx] -= 1; - - let coeffs = dvec![self.opt_mult[min0], self.opt_mult[min1]]; - - debug!(" {}", self.minorants[0][fidx]); - Ok((2 * fidx, coeffs)) - } else if mins.len() == 1 { - Ok((mins[0], dvec![1.0])) - } else { - panic!("No minorants specified to be aggregated"); - } - } - - fn move_center(&mut self, alpha: Real, d: &DVector) { - for fidx in 0..self.num_subproblems { - for i in 0..self.num_minorants_of[fidx] { - self.minorants[i][fidx].move_center(alpha, d); - } - } - } -} - -pub struct Builder; - -impl Default for Builder { - fn default() -> Self { - Builder - } -} - -impl unconstrained::Builder for Builder { - type MasterProblem = MinimalMaster; - - fn build(&mut self) -> Result { - MinimalMaster::new() - } -} Index: src/master/mod.rs ================================================================== --- src/master/mod.rs +++ src/master/mod.rs @@ -37,21 +37,10 @@ //! \hat{d})$ for some given $\hat{d} \in \mathbb{R}\^n$. pub mod boxed; pub use self::boxed::BoxedMasterProblem; -pub mod unconstrained; -pub use self::unconstrained::UnconstrainedMasterProblem; - -pub mod minimal; -pub use self::minimal::MinimalMaster; - -// pub mod grb; -// pub use master::grb::GurobiMaster; - -pub mod cpx; - pub(crate) mod primalmaster; use crate::{DVector, Minorant, Real}; use std::error::Error; use std::result::Result; @@ -175,5 +164,11 @@ type MasterProblem: MasterProblem; /// Create a new master problem. fn build(&mut self) -> Result::Err>; } + +/// The minimal bundle builder. +pub type FullMasterBuilder = boxed::Builder; + +/// The minimal bundle builder. +pub type MinimalMasterBuilder = boxed::Builder; Index: src/master/primalmaster.rs ================================================================== --- src/master/primalmaster.rs +++ src/master/primalmaster.rs @@ -16,11 +16,11 @@ */ //! A wrapper around master problems to handle primal information. use super::MasterProblem; -use crate::parallel::SubgradientExtender; +use crate::problem::SubgradientExtender; use crate::{Aggregatable, Minorant, Real}; use std::collections::HashMap; use std::ops::{Deref, DerefMut}; DELETED src/master/unconstrained.rs Index: src/master/unconstrained.rs ================================================================== --- src/master/unconstrained.rs +++ /dev/null @@ -1,139 +0,0 @@ -// Copyright (c) 2016, 2017, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -use crate::{DVector, Minorant, Real}; - -use super::SubgradientExtension; - -use std::error::Error; - -/** - * Trait for master problems without box constraints. - * - * Implementors of this trait are supposed to solve quadratic - * optimization problems of the form - * - * \\[ \min \left\\{ \hat{f}(d) + \frac{u}{2} \\| d \\|\^2 \colon - * d \in \mathbb{R}\^n \right\\}. \\] - * - * where $\hat{f}$ is a piecewise linear model, i.e. - * - * \\[ \hat{f}(d) = \max \\{ \ell_i(d) = c_i + \langle g_i, d \rangle \colon - * i=1,\dotsc,k \\} - * = \max \left\\{ \sum_{i=1}\^k \alpha_i \ell_i(d) \colon - * \alpha \in \Delta \right\\}, \\] - * - * where $\Delta := \left\\{ \alpha \in \mathbb{R}\^k \colon \sum_{i=1}\^k - * \alpha_i = 1 \right\\}$. Note, the unconstrained solver is expected - * to compute *dual* optimal solutions, i.e. the solver must compute - * optimal coefficients $\bar{\alpha}$ for the dual problem - * - * \\[ \max_{\alpha \in \Delta} \min_{d \in \mathbb{R}\^n} - * \sum_{i=1}\^k \alpha_i \ell_i(d) + \frac{u}{2} \\| d \\|\^2. \\] - */ -pub trait UnconstrainedMasterProblem: Send + 'static { - /// Unique index for a minorant. - type MinorantIndex: Copy + Eq; - - /// Error type for this master problem. - type Err: Error + Send + Sync; - - /// Return a new instance of the unconstrained master problem. - fn new() -> Result - where - Self: Sized; - - /// Return the number of subproblems. - fn num_subproblems(&self) -> usize; - - /// Set the number of subproblems (different function models.) - fn set_num_subproblems(&mut self, n: usize) -> Result<(), Self::Err>; - - /// Return the current weight. - fn weight(&self) -> Real; - - /// Set the weight of the quadratic term, must be > 0. - fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err>; - - /// Return the number of minorants of subproblem `fidx`. - fn num_minorants(&self, fidx: usize) -> usize; - - /// Compress the bundle. - /// - /// When some minorants are compressed, the callback is called with the - /// coefficients and indices of the compressed minorants and the index of - /// the new minorant. The callback may be called several times. - fn compress(&mut self, f: F) -> Result<(), Self::Err> - where - F: FnMut(Self::MinorantIndex, &mut dyn Iterator); - - /// Add a new minorant to the model. - fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result; - - /// Add or move some variables. - /// - /// The variables in `changed` have been changed, so the subgradient - /// information must be updated. Furthermore, `nnew` new variables - /// are added. - fn add_vars( - &mut self, - nnew: usize, - changed: &[usize], - extend_subgradient: &mut SubgradientExtension, - ) -> Result<(), Self::Err>; - - /// Solve the master problem. - fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<(), Self::Err>; - - /// Return the current dual optimal solution. - fn dualopt(&self) -> &DVector; - - /// Return the current dual optimal solution value. - fn dualopt_cutval(&self) -> Real; - - /// Return the multiplier associated with a minorant. - fn multiplier(&self, min: Self::MinorantIndex) -> Real; - - /// Return the multipliers associated with a subproblem. - fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box + 'a>; - - /// Return the value of the current model at the given point. - fn eval_model(&self, y: &DVector) -> Real; - - /// Aggregate the given minorants according to the current solution. - /// - /// The (indices of the) minorants to be aggregated get invalid - /// after this operation. The function returns the index of the - /// aggregated minorant along with the coefficients of the convex - /// combination. The index of the new aggregated minorant might or - /// might not be one of indices of the original minorants. - /// - /// # Error - /// The indices of the minorants `mins` must belong to subproblem `fidx`. - fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector), Self::Err>; - - /// Move the center of the master problem along $\alpha \cdot d$. - fn move_center(&mut self, alpha: Real, d: &DVector); -} - -/// A builder for creating unconstrained master problem solvers. -pub trait Builder { - /// The master problem to be build. - type MasterProblem: UnconstrainedMasterProblem; - - /// Create a new master problem. - fn build(&mut self) -> Result::Err>; -} Index: src/mcf/problem.rs ================================================================== --- src/mcf/problem.rs +++ src/mcf/problem.rs @@ -13,15 +13,15 @@ // You should have received a copy of the GNU General Public License // along with this program. If not, see // use crate::mcf; -use crate::parallel::{ +use crate::problem::{ EvalResult, FirstOrderProblem as ParallelProblem, ResultSender, Update as ParallelUpdate, UpdateSender, UpdateState as ParallelUpdateState, }; -use crate::{Aggregatable, DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation, Update, UpdateState}; +use crate::{DVector, Minorant, Real}; use itertools::izip; use log::{debug, warn}; use num_traits::Float; use threadpool::ThreadPool; @@ -183,10 +183,14 @@ rhs - lhs } } impl MMCFProblem { + pub fn num_subproblems(&self) -> usize { + self.subs.len() + } + pub fn read_mnetgen(basename: &str) -> std::result::Result { let mut buffer = String::new(); { let mut f = File::open(&format!("{}.nod", basename))?; f.read_to_string(&mut buffer)?; @@ -360,148 +364,21 @@ aggr } } -impl FirstOrderProblem for MMCFProblem { +impl ParallelProblem for MMCFProblem { type Err = Error; - type Primal = Vec; - - type EvalResult = SimpleEvaluation>; + type Primal = DVector; fn num_variables(&self) -> usize { self.active_constraints.len() } fn lower_bounds(&self) -> Option> { - Some(vec![0.0; self.active_constraints.len()]) - } - - fn upper_bounds(&self) -> Option> { - None - } - - fn num_subproblems(&self) -> usize { - if self.multimodel { - self.subs.len() - } else { - 1 - } - } - - fn evaluate(&mut self, fidx: usize, y: &[Real], _nullstep_bound: Real, _relprec: Real) -> Result { - let (objective, subg, sol) = if self.multimodel { - let (objective, subg, sol) = self.subs[fidx] - .write() - .unwrap() - .evaluate(y, self.active_constraints.iter().cloned())?; - (objective, subg, vec![sol]) - } else { - let mut objective = 0.0; - let mut subg = dvec![0.0; y.len()]; - let mut sols = Vec::with_capacity(self.subs.len()); - for sub in &mut self.subs { - let (obj, sg, sol) = sub - .write() - .unwrap() - .evaluate(y, self.active_constraints.iter().cloned())?; - objective += obj; - subg.add_scaled(1.0, &sg); - sols.push(sol); - } - (objective, subg, sols) - }; - Ok(SimpleEvaluation { - objective, - minorants: vec![( - Minorant { - constant: objective, - linear: subg, - }, - sol, - )], - }) - } - - fn update(&mut self, state: &UpdateState) -> Result> { - if self.inactive_constraints.is_empty() { - return Ok(vec![]); - } - - let nold = self.active_constraints.len(); - let subs = &self.subs; - - // if state.step != Step::Descent && !self.active_constraints.is_empty() { - // return Ok(vec![]); - // } - - let newconstraints = self - .inactive_constraints - .iter() - .map(|&cidx| { - subs.iter() - .enumerate() - .map(|(fidx, sub)| { - let primals = state.aggregated_primals(fidx); - let aggr = Aggregatable::combine(primals.into_iter().map(|(alpha, x)| (alpha, &x[0]))); - sub.read().unwrap().evaluate_constraint(&aggr, cidx) - }) - .sum::() - }) - .enumerate() - .filter_map(|(i, sg)| if sg < 1e-3 { Some(i) } else { None }) - .collect::>(); - - let inactive = &mut self.inactive_constraints; - self.active_constraints - .extend(newconstraints.into_iter().rev().map(|i| inactive.swap_remove(i))); - - // *** The following code needs `drain_filter`, which is experimental as - // of rust 1.36 *** - - // self.active_constraints - // .extend(self.inactive_constraints.drain_filter(|&mut cidx| { - // subs.iter() - // .enumerate() - // .map(|(fidx, sub)| { - // let primals = state.aggregated_primals(fidx); - // let aggr = Aggregatable::combine(primals.into_iter().map(|(alpha, x)| (alpha, &x[0]))); - // sub.read().unwrap().evaluate_constraint(&aggr, cidx) - // }) - // .sum::() < -1e-3 - // })); - - Ok(vec![ - Update::AddVariable { - lower: 0.0, - upper: Real::infinity() - }; - self.active_constraints.len() - nold - ]) - } - - fn extend_subgradient(&mut self, fidx: usize, primal: &Self::Primal, inds: &[usize]) -> Result> { - let sub = self.subs[fidx].read().unwrap(); - Ok(inds - .iter() - .map(|&i| sub.evaluate_constraint(&primal[0], self.active_constraints[i])) - .collect()) - } -} - -impl ParallelProblem for MMCFProblem { - type Err = ::Err; - - type Primal = DVector; - - fn num_variables(&self) -> usize { - FirstOrderProblem::num_variables(self) - } - - fn lower_bounds(&self) -> Option> { - FirstOrderProblem::lower_bounds(self) + Some(vec![0.0; self.num_variables()]) } fn num_subproblems(&self) -> usize { self.subs.len() } DELETED src/minorant.rs Index: src/minorant.rs ================================================================== --- src/minorant.rs +++ /dev/null @@ -1,203 +0,0 @@ -// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -//! A linear minorant. - -use crate::{DVector, Real}; - -use std::borrow::Borrow; -use std::fmt; - -/// An aggregatable object. -pub trait Aggregatable: Default { - /// Return a scaled version of `other`, i.e. `alpha * other`. - fn new_scaled(alpha: Real, other: A) -> Self - where - A: Borrow; - - /// Add a scaled version of `other` to `self`. - /// - /// This sets `self = self + alpha * other`. - fn add_scaled(&mut self, alpha: Real, other: A) - where - A: Borrow; - - /// Return $\sum\_{i=1}\^n alpha_i m_i$. - /// - /// If `aggregates` is empty return the default value. - fn combine(aggregates: I) -> Self - where - I: IntoIterator, - A: Borrow, - { - let mut it = aggregates.into_iter(); - let mut x; - if let Some((alpha, y)) = it.next() { - x = Self::new_scaled(alpha, y); - } else { - return Self::default(); - } - - for (alpha, y) in it { - x.add_scaled(alpha, y); - } - - x - } -} - -/// Implement for empty tuples. -impl Aggregatable for () { - fn new_scaled(_alpha: Real, _other: A) -> Self - where - A: Borrow, - { - } - - fn add_scaled(&mut self, _alpha: Real, _other: A) - where - A: Borrow, - { - } -} - -impl Aggregatable for Real { - fn new_scaled(alpha: Real, other: A) -> Self - where - A: Borrow, - { - alpha * other.borrow() - } - - fn add_scaled(&mut self, alpha: Real, other: A) - where - A: Borrow, - { - *self += alpha * other.borrow() - } -} - -/** - * A linear minorant of a convex function. - * - * A linear minorant of a convex function $f \colon \mathbb{R}\^n \to - * \mathbb{R}$ is a linear function of the form - * - * \\[ l \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto \langle g, x - * \rangle + c \\] - * - * such that $l(x) \le f(x)$ for all $x \in \mathbb{R}\^n$. - */ -#[derive(Clone, Debug)] -pub struct Minorant { - /// The constant term. - pub constant: Real, - - /// The linear term. - pub linear: DVector, -} - -impl fmt::Display for Minorant { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "{} + y * {}", self.constant, self.linear)?; - Ok(()) - } -} - -impl Default for Minorant { - fn default() -> Minorant { - Minorant { - constant: 0.0, - linear: dvec![], - } - } -} - -impl Minorant { - /// Return a new 0 minorant. - pub fn new(constant: Real, linear: Vec) -> Minorant { - Minorant { - constant, - linear: DVector(linear), - } - } - - /** - * Evaluate minorant at some point. - * - * This function computes $c + \langle g, x \rangle$ for this minorant - * \\[\ell \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto c + \langle g, x \rangle\\] - * and the given point $x \in \mathbb{R}\^n$. - */ - pub fn eval(&self, x: &DVector) -> Real { - self.constant + self.linear.dot(x) - } - - /** - * Move the center of the minorant. - */ - pub fn move_center(&mut self, alpha: Real, d: &DVector) { - self.constant += alpha * self.linear.dot(d); - } -} - -impl Aggregatable for Minorant { - fn new_scaled(alpha: Real, other: A) -> Self - where - A: Borrow, - { - let m = other.borrow(); - Minorant { - constant: alpha * m.constant, - linear: DVector::scaled(&m.linear, alpha), - } - } - - fn add_scaled(&mut self, alpha: Real, other: A) - where - A: Borrow, - { - let m = other.borrow(); - self.constant += alpha * m.constant; - self.linear.add_scaled(alpha, &m.linear); - } -} - -impl Aggregatable for Vec -where - T: Aggregatable, -{ - fn new_scaled(alpha: Real, other: A) -> Self - where - A: std::borrow::Borrow, - { - other - .borrow() - .iter() - .map(|y| Aggregatable::new_scaled(alpha, y)) - .collect() - } - - fn add_scaled(&mut self, alpha: Real, other: A) - where - A: std::borrow::Borrow, - { - debug_assert_eq!(self.len(), other.borrow().len(), "Vectors must have the same size"); - for (ref mut x, y) in self.iter_mut().zip(other.borrow()) { - x.add_scaled(alpha, y) - } - } -} DELETED src/parallel/masterprocess.rs Index: src/parallel/masterprocess.rs ================================================================== --- src/parallel/masterprocess.rs +++ /dev/null @@ -1,274 +0,0 @@ -/* - * Copyright (c) 2019 Frank Fischer - * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see - */ - -//! Asynchronous process solving a master problem. - -use crossbeam::channel::{unbounded as channel, Receiver, Sender}; -use log::{debug, warn}; -use std::sync::Arc; -use threadpool::ThreadPool; - -use super::problem::{FirstOrderProblem, SubgradientExtender}; -use super::solver::Error; -use crate::master::primalmaster::PrimalMaster; -use crate::master::MasterProblem; -use crate::{DVector, Minorant, Real}; - -/// Configuration information for setting up a master problem. -pub struct MasterConfig { - /// The number of subproblems. - pub num_subproblems: usize, - /// The number of variables. - pub num_vars: usize, - /// The lower bounds on the variables. - pub lower_bounds: Option, - /// The lower bounds on the variables. - pub upper_bounds: Option, -} - -/// A task for the master problem. -enum MasterTask -where - M: MasterProblem, -{ - /// Add new variables to the master problem. - AddVariables(Vec<(Option, Real, Real)>, Box>), - - /// Add a new minorant for a subfunction to the master problem. - AddMinorant(usize, Minorant, Pr), - - /// Move the center of the master problem in the given direction. - MoveCenter(Real, Arc), - - /// Start a new computation of the master problem. - Solve { center_value: Real }, - - /// Compress the bundle. - Compress, - - /// Set the weight parameter of the master problem. - SetWeight { weight: Real }, - - /// Return the current aggregated primal. - GetAggregatedPrimal { - subproblem: usize, - tx: Sender>, - }, -} - -/// The response send from a master process. -/// -/// The response contains the evaluation results of the latest -pub struct MasterResponse { - pub nxt_d: DVector, - pub nxt_mod: Real, - pub sgnorm: Real, - /// The number of internal iterations. - pub cnt_updates: usize, -} - -type ToMasterSender = Sender::Primal,

::Err, M>>; - -type ToMasterReceiver = Receiver::Primal,

::Err, M>>; - -type MasterSender = Sender>; - -pub type MasterReceiver = Receiver>; - -pub struct MasterProcess -where - P: FirstOrderProblem, - M: MasterProblem, -{ - /// The channel to transmit new tasks to the master problem. - tx: ToMasterSender, - - /// The channel to receive solutions from the master problem. - pub rx: MasterReceiver, - - phantom: std::marker::PhantomData, -} - -impl MasterProcess -where - P: FirstOrderProblem, - P::Primal: Send + 'static, - P::Err: Into> + 'static, - M: MasterProblem + Send + 'static, - M::MinorantIndex: std::hash::Hash, - M::Err: Send + 'static, -{ - pub fn start(master: M, master_config: MasterConfig, threadpool: &mut ThreadPool) -> Self { - // Create a pair of communication channels. - let (to_master_tx, to_master_rx) = channel(); - let (from_master_tx, from_master_rx) = channel(); - - // The the master process thread. - threadpool.execute(move || { - debug!("Master process started"); - let mut from_master_tx = from_master_tx; - if let Err(err) = Self::master_main(master, master_config, &mut from_master_tx, to_master_rx) { - #[allow(unused_must_use)] - { - from_master_tx.send(Err(err)); - } - } - debug!("Master proces stopped"); - }); - - MasterProcess { - tx: to_master_tx, - rx: from_master_rx, - phantom: std::marker::PhantomData, - } - } - - /// Add new variables to the master problem. - pub fn add_vars( - &mut self, - vars: Vec<(Option, Real, Real)>, - sgext: Box>, - ) -> Result<(), Error> - where - P::Err: 'static, - { - self.tx - .send(MasterTask::AddVariables(vars, sgext)) - .map_err(|err| Error::Process(err.into())) - } - - /// Add a new minorant to the master problem model. - /// - /// This adds the specified `minorant` with associated `primal` data to the - /// model of subproblem `i`. - pub fn add_minorant(&mut self, i: usize, minorant: Minorant, primal: P::Primal) -> Result<(), Error> { - self.tx - .send(MasterTask::AddMinorant(i, minorant, primal)) - .map_err(|err| Error::Process(err.into())) - } - - /// Move the center of the master problem. - /// - /// This moves the master problem's center in direction $\\alpha \\cdot d$. - pub fn move_center(&mut self, alpha: Real, d: Arc) -> Result<(), Error> { - self.tx - .send(MasterTask::MoveCenter(alpha, d)) - .map_err(|err| Error::Process(err.into())) - } - - /// Solve the master problem. - /// - /// The current function value in the center `center_value`. - /// Once the master problem is solved the process will send a - /// [`MasterResponse`] message to the `tx` channel. - pub fn solve(&mut self, center_value: Real) -> Result<(), Error> { - self.tx - .send(MasterTask::Solve { center_value }) - .map_err(|err| Error::Process(err.into())) - } - - /// Compresses the model. - pub fn compress(&mut self) -> Result<(), Error> { - self.tx - .send(MasterTask::Compress) - .map_err(|err| Error::Process(err.into())) - } - - /// Sets the new weight of the proximal term in the master problem. - pub fn set_weight(&mut self, weight: Real) -> Result<(), Error> { - self.tx - .send(MasterTask::SetWeight { weight }) - .map_err(|err| Error::Process(err.into())) - } - - /// Get the current aggregated primal for a certain subproblem. - pub fn get_aggregated_primal(&self, subproblem: usize) -> Result> { - let (tx, rx) = channel(); - self.tx - .send(MasterTask::GetAggregatedPrimal { subproblem, tx }) - .map_err(|err| Error::Process(err.into()))?; - rx.recv() - .map_err(|err| Error::Process(err.into()))? - .map_err(|err| Error::Master(err.into())) - } - - /// The main loop of the master process. - fn master_main( - master: M, - master_config: MasterConfig, - tx: &mut MasterSender, - rx: ToMasterReceiver, - ) -> Result<(), M::Err> { - let mut master = PrimalMaster::<_, P::Primal>::new(master); - - // Initialize the master problem. - master.set_num_subproblems(master_config.num_subproblems)?; - master.set_vars( - master_config.num_vars, - master_config.lower_bounds, - master_config.upper_bounds, - )?; - - // The main iteration: wait for new tasks. - for m in rx { - match m { - MasterTask::AddVariables(vars, sgext) => { - debug!("master: add {} variables to the subproblem", vars.len()); - master.add_vars(vars, sgext)?; - } - MasterTask::AddMinorant(i, m, primal) => { - debug!("master: add minorant to subproblem {}", i); - master.add_minorant(i, m, primal)?; - } - MasterTask::MoveCenter(alpha, d) => { - debug!("master: move center"); - master.move_center(alpha, &d); - } - MasterTask::Compress => { - debug!("Compress bundle"); - master.compress()?; - } - MasterTask::Solve { center_value } => { - debug!("master: solve with center_value {}", center_value); - master.solve(center_value)?; - let master_response = MasterResponse { - nxt_d: master.get_primopt(), - nxt_mod: master.get_primoptval(), - sgnorm: master.get_dualoptnorm2().sqrt(), - cnt_updates: master.cnt_updates(), - }; - if let Err(err) = tx.send(Ok(master_response)) { - warn!("Master process cancelled because of channel error: {}", err); - break; - } - } - MasterTask::SetWeight { weight } => { - debug!("master: set weight {}", weight); - master.set_weight(weight)?; - } - MasterTask::GetAggregatedPrimal { subproblem, tx } => { - debug!("master: get aggregated primal for {}", subproblem); - if tx.send(master.aggregated_primal(subproblem)).is_err() { - warn!("Sending of aggregated primal for {} failed", subproblem); - }; - } - }; - } - - Ok(()) - } -} DELETED src/parallel/mod.rs Index: src/parallel/mod.rs ================================================================== --- src/parallel/mod.rs +++ /dev/null @@ -1,36 +0,0 @@ -/* - * Copyright (c) 2019 Frank Fischer - * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see - */ - -//! Implementation of a asynchronous parallel proximal bundle method. - -mod problem; -pub use self::problem::{ - EvalResult, FirstOrderProblem, ResultSender, SubgradientExtender, Update, UpdateSender, UpdateState, -}; - -mod solver; -pub use self::solver::Solver; - -/// The default bundle solver with general master problem. -pub type DefaultSolver

= - Solver; - -/// A bundle solver with a minimal cutting plane model. -pub type NoBundleSolver

= - Solver; - -mod masterprocess; DELETED src/parallel/problem.rs Index: src/parallel/problem.rs ================================================================== --- src/parallel/problem.rs +++ /dev/null @@ -1,190 +0,0 @@ -/* - * Copyright (c) 2019 Frank Fischer - * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see - */ - -//! An asynchronous first-order oracle. - -use crate::{Aggregatable, DVector, Minorant, Real}; -use crossbeam::channel::Sender; -use std::sync::Arc; - -/// Evaluation result. -/// -/// The result of an evaluation is new information to be made -/// available to the solver and the master problem. There are -/// essentially two types of information: -/// -/// 1. The (exact) function value of a sub-function at some point. -/// 2. A minorant of some sub-function. -#[derive(Debug)] -pub enum EvalResult { - /// The objective value at some point. - ObjectiveValue { index: I, value: Real }, - /// A minorant with an associated primal. - Minorant { index: I, minorant: Minorant, primal: P }, -} - -/// Channel to send evaluation results to. -pub type ResultSender = Sender, E>>; - -/// Problem update information. -/// -/// The solver calls the `update` method of the problem regularly. -/// This method can modify the problem by adding (or moving) -/// variables. The possible updates are encoded in this type. -pub enum Update { - /// Add new variables with bounds. - /// - /// The initial value of each variable will be the feasible value - /// closest to 0. - AddVariables { - index: I, - bounds: Vec<(Real, Real)>, - sgext: Box>, - }, -} - -/// The subgradient extender is a callback used to update existing minorants -/// given their associated primal data. -pub type SubgradientExtender = dyn FnMut(usize, &P, &[usize]) -> Result + Send; - -/// This trait provides information available in the -/// [`FirstOrderProblem::update`] method. -pub trait UpdateState

{ - /// Whether the last step was a descent step. - fn was_descent(&self) -> bool; - - /// Whether the last step was a null step. - fn was_null(&self) -> bool { - !self.was_descent() - } - - /// The (old) current center of stability. - fn center(&self) -> Arc; - - /// The candidate point. - /// - /// After a descent step, i.e. if [`UpdateState::was_descent`] is `true`, - /// this is the new center of stability. - fn candidate(&self) -> Arc; - - /// The current aggregated primal information. - /// - /// Return the aggregated primal information for the given subproblem. - fn aggregated_primal(&self, i: usize) -> P; -} - -/// Channel to send problem updates to. -pub type UpdateSender = Sender, E>>; - -/// Trait for implementing a first-order problem description. -/// -/// All computations made by an implementation are supposed to -/// be asynchronous. Hence, the interface is slightly different -/// compared with [`crate::FirstOrderProblem`]. -pub trait FirstOrderProblem { - /// Error raised by this oracle. - type Err; - - /// The primal information associated with a minorant. - type Primal: Aggregatable + Send + 'static; - - /// Return the number of variables. - fn num_variables(&self) -> usize; - - /// Return the lower bounds on the variables. - /// - /// If no lower bounds a specified, $-\infty$ is assumed. - /// - /// The lower bounds must be less then or equal the upper bounds. - fn lower_bounds(&self) -> Option> { - None - } - - /** - * Return the upper bounds on the variables. - * - * If no lower bounds a specified, $+\infty$ is assumed. - * - * The upper bounds must be greater than or equal the upper bounds. - */ - fn upper_bounds(&self) -> Option> { - None - } - - /// Return the number of subproblems. - fn num_subproblems(&self) -> usize { - 1 - } - - /// Start background processes. - /// - /// This method is called right before the solver starts the solution process. - /// It can be used to setup any background tasks required for the evaluation - /// of the subfunctions. - /// - /// Remember that background processes should be cleanup when the problem - /// is deleted (e.g. by implementing the [`Drop`] trait). - /// - /// The default implementation does nothing. - fn start(&mut self) {} - - /// Stop background processes. - /// - /// This method is called right after the solver stops the solution process. - /// It can be used to stop any background tasks required for the evaluation - /// of the subfunctions. - /// - /// A correct implementation of should cleanup all processes from the [`Drop`] - /// thread. - /// - /// The default implementation does nothing. - fn stop(&mut self) {} - - /// Start the evaluation of the i^th subproblem at the given point. - /// - /// The results of the evaluation should be passed to the provided channel. - /// In order to work correctly, the results must contain (an upper bound on) - /// the objective value at $y$ as well as at least one subgradient centered - /// at $y$ eventually. - fn evaluate( - &mut self, - i: usize, - y: Arc, - index: I, - tx: ResultSender, - ) -> Result<(), Self::Err>; - - /// Called to update the problem. - /// - /// This method is called regularly by the solver. The problem should send problem update - /// information (e.g. adding new variables) to the provided channel. - /// - /// The updates might be generated asynchronously. - /// - /// The default implementation does nothing. - fn update( - &mut self, - _state: &U, - _index: I, - _tx: UpdateSender, - ) -> Result<(), Self::Err> - where - U: UpdateState, - { - Ok(()) - } -} DELETED src/parallel/solver.rs Index: src/parallel/solver.rs ================================================================== --- src/parallel/solver.rs +++ /dev/null @@ -1,856 +0,0 @@ -/* - * Copyright (c) 2019 Frank Fischer - * - * This program is free software: you can redistribute it and/or - * modify it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 3 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, but - * WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - * General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see - */ - -//! An asynchronous parallel bundle solver. - -use crossbeam::channel::{select, unbounded as channel, Receiver, Sender}; -use log::{debug, info}; -use num_cpus; -use num_traits::Float; -use std::sync::Arc; -use std::time::Instant; -use threadpool::ThreadPool; - -use crate::{DVector, Real}; - -use super::masterprocess::{MasterConfig, MasterProcess, MasterResponse}; -use super::problem::{EvalResult, FirstOrderProblem, Update, UpdateState}; -use crate::master::{self, MasterProblem}; -use crate::solver::{SolverParams, Step}; -use crate::terminator::{StandardTerminatable, StandardTerminator, Terminator}; -use crate::weighter::{HKWeightable, HKWeighter, Weighter}; - -/// The default iteration limit. -pub const DEFAULT_ITERATION_LIMIT: usize = 10_000; - -/// Error raised by the parallel bundle [`Solver`]. -#[derive(Debug)] -pub enum Error { - /// An error raised when creating a new master problem solver. - BuildMaster(Box), - /// An error raised by the master problem process. - Master(Box), - /// The iteration limit has been reached. - IterationLimit { limit: usize }, - /// An error raised by a subproblem evaluation. - Evaluation(E), - /// An error raised subproblem update. - Update(E), - /// The dimension of some data is wrong. - Dimension(String), - /// Invalid bounds for a variable. - InvalidBounds { lower: Real, upper: Real }, - /// The value of a variable is outside its bounds. - ViolatedBounds { lower: Real, upper: Real, value: Real }, - /// The variable index is out of bounds. - InvalidVariable { index: usize, nvars: usize }, - /// An error occurred in a subprocess. - Process(Box), - /// A method requiring an initialized solver has been called. - NotInitialized, - /// The problem has not been solved yet. - NotSolved, -} - -impl std::fmt::Display for Error -where - E: std::fmt::Display, -{ - fn fmt(&self, fmt: &mut std::fmt::Formatter) -> std::result::Result<(), std::fmt::Error> { - use Error::*; - match self { - BuildMaster(err) => writeln!(fmt, "Cannot create master problem solver: {}", err), - Master(err) => writeln!(fmt, "Error in master problem: {}", err), - IterationLimit { limit } => writeln!(fmt, "The iteration limit has been reached: {}", limit), - Evaluation(err) => writeln!(fmt, "Error in subproblem evaluation: {}", err), - Update(err) => writeln!(fmt, "Error in subproblem update: {}", err), - Dimension(what) => writeln!(fmt, "Wrong dimension for {}", what), - InvalidBounds { lower, upper } => write!(fmt, "Invalid bounds, lower:{}, upper:{}", lower, upper), - ViolatedBounds { lower, upper, value } => write!( - fmt, - "Violated bounds, lower:{}, upper:{}, value:{}", - lower, upper, value - ), - InvalidVariable { index, nvars } => { - write!(fmt, "Variable index out of bounds, got:{} must be < {}", index, nvars) - } - Process(err) => writeln!(fmt, "Error in subprocess: {}", err), - NotInitialized => writeln!(fmt, "The solver must be initialized (called Solver::init()?)"), - NotSolved => writeln!(fmt, "The problem has not been solved yet"), - } - } -} - -impl std::error::Error for Error -where - E: std::error::Error + 'static, -{ - fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { - use Error::*; - match self { - BuildMaster(err) => Some(err.as_ref()), - Master(err) => Some(err.as_ref()), - Evaluation(err) => Some(err), - Process(err) => Some(err.as_ref()), - _ => None, - } - } -} - -type ClientSender

= - Sender::Primal>,

::Err>>; - -type ClientReceiver

= - Receiver::Primal>,

::Err>>; - -/// Parameters for tuning the solver. -pub type Parameters = SolverParams; - -pub struct SolverData { - /// Current center of stability. - cur_y: DVector, - - /// Function value in the current point. - cur_val: Real, - - /// Function value at the current candidate. - nxt_val: Real, - - /// Model value at the current candidate. - nxt_mod: Real, - - /// The value of the new minorant in the current center. - new_cutval: Real, - - /// The current expected progress. - /// - /// This value is actually `cur_val - nxt_val`. We store it separately only - /// for debugging purposes because after a descent step `cur_val` will be - /// changed and we could not see the "old" expected progress anymore that - /// led to the descent step. - expected_progress: Real, - - /// Norm of current aggregated subgradient. - sgnorm: Real, - - /// The currently used master problem weight. - cur_weight: Real, -} - -impl SolverData { - /// Reset solver data to initial values. - /// - /// This means that almost everything is set to +infinity so that - /// a null-step is forced after the first evaluation. - fn init(&mut self, y: DVector) { - self.cur_y = y; - self.cur_val = Real::infinity(); - self.nxt_val = Real::infinity(); - self.nxt_mod = -Real::infinity(); - self.new_cutval = -Real::infinity(); - self.expected_progress = Real::infinity(); - self.sgnorm = Real::infinity(); - self.cur_weight = 1.0; - } -} - -impl StandardTerminatable for SolverData { - fn center_value(&self) -> Real { - self.cur_val - } - - fn expected_progress(&self) -> Real { - self.expected_progress - } -} - -impl HKWeightable for SolverData { - fn current_weight(&self) -> Real { - self.cur_weight - } - - fn center(&self) -> &DVector { - &self.cur_y - } - - fn center_value(&self) -> Real { - self.cur_val - } - - fn candidate_value(&self) -> Real { - self.nxt_val - } - - fn candidate_model(&self) -> Real { - self.nxt_mod - } - - fn new_cutvalue(&self) -> Real { - self.new_cutval - } - - fn sgnorm(&self) -> Real { - self.sgnorm - } -} - -/// Internal data used during the main iteration loop. -struct IterData { - /// Maximal number of iterations. - max_iter: usize, - cnt_iter: usize, - cnt_updates: usize, - nxt_ubs: Vec, - cnt_remaining_ubs: usize, - nxt_cutvals: Vec, - cnt_remaining_mins: usize, - nxt_d: Arc, - nxt_y: Arc, - /// True if the problem has been updated after the last evaluation. - updated: bool, -} - -impl IterData { - fn new(num_subproblems: usize, num_variables: usize, max_iter: usize) -> Self { - IterData { - max_iter, - cnt_iter: 0, - cnt_updates: 0, - nxt_ubs: vec![Real::infinity(); num_subproblems], - cnt_remaining_ubs: num_subproblems, - nxt_cutvals: vec![-Real::infinity(); num_subproblems], - cnt_remaining_mins: num_subproblems, - nxt_d: Arc::new(dvec![0.0; num_variables]), - nxt_y: Arc::new(dvec![]), - updated: true, - } - } -} - -/// Data providing access for updating the problem. -struct UpdateData<'a, P, M> -where - P: FirstOrderProblem, - M: MasterProblem, -{ - /// Type of step. - step: Step, - - /// Current center of stability. - cur_y: &'a DVector, - - /// Current candidate. - nxt_y: &'a Arc, - - /// The master process. - master_proc: &'a MasterProcess, -} - -impl<'a, P, M> UpdateState for UpdateData<'a, P, M> -where - P: FirstOrderProblem, - P::Err: Into> + 'static, - M: MasterProblem, - M::MinorantIndex: std::hash::Hash, -{ - fn was_descent(&self) -> bool { - self.step == Step::Descent - } - - fn center(&self) -> Arc { - Arc::new(self.cur_y.clone()) - } - - fn candidate(&self) -> Arc { - self.nxt_y.clone() - } - - fn aggregated_primal(&self, i: usize) -> P::Primal { - self.master_proc - .get_aggregated_primal(i) - .map_err(|_| "get_aggregated_primal".to_string()) - .expect("Cannot get aggregated primal from master process") - } -} - -/// Implementation of a parallel bundle method. -pub struct Solver -where - P: FirstOrderProblem, - M: master::Builder, -{ - /// Parameters for the solver. - pub params: Parameters, - - /// Termination predicate. - pub terminator: T, - - /// Weighter heuristic. - pub weighter: W, - - /// The threadpool of the solver. - pub threadpool: ThreadPool, - - /// The master problem builder. - pub master: M, - - /// The first order problem. - problem: P, - - /// The algorithm data. - data: SolverData, - - /// The master problem process. - master_proc: Option>, - - /// The channel to receive the evaluation results from subproblems. - client_tx: Option>, - - /// The channel to receive the evaluation results from subproblems. - client_rx: Option>, - - /// Number of descent steps. - cnt_descent: usize, - - /// Number of null steps. - cnt_null: usize, - - /// Number of function evaluation. - cnt_evals: usize, - - /// Time when the solution process started. - /// - /// This is actually the time of the last call to `Solver::init`. - start_time: Instant, -} - -impl Solver -where - P: FirstOrderProblem, - P::Err: Into> + 'static, - T: Terminator + Default, - W: Weighter + Default, - M: master::Builder, - M::MasterProblem: MasterProblem, - ::MinorantIndex: std::hash::Hash, -{ - /// Create a new parallel bundle solver. - pub fn new(problem: P) -> Self - where - M: Default, - { - Solver { - params: Parameters::default(), - terminator: Default::default(), - weighter: Default::default(), - problem, - data: SolverData { - cur_y: dvec![], - cur_val: 0.0, - nxt_val: 0.0, - nxt_mod: 0.0, - new_cutval: 0.0, - expected_progress: 0.0, - sgnorm: 0.0, - cur_weight: 1.0, - }, - - threadpool: ThreadPool::with_name("Parallel bundle solver".to_string(), num_cpus::get()), - master: M::default(), - master_proc: None, - client_tx: None, - client_rx: None, - - cnt_descent: 0, - cnt_null: 0, - cnt_evals: 0, - - start_time: Instant::now(), - } - } - - /// Create a new parallel bundle solver. - pub fn with_master(problem: P, master: M) -> Self { - Solver { - params: Parameters::default(), - terminator: Default::default(), - weighter: Default::default(), - problem, - data: SolverData { - cur_y: dvec![], - cur_val: 0.0, - nxt_val: 0.0, - nxt_mod: 0.0, - new_cutval: 0.0, - expected_progress: 0.0, - sgnorm: 0.0, - cur_weight: 1.0, - }, - - threadpool: ThreadPool::with_name("Parallel bundle solver".to_string(), num_cpus::get()), - master, - master_proc: None, - client_tx: None, - client_rx: None, - - cnt_descent: 0, - cnt_null: 0, - cnt_evals: 0, - - start_time: Instant::now(), - } - } - - /// Return the underlying threadpool. - /// - /// In order to use the same threadpool for concurrent processes, - /// just clone the returned `ThreadPool`. - pub fn threadpool(&self) -> &ThreadPool { - &self.threadpool - } - - /// Set the threadpool. - /// - /// This function allows to use a specific threadpool for all processes - /// spawned by the solver. Note that this does not involve any threads - /// used by the problem because the solver is not responsible for executing - /// the evaluation process of the subproblems. However, the problem might - /// use the same threadpool as the solver. - pub fn set_threadpool(&mut self, threadpool: ThreadPool) { - self.threadpool = threadpool; - } - - /// Return the current problem associated with the solver. - pub fn problem(&self) -> &P { - &self.problem - } - - /// Initialize the solver. - /// - /// This will reset the internal data structures so that a new fresh - /// solution process can be started. - /// - /// It will also setup all worker processes. - /// - /// This function is automatically called by [`Solver::solve`]. - pub fn init(&mut self) -> Result<(), Error> { - debug!("Initialize solver"); - - let n = self.problem.num_variables(); - let m = self.problem.num_subproblems(); - - self.data.init(dvec![0.0; n]); - self.cnt_descent = 0; - self.cnt_null = 0; - self.cnt_evals = 0; - - let (tx, rx) = channel(); - self.client_tx = Some(tx); - self.client_rx = Some(rx); - - let master_config = MasterConfig { - num_subproblems: m, - num_vars: n, - lower_bounds: self.problem.lower_bounds().map(DVector), - upper_bounds: self.problem.upper_bounds().map(DVector), - }; - - if master_config - .lower_bounds - .as_ref() - .map(|lb| lb.len() != n) - .unwrap_or(false) - { - return Err(Error::Dimension("lower bounds".to_string())); - } - if master_config - .upper_bounds - .as_ref() - .map(|ub| ub.len() != n) - .unwrap_or(false) - { - return Err(Error::Dimension("upper bounds".to_string())); - } - - debug!("Start master process"); - self.master_proc = Some(MasterProcess::start( - self.master.build().map_err(|err| Error::BuildMaster(err.into()))?, - master_config, - &mut self.threadpool, - )); - - debug!("Initial problem evaluation"); - // We need an initial evaluation of all oracles for the first center. - let y = Arc::new(self.data.cur_y.clone()); - for i in 0..m { - self.problem - .evaluate(i, y.clone(), i, self.client_tx.clone().unwrap()) - .map_err(Error::Evaluation)?; - } - - debug!("Initialization complete"); - - self.start_time = Instant::now(); - - Ok(()) - } - - /// Solve the problem with the default maximal iteration limit [`DEFAULT_ITERATION_LIMIT`]. - pub fn solve(&mut self) -> Result<(), Error> { - self.solve_with_limit(DEFAULT_ITERATION_LIMIT) - } - - /// Solve the problem with a maximal iteration limit. - pub fn solve_with_limit(&mut self, limit: usize) -> Result<(), Error> { - // First initialize the internal data structures. - self.init()?; - - if self.solve_iter(limit)? { - Ok(()) - } else { - Err(Error::IterationLimit { limit }) - } - } - - /// Solve the problem but stop after at most `niter` iterations. - /// - /// The function returns `Ok(true)` if the termination criterion - /// has been satisfied. Otherwise it returns `Ok(false)` or an - /// error code. - /// - /// If this function is called again, the solution process is - /// continued from the previous point. Because of this one *must* - /// call `init()` before the first call to this function. - pub fn solve_iter(&mut self, niter: usize) -> Result> { - debug!("Start solving up to {} iterations", niter); - - let mut itdata = IterData::new(self.problem.num_subproblems(), self.problem.num_variables(), niter); - - loop { - select! { - recv(self.client_rx.as_ref().ok_or(Error::NotInitialized)?) -> msg => { - let msg = msg - .map_err(|err| Error::Process(err.into()))? - .map_err(Error::Evaluation)?; - if self.handle_client_response(msg, &mut itdata)? { - return Ok(false); - } - }, - recv(self.master_proc.as_ref().ok_or(Error::NotInitialized)?.rx) -> msg => { - debug!("Receive master response"); - // Receive result (new candidate) from the master - let master_res = msg - .map_err(|err| Error::Process(err.into()))? - .map_err(|err| Error::Master(err.into()))?; - - if self.handle_master_response(master_res, &mut itdata)? { - return Ok(true); - } - }, - } - } - } - - /// Handle a response from a subproblem evaluation. - /// - /// The function returns `Ok(true)` if the final iteration count has been reached. - fn handle_client_response( - &mut self, - msg: EvalResult::Primal>, - itdata: &mut IterData, - ) -> Result> { - let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; - match msg { - EvalResult::ObjectiveValue { index, value } => { - debug!("Receive objective from subproblem {}: {}", index, value); - if itdata.nxt_ubs[index].is_infinite() { - itdata.cnt_remaining_ubs -= 1; - } - itdata.nxt_ubs[index] = itdata.nxt_ubs[index].min(value); - } - EvalResult::Minorant { - index, - mut minorant, - primal, - } => { - debug!("Receive minorant from subproblem {}", index); - if itdata.nxt_cutvals[index].is_infinite() { - itdata.cnt_remaining_mins -= 1; - } - // move center of minorant to cur_y - minorant.move_center(-1.0, &itdata.nxt_d); - itdata.nxt_cutvals[index] = itdata.nxt_cutvals[index].max(minorant.constant); - // add minorant to master problem - master.add_minorant(index, minorant, primal)?; - } - } - - if itdata.cnt_remaining_ubs > 0 || itdata.cnt_remaining_mins > 0 { - // Haven't received data from all subproblems, yet. - return Ok(false); - } - - // All subproblems have been evaluated, do a step. - let nxt_ub = itdata.nxt_ubs.iter().sum::(); - let descent_bnd = Self::get_descent_bound(self.params.acceptance_factor, &self.data); - - self.data.nxt_val = nxt_ub; - self.data.new_cutval = itdata.nxt_cutvals.iter().sum::(); - - debug!("Step"); - debug!(" cur_val ={}", self.data.cur_val); - debug!(" nxt_mod ={}", self.data.nxt_mod); - debug!(" nxt_ub ={}", nxt_ub); - debug!(" descent_bnd={}", descent_bnd); - - itdata.updated = false; - let step; - if self.data.cur_val.is_infinite() { - // This is the first evaluation. We effectively get - // the function value at the current center but - // we do not have a model estimate yet. Hence, we do not know - // a good guess for the weight. - step = Step::Descent; - self.data.cur_val = nxt_ub; - self.data.cur_weight = Real::infinity(); - master.set_weight(1.0)?; - - itdata.updated = true; - - debug!("First Step"); - debug!(" cur_val={}", self.data.cur_val); - debug!(" cur_y={}", self.data.cur_y); - } else if nxt_ub <= descent_bnd { - step = Step::Descent; - self.cnt_descent += 1; - - // Note that we must update the weight *before* we - // change the internal data, so the old information - // that caused the descent step is still available. - self.data.cur_weight = self.weighter.descent_weight(&self.data); - self.data.cur_y = itdata.nxt_y.as_ref().clone(); - self.data.cur_val = nxt_ub; - - master.move_center(1.0, itdata.nxt_d.clone())?; - master.set_weight(self.data.cur_weight)?; - - debug!("Descent Step"); - debug!(" dir ={}", itdata.nxt_d); - debug!(" newy={}", self.data.cur_y); - } else { - step = Step::Null; - self.cnt_null += 1; - self.data.cur_weight = self.weighter.null_weight(&self.data); - master.set_weight(self.data.cur_weight)?; - } - - Self::show_info( - &self.start_time, - step, - &self.data, - self.cnt_descent, - self.cnt_null, - itdata.cnt_updates, - ); - itdata.cnt_iter += 1; - - // Update problem. - if Self::update_problem(&mut self.problem, step, &mut self.data, itdata, master)? { - itdata.updated = true; - } - - // Compute the new candidate. The main loop will wait for the result of - // this solution process of the master problem. - master.solve(self.data.cur_val)?; - - Ok(itdata.cnt_iter >= itdata.max_iter) - } - - fn handle_master_response( - &mut self, - master_res: MasterResponse, - itdata: &mut IterData, - ) -> Result> { - let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; - - self.data.nxt_mod = master_res.nxt_mod; - self.data.sgnorm = master_res.sgnorm; - self.data.expected_progress = self.data.cur_val - self.data.nxt_mod; - itdata.cnt_updates = master_res.cnt_updates; - - // If this is the very first solution of the model, - // we use its result as to make a good guess for the initial weight - // of the proximal term and resolve. - if self.data.cur_weight.is_infinite() { - self.data.cur_weight = self.weighter.initial_weight(&self.data); - master.set_weight(self.data.cur_weight)?; - master.solve(self.data.cur_val)?; - return Ok(false); - } - - if self.terminator.terminate(&self.data) && !itdata.updated { - Self::show_info( - &self.start_time, - Step::Term, - &self.data, - self.cnt_descent, - self.cnt_null, - itdata.cnt_updates, - ); - info!("Termination criterion satisfied"); - return Ok(true); - } - - // Compress bundle - master.compress()?; - - // Compute new candidate. - let mut next_y = dvec![]; - itdata.nxt_d = Arc::new(master_res.nxt_d); - next_y.add(&self.data.cur_y, &itdata.nxt_d); - itdata.nxt_y = Arc::new(next_y); - - // Reset evaluation data. - itdata.nxt_ubs.clear(); - itdata.nxt_ubs.resize(self.problem.num_subproblems(), Real::infinity()); - itdata.cnt_remaining_ubs = self.problem.num_subproblems(); - itdata.nxt_cutvals.clear(); - itdata - .nxt_cutvals - .resize(self.problem.num_subproblems(), -Real::infinity()); - itdata.cnt_remaining_mins = self.problem.num_subproblems(); - - // Start evaluation of all subproblems at the new candidate. - let client_tx = self.client_tx.as_ref().ok_or(Error::NotInitialized)?; - for i in 0..self.problem.num_subproblems() { - self.problem - .evaluate(i, itdata.nxt_y.clone(), i, client_tx.clone()) - .map_err(Error::Evaluation)?; - } - Ok(false) - } - - fn update_problem( - problem: &mut P, - step: Step, - data: &mut SolverData, - itdata: &mut IterData, - master_proc: &mut MasterProcess, - ) -> Result> { - let (update_tx, update_rx) = channel(); - problem - .update( - &UpdateData { - cur_y: &data.cur_y, - nxt_y: &itdata.nxt_y, - step, - master_proc, - }, - itdata.cnt_iter, - update_tx, - ) - .map_err(Error::Update)?; - - let mut have_update = false; - for update in update_rx { - let update = update.map_err(Error::Update)?; - have_update = true; - match update { - Update::AddVariables { bounds, sgext, .. } => { - let mut newvars = Vec::with_capacity(bounds.len()); - for (lower, upper) in bounds { - if lower > upper { - return Err(Error::InvalidBounds { lower, upper }); - } - let value = if lower > 0.0 { - lower - } else if upper < 0.0 { - upper - } else { - 0.0 - }; - //self.bounds.push((lower, upper)); - newvars.push((None, lower - value, upper - value, value)); - } - if !newvars.is_empty() { - // modify moved variables - for (index, val) in newvars.iter().filter_map(|v| v.0.map(|i| (i, v.3))) { - data.cur_y[index] = val; - } - - // add new variables - data.cur_y.extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3)); - - master_proc.add_vars(newvars.iter().map(|v| (v.0, v.1, v.2)).collect(), sgext)?; - } - } - } - } - - Ok(have_update) - } - - /// Return the bound the function value must be below of to enforce a descent step. - /// - /// If the oracle guarantees that $f(\bar{y}) \le$ this bound, the - /// bundle method will perform a descent step. - /// - /// This value is $f(\hat{y}) + \varrho \cdot \Delta$ where - /// $\Delta = f(\hat{y}) - \hat{f}(\bar{y})$ is the expected - /// progress and $\varrho$ is the `acceptance_factor`. - fn get_descent_bound(acceptance_factor: Real, data: &SolverData) -> Real { - data.cur_val - acceptance_factor * (data.cur_val - data.nxt_mod) - } - - fn show_info( - start_time: &Instant, - step: Step, - data: &SolverData, - cnt_descent: usize, - cnt_null: usize, - cnt_updates: usize, - ) { - let time = start_time.elapsed(); - info!( - "{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4} {:4} {:4}{:1} {:9.4} {:9.4} \ - {:12.6e}({:12.6e}) {:12.6e}", - if step == Step::Term { "_endit" } else { "endit " }, - time.as_secs() / 3600, - (time.as_secs() / 60) % 60, - time.as_secs() % 60, - time.subsec_nanos() / 10_000_000, - cnt_descent, - cnt_descent + cnt_null, - cnt_updates, - if step == Step::Descent { "*" } else { " " }, - data.cur_weight, - data.expected_progress(), - data.nxt_mod, - data.nxt_val, - data.cur_val - ); - } - - /// Return the aggregated primal of the given subproblem. - pub fn aggregated_primal(&self, subproblem: usize) -> Result> { - Ok(self - .master_proc - .as_ref() - .ok_or(Error::NotSolved)? - .get_aggregated_primal(subproblem)?) - } -} ADDED src/problem.rs Index: src/problem.rs ================================================================== --- /dev/null +++ src/problem.rs @@ -0,0 +1,190 @@ +/* + * Copyright (c) 2019 Frank Fischer + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see + */ + +//! An asynchronous first-order oracle. + +use crate::{Aggregatable, DVector, Minorant, Real}; +use crossbeam::channel::Sender; +use std::sync::Arc; + +/// Evaluation result. +/// +/// The result of an evaluation is new information to be made +/// available to the solver and the master problem. There are +/// essentially two types of information: +/// +/// 1. The (exact) function value of a sub-function at some point. +/// 2. A minorant of some sub-function. +#[derive(Debug)] +pub enum EvalResult { + /// The objective value at some point. + ObjectiveValue { index: I, value: Real }, + /// A minorant with an associated primal. + Minorant { index: I, minorant: Minorant, primal: P }, +} + +/// Channel to send evaluation results to. +pub type ResultSender = Sender, E>>; + +/// Problem update information. +/// +/// The solver calls the `update` method of the problem regularly. +/// This method can modify the problem by adding (or moving) +/// variables. The possible updates are encoded in this type. +pub enum Update { + /// Add new variables with bounds. + /// + /// The initial value of each variable will be the feasible value + /// closest to 0. + AddVariables { + index: I, + bounds: Vec<(Real, Real)>, + sgext: Box>, + }, +} + +/// The subgradient extender is a callback used to update existing minorants +/// given their associated primal data. +pub type SubgradientExtender = dyn FnMut(usize, &P, &[usize]) -> Result + Send; + +/// This trait provides information available in the +/// [`FirstOrderProblem::update`] method. +pub trait UpdateState

{ + /// Whether the last step was a descent step. + fn was_descent(&self) -> bool; + + /// Whether the last step was a null step. + fn was_null(&self) -> bool { + !self.was_descent() + } + + /// The (old) current center of stability. + fn center(&self) -> Arc; + + /// The candidate point. + /// + /// After a descent step, i.e. if [`UpdateState::was_descent`] is `true`, + /// this is the new center of stability. + fn candidate(&self) -> Arc; + + /// The current aggregated primal information. + /// + /// Return the aggregated primal information for the given subproblem. + fn aggregated_primal(&self, i: usize) -> P; +} + +/// Channel to send problem updates to. +pub type UpdateSender = Sender, E>>; + +/// Trait for implementing a first-order problem description. +/// +/// All computations made by an implementation are supposed to +/// be asynchronous. Hence, the interface is slightly different +/// compared with [`crate::FirstOrderProblem`]. +pub trait FirstOrderProblem { + /// Error raised by this oracle. + type Err; + + /// The primal information associated with a minorant. + type Primal: Aggregatable + Send + 'static; + + /// Return the number of variables. + fn num_variables(&self) -> usize; + + /// Return the lower bounds on the variables. + /// + /// If no lower bounds a specified, $-\infty$ is assumed. + /// + /// The lower bounds must be less then or equal the upper bounds. + fn lower_bounds(&self) -> Option> { + None + } + + /** + * Return the upper bounds on the variables. + * + * If no lower bounds a specified, $+\infty$ is assumed. + * + * The upper bounds must be greater than or equal the upper bounds. + */ + fn upper_bounds(&self) -> Option> { + None + } + + /// Return the number of subproblems. + fn num_subproblems(&self) -> usize { + 1 + } + + /// Start background processes. + /// + /// This method is called right before the solver starts the solution process. + /// It can be used to setup any background tasks required for the evaluation + /// of the subfunctions. + /// + /// Remember that background processes should be cleanup when the problem + /// is deleted (e.g. by implementing the [`Drop`] trait). + /// + /// The default implementation does nothing. + fn start(&mut self) {} + + /// Stop background processes. + /// + /// This method is called right after the solver stops the solution process. + /// It can be used to stop any background tasks required for the evaluation + /// of the subfunctions. + /// + /// A correct implementation of should cleanup all processes from the [`Drop`] + /// thread. + /// + /// The default implementation does nothing. + fn stop(&mut self) {} + + /// Start the evaluation of the i^th subproblem at the given point. + /// + /// The results of the evaluation should be passed to the provided channel. + /// In order to work correctly, the results must contain (an upper bound on) + /// the objective value at $y$ as well as at least one subgradient centered + /// at $y$ eventually. + fn evaluate( + &mut self, + i: usize, + y: Arc, + index: I, + tx: ResultSender, + ) -> Result<(), Self::Err>; + + /// Called to update the problem. + /// + /// This method is called regularly by the solver. The problem should send problem update + /// information (e.g. adding new variables) to the provided channel. + /// + /// The updates might be generated asynchronously. + /// + /// The default implementation does nothing. + fn update( + &mut self, + _state: &U, + _index: I, + _tx: UpdateSender, + ) -> Result<(), Self::Err> + where + U: UpdateState, + { + Ok(()) + } +} Index: src/solver.rs ================================================================== --- src/solver.rs +++ src/solver.rs @@ -1,1068 +1,23 @@ -// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -//! The main bundle method solver. - -use crate::{Aggregatable, DVector, Real}; -use crate::{Evaluation, FirstOrderProblem, Update}; - -use crate::master::{self, MasterProblem}; -use crate::terminator::{StandardTerminatable, StandardTerminator, Terminator}; -use crate::weighter::{HKWeightable, HKWeighter, Weighter}; - -use log::{debug, info, warn}; - -use std::error::Error; -use std::f64::{INFINITY, NEG_INFINITY}; -use std::fmt; -use std::mem::swap; -use std::time::Instant; - -/// A solver error. -#[derive(Debug)] -pub enum SolverError { - /// An error occurred during oracle evaluation. - Evaluation(E), - /// An error occurred during oracle update. - Update(E), - /// An error has been raised by the master problem. - BuildMaster(Box), - /// An error has been raised by the master problem. - Master(MErr), - /// The oracle did not return a minorant. - NoMinorant, - /// The dimension of some data is wrong. - Dimension, - /// Some parameter has an invalid value. - Parameter(ParameterError), - /// The lower bound of a variable is larger than the upper bound. - InvalidBounds { lower: Real, upper: Real }, - /// The value of a variable is outside its bounds. - ViolatedBounds { lower: Real, upper: Real, value: Real }, - /// The variable index is out of bounds. - InvalidVariable { index: usize, nvars: usize }, - /// Iteration limit has been reached. - IterationLimit { limit: usize }, -} - -impl fmt::Display for SolverError -where - E: fmt::Display, - MErr: fmt::Display, -{ - fn fmt(&self, fmt: &mut fmt::Formatter) -> std::result::Result<(), fmt::Error> { - use self::SolverError::*; - match self { - Evaluation(err) => write!(fmt, "Oracle evaluation failed: {}", err), - Update(err) => write!(fmt, "Oracle update failed: {}", err), - BuildMaster(err) => write!(fmt, "Creation of master problem failed: {}", err), - Master(err) => write!(fmt, "Master problem failed: {}", err), - NoMinorant => write!(fmt, "The oracle did not return a minorant"), - Dimension => write!(fmt, "Dimension of lower bounds does not match number of variables"), - Parameter(msg) => write!(fmt, "Parameter error: {}", msg), - InvalidBounds { lower, upper } => write!(fmt, "Invalid bounds, lower:{}, upper:{}", lower, upper), - ViolatedBounds { lower, upper, value } => write!( - fmt, - "Violated bounds, lower:{}, upper:{}, value:{}", - lower, upper, value - ), - InvalidVariable { index, nvars } => { - write!(fmt, "Variable index out of bounds, got:{} must be < {}", index, nvars) - } - IterationLimit { limit } => write!(fmt, "The iteration limit of {} has been reached.", limit), - } - } -} - -impl Error for SolverError -where - E: Error + 'static, - MErr: Error + 'static, -{ - fn source(&self) -> Option<&(dyn Error + 'static)> { - match self { - SolverError::Evaluation(err) => Some(err), - SolverError::Update(err) => Some(err), - SolverError::Master(err) => Some(err), - SolverError::BuildMaster(err) => Some(err.as_ref()), - _ => None, - } - } -} - -impl From for SolverError { - fn from(err: MErr) -> SolverError { - SolverError::Master(err) - } -} - -/** - * The current state of the bundle method. - * - * Captures the current state of the bundle method during the run of - * the algorithm. This state is passed to certain callbacks like - * Terminator or Weighter so that they can compute their result - * depending on the state. - */ -pub struct BundleState<'a> { - /// Current center of stability. - pub cur_y: &'a DVector, - - /// Function value in current center. - pub cur_val: Real, - - /// Current candidate, point of last evaluation. - pub nxt_y: &'a DVector, - - /// Function value in candidate. - pub nxt_val: Real, - - /// Model value in candidate. - pub nxt_mod: Real, - - /// Cut value of new subgradient in current center. - pub new_cutval: Real, - - /// The current aggregated subgradient norm. - pub sgnorm: Real, - - /// The expected progress of the current model. - pub expected_progress: Real, - - /// Currently used weight of quadratic term. - pub weight: Real, - - /** - * The type of the current step. - * - * If the current step is Step::Term, the weighter should be reset. - */ - pub step: Step, -} - -macro_rules! current_state { - ($slf:ident, $step:expr) => { - BundleState { - cur_y: &$slf.cur_y, - cur_val: $slf.cur_val, - nxt_y: &$slf.nxt_y, - nxt_mod: $slf.nxt_mod, - nxt_val: $slf.nxt_val, - new_cutval: $slf.new_cutval, - sgnorm: $slf.sgnorm, - weight: $slf.master.weight(), - step: $step, - expected_progress: $slf.expected_progress, - } - }; -} - -impl<'a> HKWeightable for BundleState<'a> { - fn current_weight(&self) -> Real { - self.weight - } - - fn center(&self) -> &DVector { - self.cur_y - } - - fn center_value(&self) -> Real { - self.cur_val - } - - fn candidate_value(&self) -> Real { - self.nxt_val - } - - fn candidate_model(&self) -> Real { - self.nxt_mod - } - - fn new_cutvalue(&self) -> Real { - self.new_cutval - } - - fn sgnorm(&self) -> Real { - self.sgnorm - } -} - -impl<'a> StandardTerminatable for BundleState<'a> { - fn expected_progress(&self) -> Real { - self.expected_progress - } - - fn center_value(&self) -> Real { - self.cur_val - } -} - -/// An invalid value for some parameter has been passes. -#[derive(Debug)] -pub struct ParameterError(String); - -impl fmt::Display for ParameterError { - fn fmt(&self, fmt: &mut fmt::Formatter) -> std::result::Result<(), fmt::Error> { - write!(fmt, "{}", self.0) - } -} - -impl Error for ParameterError {} - -/// Parameters for tuning the solver. -#[derive(Clone, Debug)] -pub struct SolverParams { - /// Maximal individual bundle size. - pub max_bundle_size: usize, - - /** - * Factor for doing a descent step. - * - * If the proportion of actual decrease to predicted decrease is - * at least that high, a descent step will be done. - * - * Must be in (0,1). - */ - pub acceptance_factor: Real, - - /** - * Factor for doing a null step. - * - * Factor that guarantees a null step. This factor is used to - * compute a bound for the function oracle, that guarantees a null - * step. If the function is evaluated by some iterative method that ensures - * an objective value that is at least as large as this bound, the - * oracle can stop returning an appropriate $\varepsilon$-subgradient. - * - * Must be in (0, acceptance_factor). - */ - pub nullstep_factor: Real, -} - -impl SolverParams { - /// Verify that all parameters are valid. - fn check(&self) -> std::result::Result<(), ParameterError> { - if self.max_bundle_size < 2 { - Err(ParameterError(format!( - "max_bundle_size must be >= 2 (got: {})", - self.max_bundle_size - ))) - } else if self.acceptance_factor <= 0.0 || self.acceptance_factor >= 1.0 { - Err(ParameterError(format!( - "acceptance_factor must be in (0,1) (got: {})", - self.acceptance_factor - ))) - } else if self.nullstep_factor <= 0.0 || self.nullstep_factor > self.acceptance_factor { - Err(ParameterError(format!( - "nullstep_factor must be in (0,acceptance_factor] (got: {}, acceptance_factor:{})", - self.nullstep_factor, self.acceptance_factor - ))) - } else { - Ok(()) - } - } -} - -impl Default for SolverParams { - fn default() -> SolverParams { - SolverParams { - max_bundle_size: 50, - - nullstep_factor: 0.1, - acceptance_factor: 0.1, - } - } -} - -/// The step type that has been performed. -#[derive(Clone, Copy, PartialEq, Eq, Debug)] -pub enum Step { - /// A null step has been performed. - Null, - /// A descent step has been performed. - Descent, - /// No step but the algorithm has been terminated. - Term, -} - -/// Information about a minorant. -#[derive(Debug, Clone)] -struct MinorantInfo { - /// The minorant's index in the master problem - index: usize, - /// Current multiplier. - multiplier: Real, -} - -/// Information about the last iteration. -#[derive(Debug, Clone, Copy, PartialEq)] -pub enum IterationInfo { - NewMinorantTooHigh { new: Real, old: Real }, - UpperBoundNullStep, - ShallowCut, -} - -/// State information for the update callback. -pub struct UpdateState<'a, Pr: 'a> { - /// Current model minorants. - minorants: &'a [Vec], - /// The primals. - primals: &'a Vec>, - /// The last step type. - pub step: Step, - /// Iteration information. - pub iteration_info: &'a [IterationInfo], - /// The current candidate. If the step was a descent step, this is - /// the new center. - pub nxt_y: &'a DVector, - /// The center. IF the step was a descent step, this is the old - /// center. - pub cur_y: &'a DVector, -} - -impl<'a, Pr: 'a> UpdateState<'a, Pr> { - pub fn aggregated_primals(&self, subproblem: usize) -> Vec<(Real, &Pr)> { - self.minorants[subproblem] - .iter() - .map(|m| (m.multiplier, self.primals[m.index].as_ref().unwrap())) - .collect() - } - - /// Return the last primal for a given subproblem. - /// - /// This is the last primal generated by the oracle. - pub fn last_primal(&self, fidx: usize) -> Option<&Pr> { - self.minorants[fidx].last().and_then(|m| self.primals[m.index].as_ref()) - } -} - -/// The default builder. -pub type FullMasterBuilder = master::boxed::Builder; - -/** - * Implementation of a bundle method. - */ -pub struct Solver -where - P: FirstOrderProblem, - M: master::Builder, -{ - /// The first order problem description. - problem: P, - - /// The solver parameter. - pub params: SolverParams, - - /// Termination predicate. - pub terminator: T, - - /// Weighter heuristic. - pub weighter: W, - - /// Lower and upper bounds of all variables. - bounds: Vec<(Real, Real)>, - - /// Current center of stability. - cur_y: DVector, - - /// Function value in current point. - cur_val: Real, - - /// Model value in current point. - cur_mod: Real, - - /// Vector of subproblem function values in current point. - cur_vals: DVector, - - /// Vector of model values in current point. - cur_mods: DVector, - - /** - * Whether the data of the current center is valid. - * - * This variable is set to false of the problem data changes so - * the function is re-evaluated at the center. - */ - cur_valid: bool, - - /// Direction from current center to candidate. - nxt_d: DVector, - - /// Current candidate point. - nxt_y: DVector, - - /// (Upper bound on) function value in candidate. - nxt_val: Real, - - /// Model value in candidate. - nxt_mod: Real, - - /// DVector of subproblem function values in candidate. - nxt_vals: DVector, - - /// Vector of model values in candidate point. - nxt_mods: DVector, - - /// Cut value of new subgradient in current center. - new_cutval: Real, - - /// Norm of current aggregated subgradient. - sgnorm: Real, - - /// Expected progress. - expected_progress: Real, - - /// Number of descent steps. - cnt_descent: usize, - - /// Number of null steps. - cnt_null: usize, - - /** - * Time when the solution process started. - * - * This is actually the time of the last call to `Solver::init`. - */ - start_time: Instant, - - /// The master problem. - master: M::MasterProblem, - - /// The active minorant indices for each subproblem. - minorants: Vec>, - - /// The primals associated with each global minorant index. - primals: Vec>, - - /// Accumulated information about the last iteration. - iterinfos: Vec, -} - -pub type Result = std::result::Result< - T, - SolverError<

::Err, <::MasterProblem as MasterProblem>::Err>, ->; - -impl Solver -where - P: FirstOrderProblem, - P::Err: Into>, - T: for<'a> Terminator> + Default, - W: for<'a> Weighter> + Default, - M: master::Builder + Default, - M::MasterProblem: MasterProblem, -{ - /** - * Create a new solver for the given problem. - * - * Note that the solver owns the problem, so you cannot use the - * same problem description elsewhere as long as it is assigned to - * the solver. However, it is possible to get a reference to the - * internally stored problem using `Solver::problem()`. - */ - #[allow(clippy::type_complexity)] - pub fn new_params(problem: P, params: SolverParams) -> Result, P, M> { - Ok(Solver { - problem, - params, - terminator: T::default(), - weighter: W::default(), - bounds: vec![], - cur_y: dvec![], - cur_val: 0.0, - cur_mod: 0.0, - cur_vals: dvec![], - cur_mods: dvec![], - cur_valid: false, - nxt_d: dvec![], - nxt_y: dvec![], - nxt_val: 0.0, - nxt_mod: 0.0, - nxt_vals: dvec![], - nxt_mods: dvec![], - new_cutval: 0.0, - sgnorm: 0.0, - expected_progress: 0.0, - cnt_descent: 0, - cnt_null: 0, - start_time: Instant::now(), - master: M::default().build().map_err(|e| SolverError::BuildMaster(e.into()))?, - minorants: vec![], - primals: vec![], - iterinfos: vec![], - }) - } - - /// A new solver with default parameter. - #[allow(clippy::type_complexity)] - pub fn new(problem: P) -> Result, P, M> { - Solver::new_params(problem, SolverParams::default()) - } - - /** - * Set the first order problem description associated with this - * solver. - * - * Note that the solver owns the problem, so you cannot use the - * same problem description elsewhere as long as it is assigned to - * the solver. However, it is possible to get a reference to the - * internally stored problem using `Solver::problem()`. - */ - pub fn set_problem(&mut self, problem: P) { - self.problem = problem; - } - - /// Returns a reference to the solver's current problem. - pub fn problem(&self) -> &P { - &self.problem - } - - /// Initialize the solver. - pub fn init(&mut self) -> Result<(), P, M> { - self.params.check().map_err(SolverError::Parameter)?; - if self.cur_y.len() != self.problem.num_variables() { - self.cur_valid = false; - self.cur_y.init0(self.problem.num_variables()); - } - - let lb = self.problem.lower_bounds(); - let ub = self.problem.upper_bounds(); - self.bounds.clear(); - self.bounds.reserve(self.cur_y.len()); - for i in 0..self.cur_y.len() { - let lb_i = lb.as_ref().map(|x| x[i]).unwrap_or(NEG_INFINITY); - let ub_i = ub.as_ref().map(|x| x[i]).unwrap_or(INFINITY); - if lb_i > ub_i { - return Err(SolverError::InvalidBounds { - lower: lb_i, - upper: ub_i, - }); - } - if self.cur_y[i] < lb_i { - self.cur_valid = false; - self.cur_y[i] = lb_i; - } else if self.cur_y[i] > ub_i { - self.cur_valid = false; - self.cur_y[i] = ub_i; - } - self.bounds.push((lb_i, ub_i)); - } - - let m = self.problem.num_subproblems(); - self.cur_vals.init0(m); - self.cur_mods.init0(m); - self.nxt_vals.init0(m); - self.nxt_mods.init0(m); - - self.start_time = Instant::now(); - - Ok(()) - } - - /// Solve the problem with at most 10_000 iterations. - /// - /// Use `solve_with_limit` for an explicit iteration limit. - pub fn solve(&mut self) -> Result<(), P, M> { - const LIMIT: usize = 10_000; - self.solve_with_limit(LIMIT) - } - - /// Solve the problem with explicit iteration limit. - pub fn solve_with_limit(&mut self, iter_limit: usize) -> Result<(), P, M> { - // First initialize the internal data structures. - self.init()?; - - if self.solve_iter(iter_limit)? { - Ok(()) - } else { - Err(SolverError::IterationLimit { limit: iter_limit }) - } - } - - /// Solve the problem but stop after `niter` iterations. - /// - /// The function returns `Ok(true)` if the termination criterion - /// has been satisfied. Otherwise it returns `Ok(false)` or an - /// error code. - /// - /// If this function is called again, the solution process is - /// continued from the previous point. Because of this one must - /// call `init()` before the first call to this function. - pub fn solve_iter(&mut self, niter: usize) -> Result { - for _ in 0..niter { - let mut term = self.step()?; - let changed = self.update_problem(term)?; - // do not stop if the problem has been changed - if changed && term == Step::Term { - term = Step::Null - } - self.show_info(term); - if term == Step::Term { - return Ok(true); - } - } - Ok(false) - } - - /// Called to update the problem. - /// - /// Calling this function typically triggers the problem to - /// separate new constraints depending on the current solution. - fn update_problem(&mut self, term: Step) -> Result { - let updates = { - let state = UpdateState { - minorants: &self.minorants, - primals: &self.primals, - step: term, - iteration_info: &self.iterinfos, - // this is a dirty trick: when updating the center, we - // simply swapped the `cur_*` fields with the `nxt_*` - // fields - cur_y: if term == Step::Descent { - &self.nxt_y - } else { - &self.cur_y - }, - nxt_y: if term == Step::Descent { - &self.cur_y - } else { - &self.nxt_y - }, - }; - self.problem.update(&state).map_err(SolverError::Update)? - }; - - let mut newvars = Vec::with_capacity(updates.len()); - for u in updates { - match u { - Update::AddVariable { lower, upper } => { - if lower > upper { - return Err(SolverError::InvalidBounds { lower, upper }); - } - let value = if lower > 0.0 { - lower - } else if upper < 0.0 { - upper - } else { - 0.0 - }; - self.bounds.push((lower, upper)); - newvars.push((None, lower - value, upper - value, value)); - } - Update::AddVariableValue { lower, upper, value } => { - if lower > upper { - return Err(SolverError::InvalidBounds { lower, upper }); - } - if value < lower || value > upper { - return Err(SolverError::ViolatedBounds { lower, upper, value }); - } - self.bounds.push((lower, upper)); - newvars.push((None, lower - value, upper - value, value)); - } - Update::MoveVariable { index, value } => { - if index >= self.bounds.len() { - return Err(SolverError::InvalidVariable { - index, - nvars: self.bounds.len(), - }); - } - let (lower, upper) = self.bounds[index]; - if value < lower || value > upper { - return Err(SolverError::ViolatedBounds { lower, upper, value }); - } - newvars.push((Some(index), lower - value, upper - value, value)); - } - } - } - - if !newvars.is_empty() { - let problem = &mut self.problem; - let primals = &self.primals; - self.master.add_vars( - &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::>(), - &mut |fidx, minidx, vars| { - problem - .extend_subgradient(fidx, primals[minidx].as_ref().unwrap(), vars) - .map(DVector) - .map_err(|e| e.into()) - }, - )?; - // modify moved variables - for (index, val) in newvars.iter().filter_map(|v| v.0.map(|i| (i, v.3))) { - self.cur_y[index] = val; - self.nxt_y[index] = val; - self.nxt_d[index] = 0.0; - } - // add new variables - self.cur_y.extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3)); - self.nxt_y.extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3)); - self.nxt_d.resize(self.nxt_y.len(), 0.0); - Ok(true) - } else { - Ok(false) - } - } - - /// Return the current aggregated primal information for a subproblem. - /// - /// This function returns all currently used minorants $x_i$ along - /// with their coefficients $\alpha_i$. The aggregated primal can - /// be computed by combining the minorants $\bar{x} = - /// \sum_{i=1}\^m \alpha_i x_i$. - pub fn aggregated_primals(&self, subproblem: usize) -> P::Primal { - Aggregatable::combine( - self.minorants[subproblem] - .iter() - .map(|m| (m.multiplier, self.primals[m.index].as_ref().unwrap())), - ) - } - - fn show_info(&self, step: Step) { - let time = self.start_time.elapsed(); - info!( - "{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4} {:4} {:4}{:1} {:9.4} {:9.4} \ - {:12.6e}({:12.6e}) {:12.6e}", - if step == Step::Term { "_endit" } else { "endit " }, - time.as_secs() / 3600, - (time.as_secs() / 60) % 60, - time.as_secs() % 60, - time.subsec_nanos() / 10_000_000, - self.cnt_descent, - self.cnt_descent + self.cnt_null, - self.master.cnt_updates(), - if step == Step::Descent { "*" } else { " " }, - self.master.weight(), - self.expected_progress, - self.nxt_mod, - self.nxt_val, - self.cur_val - ); - } - - /// Return the current center of stability. - pub fn center(&self) -> &[Real] { - &self.cur_y - } - - /// Return the last candidate point. - pub fn candidate(&self) -> &[Real] { - &self.nxt_y - } - - /** - * Initializes the master problem. - * - * The oracle is evaluated once at the initial center and the - * master problem is initialized with the returned subgradient - * information. - */ - fn init_master(&mut self) -> Result<(), P, M> { - let m = self.problem.num_subproblems(); - - let lb = self.problem.lower_bounds().map(DVector); - let ub = self.problem.upper_bounds().map(DVector); - - if lb - .as_ref() - .map(|lb| lb.len() != self.problem.num_variables()) - .unwrap_or(false) - { - return Err(SolverError::Dimension); - } - if ub - .as_ref() - .map(|ub| ub.len() != self.problem.num_variables()) - .unwrap_or(false) - { - return Err(SolverError::Dimension); - } - - self.master.set_num_subproblems(m)?; - self.master.set_vars(self.problem.num_variables(), lb, ub)?; - - self.minorants = (0..m).map(|_| vec![]).collect(); - - self.cur_val = 0.0; - for i in 0..m { - let result = self - .problem - .evaluate(i, &self.cur_y, INFINITY, 0.0) - .map_err(SolverError::Evaluation)?; - self.cur_vals[i] = result.objective(); - self.cur_val += self.cur_vals[i]; - - let mut minorants = result.into_iter(); - if let Some((minorant, primal)) = minorants.next() { - self.cur_mods[i] = minorant.constant; - self.cur_mod += self.cur_mods[i]; - let minidx = self.master.add_minorant(i, minorant)?; - self.minorants[i].push(MinorantInfo { - index: minidx, - multiplier: 0.0, - }); - if minidx >= self.primals.len() { - self.primals.resize_with(minidx + 1, || None); - } - self.primals[minidx] = Some(primal); - } else { - return Err(SolverError::NoMinorant); - } - } - - self.cur_valid = true; - - // Solve the master problem once to compute the initial - // subgradient. - // - // We could compute that subgradient directly by - // adding up the initial minorants, but this would not include - // the eta terms. However, this is a heuristic anyway because - // we assume an initial weight of 1.0, which, in general, will - // *not* be the initial weight for the first iteration. - self.master.set_weight(1.0)?; - self.master.solve(self.cur_val)?; - self.sgnorm = self.master.get_dualoptnorm2().sqrt(); - - // Compute the real initial weight. - let state = current_state!(self, Step::Term); - let new_weight = self.weighter.initial_weight(&state); - self.master.set_weight(new_weight)?; - - debug!("Init master completed"); - - Ok(()) - } - - /// Solve the model (i.e. master problem) to compute the next candidate. - fn solve_model(&mut self) -> Result<(), P, M> { - self.master.solve(self.cur_val)?; - self.nxt_d = self.master.get_primopt(); - self.nxt_y.add(&self.cur_y, &self.nxt_d); - self.nxt_mod = self.master.get_primoptval(); - self.sgnorm = self.master.get_dualoptnorm2().sqrt(); - self.expected_progress = self.cur_val - self.nxt_mod; - - // update multiplier from master solution - for i in 0..self.problem.num_subproblems() { - for m in &mut self.minorants[i] { - m.multiplier = self.master.multiplier(m.index); - } - } - - debug!("Model result"); - debug!(" cur_val ={}", self.cur_val); - debug!(" nxt_mod ={}", self.nxt_mod); - debug!(" expected={}", self.expected_progress); - Ok(()) - } - - /// Reduce size of bundle. - fn compress_bundle(&mut self) -> Result<(), P, M> { - for i in 0..self.problem.num_subproblems() { - let n = self.master.num_minorants(i); - if n >= self.params.max_bundle_size { - // aggregate minorants with smallest coefficients - self.minorants[i].sort_by_key(|m| -((1e6 * m.multiplier) as isize)); - let aggr = self.minorants[i].split_off(self.params.max_bundle_size - 2); - let aggr_sum = aggr.iter().map(|m| m.multiplier).sum(); - let (aggr_mins, aggr_primals): (Vec<_>, Vec<_>) = aggr - .into_iter() - .map(|m| (m.index, self.primals[m.index].take().unwrap())) - .unzip(); - let (aggr_min, aggr_coeffs) = self.master.aggregate(i, &aggr_mins)?; - // append aggregated minorant - self.minorants[i].push(MinorantInfo { - index: aggr_min, - multiplier: aggr_sum, - }); - self.primals[aggr_min] = Some(Aggregatable::combine(aggr_coeffs.into_iter().zip(&aggr_primals))); - } - } - Ok(()) - } - - /// Perform a descent step. - fn descent_step(&mut self) -> Result<(), P, M> { - let new_weight = self.weighter.descent_weight(¤t_state!(self, Step::Descent)); - self.master.set_weight(new_weight)?; - self.cnt_descent += 1; - swap(&mut self.cur_y, &mut self.nxt_y); - swap(&mut self.cur_val, &mut self.nxt_val); - swap(&mut self.cur_mod, &mut self.nxt_mod); - swap(&mut self.cur_vals, &mut self.nxt_vals); - swap(&mut self.cur_mods, &mut self.nxt_mods); - self.master.move_center(1.0, &self.nxt_d); - debug!("Descent Step"); - debug!(" dir ={}", self.nxt_d); - debug!(" newy={}", self.cur_y); - Ok(()) - } - - /// Perform a null step. - fn null_step(&mut self) -> Result<(), P, M> { - let new_weight = self.weighter.null_weight(¤t_state!(self, Step::Null)); - self.master.set_weight(new_weight)?; - self.cnt_null += 1; - debug!("Null Step"); - Ok(()) - } - - /// Perform one bundle iteration. - #[allow(clippy::collapsible_if)] - pub fn step(&mut self) -> Result { - self.iterinfos.clear(); - - if !self.cur_valid { - // current point needs new evaluation - self.init_master()?; - } - - self.solve_model()?; - if self.terminator.terminate(¤t_state!(self, Step::Term)) { - return Ok(Step::Term); - } - - let m = self.problem.num_subproblems(); - let descent_bnd = self.get_descent_bound(); - let nullstep_bnd = if m == 1 { self.get_nullstep_bound() } else { INFINITY }; - let relprec = if m == 1 { self.get_relative_precision() } else { 0.0 }; - - self.compress_bundle()?; - - let mut nxt_lb = 0.0; - let mut nxt_ub = 0.0; - self.new_cutval = 0.0; - for fidx in 0..self.problem.num_subproblems() { - let result = self - .problem - .evaluate(fidx, &self.nxt_y, nullstep_bnd, relprec) - .map_err(SolverError::Evaluation)?; - let fun_ub = result.objective(); - - let mut minorants = result.into_iter(); - let mut nxt_minorant; - let nxt_primal; - match minorants.next() { - Some((m, p)) => { - nxt_minorant = m; - nxt_primal = p; - } - None => return Err(SolverError::NoMinorant), - } - let fun_lb = nxt_minorant.constant; - - nxt_lb += fun_lb; - nxt_ub += fun_ub; - self.nxt_vals[fidx] = fun_ub; - - // move center of minorant to cur_y - nxt_minorant.move_center(-1.0, &self.nxt_d); - self.new_cutval += nxt_minorant.constant; - let minidx = self.master.add_minorant(fidx, nxt_minorant)?; - self.minorants[fidx].push(MinorantInfo { - index: minidx, - multiplier: 0.0, - }); - if minidx >= self.primals.len() { - self.primals.resize_with(minidx + 1, || None); - } - self.primals[minidx] = Some(nxt_primal); - } - - if self.new_cutval > self.cur_val + 1e-3 { - warn!( - "New minorant has higher value in center new:{} old:{}", - self.new_cutval, self.cur_val - ); - self.cur_val = self.new_cutval; - self.iterinfos.push(IterationInfo::NewMinorantTooHigh { - new: self.new_cutval, - old: self.cur_val, - }); - } - - self.nxt_val = nxt_ub; - - // check for potential problems with relative precision of all kinds - if nxt_lb <= descent_bnd { - // lower bound gives descent step - if nxt_ub > descent_bnd { - // upper bound will produce null-step - if self.cur_val - nxt_lb > (self.cur_val - self.nxt_mod) * self.params.nullstep_factor.max(0.5) { - warn!("Relative precision of returned objective interval enforces null-step."); - self.iterinfos.push(IterationInfo::UpperBoundNullStep); - } - } - } else if self.cur_val - nxt_lb > 0.8 * (self.cur_val - self.nxt_mod) { - // TODO: double check with ConicBundle if this test makes sense. - // lower bound gives already a null step - // subgradient won't yield much improvement - warn!("Shallow cut (subgradient won't yield much improvement)"); - self.iterinfos.push(IterationInfo::ShallowCut); - } - - debug!("Step"); - debug!(" cur_val ={}", self.cur_val); - debug!(" nxt_mod ={}", self.nxt_mod); - debug!(" nxt_ub ={}", self.nxt_val); - debug!(" descent_bnd={}", descent_bnd); - - // do a descent step or null step - if nxt_ub <= descent_bnd { - self.descent_step()?; - Ok(Step::Descent) - } else { - self.null_step()?; - Ok(Step::Null) - } - } - - /** - * Return the bound on the function value that enforces a - * nullstep. - * - * If the oracle guarantees that $f(\bar{y}) \ge$ this bound, the - * bundle method will perform a nullstep. - * - * This value is $f(\hat{y}) + \varrho' \cdot \Delta$ where - * $\Delta = f(\hat{y}) - \hat{f}(\bar{y})$ is the expected - * progress and $\varrho'$ is the `nullstep_factor`. - */ - fn get_nullstep_bound(&self) -> Real { - self.cur_val - self.params.nullstep_factor * (self.cur_val - self.nxt_mod) - } - - /** - * Return the bound the function value must be below of to enforce a descent step. - * - * If the oracle guarantees that $f(\bar{y}) \le$ this bound, the - * bundle method will perform a descent step. - * - * This value is $f(\hat{y}) + \varrho \cdot \Delta$ where - * $\Delta = f(\hat{y}) - \hat{f}(\bar{y})$ is the expected - * progress and $\varrho$ is the `acceptance_factor`. - */ - fn get_descent_bound(&self) -> Real { - self.cur_val - self.params.acceptance_factor * (self.cur_val - self.nxt_mod) - } - - /** - * Return the required relative precision for the computation. - */ - fn get_relative_precision(&self) -> Real { - (0.1 * (self.cur_val - self.get_nullstep_bound()) / (self.cur_val.abs() + 1.0)).min(1e-3) - } -} +/* + * Copyright (c) 2019 Frank Fischer + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see + */ + +//! The basic solver implementation. + +pub mod sync; +pub use sync::{DefaultSolver, NoBundleSolver}; + +mod masterprocess; ADDED src/solver/masterprocess.rs Index: src/solver/masterprocess.rs ================================================================== --- /dev/null +++ src/solver/masterprocess.rs @@ -0,0 +1,274 @@ +/* + * Copyright (c) 2019 Frank Fischer + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see + */ + +//! Asynchronous process solving a master problem. + +use crossbeam::channel::{unbounded as channel, Receiver, Sender}; +use log::{debug, warn}; +use std::sync::Arc; +use threadpool::ThreadPool; + +use super::sync::Error; +use crate::master::primalmaster::PrimalMaster; +use crate::master::MasterProblem; +use crate::problem::{FirstOrderProblem, SubgradientExtender}; +use crate::{DVector, Minorant, Real}; + +/// Configuration information for setting up a master problem. +pub struct MasterConfig { + /// The number of subproblems. + pub num_subproblems: usize, + /// The number of variables. + pub num_vars: usize, + /// The lower bounds on the variables. + pub lower_bounds: Option, + /// The lower bounds on the variables. + pub upper_bounds: Option, +} + +/// A task for the master problem. +enum MasterTask +where + M: MasterProblem, +{ + /// Add new variables to the master problem. + AddVariables(Vec<(Option, Real, Real)>, Box>), + + /// Add a new minorant for a subfunction to the master problem. + AddMinorant(usize, Minorant, Pr), + + /// Move the center of the master problem in the given direction. + MoveCenter(Real, Arc), + + /// Start a new computation of the master problem. + Solve { center_value: Real }, + + /// Compress the bundle. + Compress, + + /// Set the weight parameter of the master problem. + SetWeight { weight: Real }, + + /// Return the current aggregated primal. + GetAggregatedPrimal { + subproblem: usize, + tx: Sender>, + }, +} + +/// The response send from a master process. +/// +/// The response contains the evaluation results of the latest +pub struct MasterResponse { + pub nxt_d: DVector, + pub nxt_mod: Real, + pub sgnorm: Real, + /// The number of internal iterations. + pub cnt_updates: usize, +} + +type ToMasterSender = Sender::Primal,

::Err, M>>; + +type ToMasterReceiver = Receiver::Primal,

::Err, M>>; + +type MasterSender = Sender>; + +pub type MasterReceiver = Receiver>; + +pub struct MasterProcess +where + P: FirstOrderProblem, + M: MasterProblem, +{ + /// The channel to transmit new tasks to the master problem. + tx: ToMasterSender, + + /// The channel to receive solutions from the master problem. + pub rx: MasterReceiver, + + phantom: std::marker::PhantomData, +} + +impl MasterProcess +where + P: FirstOrderProblem, + P::Primal: Send + 'static, + P::Err: Into> + 'static, + M: MasterProblem + Send + 'static, + M::MinorantIndex: std::hash::Hash, + M::Err: Send + 'static, +{ + pub fn start(master: M, master_config: MasterConfig, threadpool: &mut ThreadPool) -> Self { + // Create a pair of communication channels. + let (to_master_tx, to_master_rx) = channel(); + let (from_master_tx, from_master_rx) = channel(); + + // The the master process thread. + threadpool.execute(move || { + debug!("Master process started"); + let mut from_master_tx = from_master_tx; + if let Err(err) = Self::master_main(master, master_config, &mut from_master_tx, to_master_rx) { + #[allow(unused_must_use)] + { + from_master_tx.send(Err(err)); + } + } + debug!("Master proces stopped"); + }); + + MasterProcess { + tx: to_master_tx, + rx: from_master_rx, + phantom: std::marker::PhantomData, + } + } + + /// Add new variables to the master problem. + pub fn add_vars( + &mut self, + vars: Vec<(Option, Real, Real)>, + sgext: Box>, + ) -> Result<(), Error> + where + P::Err: 'static, + { + self.tx + .send(MasterTask::AddVariables(vars, sgext)) + .map_err(|err| Error::Process(err.into())) + } + + /// Add a new minorant to the master problem model. + /// + /// This adds the specified `minorant` with associated `primal` data to the + /// model of subproblem `i`. + pub fn add_minorant(&mut self, i: usize, minorant: Minorant, primal: P::Primal) -> Result<(), Error> { + self.tx + .send(MasterTask::AddMinorant(i, minorant, primal)) + .map_err(|err| Error::Process(err.into())) + } + + /// Move the center of the master problem. + /// + /// This moves the master problem's center in direction $\\alpha \\cdot d$. + pub fn move_center(&mut self, alpha: Real, d: Arc) -> Result<(), Error> { + self.tx + .send(MasterTask::MoveCenter(alpha, d)) + .map_err(|err| Error::Process(err.into())) + } + + /// Solve the master problem. + /// + /// The current function value in the center `center_value`. + /// Once the master problem is solved the process will send a + /// [`MasterResponse`] message to the `tx` channel. + pub fn solve(&mut self, center_value: Real) -> Result<(), Error> { + self.tx + .send(MasterTask::Solve { center_value }) + .map_err(|err| Error::Process(err.into())) + } + + /// Compresses the model. + pub fn compress(&mut self) -> Result<(), Error> { + self.tx + .send(MasterTask::Compress) + .map_err(|err| Error::Process(err.into())) + } + + /// Sets the new weight of the proximal term in the master problem. + pub fn set_weight(&mut self, weight: Real) -> Result<(), Error> { + self.tx + .send(MasterTask::SetWeight { weight }) + .map_err(|err| Error::Process(err.into())) + } + + /// Get the current aggregated primal for a certain subproblem. + pub fn get_aggregated_primal(&self, subproblem: usize) -> Result> { + let (tx, rx) = channel(); + self.tx + .send(MasterTask::GetAggregatedPrimal { subproblem, tx }) + .map_err(|err| Error::Process(err.into()))?; + rx.recv() + .map_err(|err| Error::Process(err.into()))? + .map_err(|err| Error::Master(err.into())) + } + + /// The main loop of the master process. + fn master_main( + master: M, + master_config: MasterConfig, + tx: &mut MasterSender, + rx: ToMasterReceiver, + ) -> Result<(), M::Err> { + let mut master = PrimalMaster::<_, P::Primal>::new(master); + + // Initialize the master problem. + master.set_num_subproblems(master_config.num_subproblems)?; + master.set_vars( + master_config.num_vars, + master_config.lower_bounds, + master_config.upper_bounds, + )?; + + // The main iteration: wait for new tasks. + for m in rx { + match m { + MasterTask::AddVariables(vars, sgext) => { + debug!("master: add {} variables to the subproblem", vars.len()); + master.add_vars(vars, sgext)?; + } + MasterTask::AddMinorant(i, m, primal) => { + debug!("master: add minorant to subproblem {}", i); + master.add_minorant(i, m, primal)?; + } + MasterTask::MoveCenter(alpha, d) => { + debug!("master: move center"); + master.move_center(alpha, &d); + } + MasterTask::Compress => { + debug!("Compress bundle"); + master.compress()?; + } + MasterTask::Solve { center_value } => { + debug!("master: solve with center_value {}", center_value); + master.solve(center_value)?; + let master_response = MasterResponse { + nxt_d: master.get_primopt(), + nxt_mod: master.get_primoptval(), + sgnorm: master.get_dualoptnorm2().sqrt(), + cnt_updates: master.cnt_updates(), + }; + if let Err(err) = tx.send(Ok(master_response)) { + warn!("Master process cancelled because of channel error: {}", err); + break; + } + } + MasterTask::SetWeight { weight } => { + debug!("master: set weight {}", weight); + master.set_weight(weight)?; + } + MasterTask::GetAggregatedPrimal { subproblem, tx } => { + debug!("master: get aggregated primal for {}", subproblem); + if tx.send(master.aggregated_primal(subproblem)).is_err() { + warn!("Sending of aggregated primal for {} failed", subproblem); + }; + } + }; + } + + Ok(()) + } +} ADDED src/solver/sync.rs Index: src/solver/sync.rs ================================================================== --- /dev/null +++ src/solver/sync.rs @@ -0,0 +1,898 @@ +/* + * Copyright (c) 2019 Frank Fischer + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 3 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see + */ + +//! An asynchronous parallel bundle solver. + +use crossbeam::channel::{select, unbounded as channel, Receiver, Sender}; +use log::{debug, info}; +use num_cpus; +use num_traits::Float; +use std::sync::Arc; +use std::time::Instant; +use threadpool::ThreadPool; + +use crate::{DVector, Real}; + +use super::masterprocess::{MasterConfig, MasterProcess, MasterResponse}; +use crate::master::{self, MasterProblem}; +use crate::problem::{EvalResult, FirstOrderProblem, Update, UpdateState}; +use crate::terminator::{StandardTerminatable, StandardTerminator, Terminator}; +use crate::weighter::{HKWeightable, HKWeighter, Weighter}; + +/// The default iteration limit. +pub const DEFAULT_ITERATION_LIMIT: usize = 10_000; + +/// The default solver. +pub type DefaultSolver

= Solver; + +/// The minimal bundle solver. +pub type NoBundleSolver

= Solver; + +/// Error raised by the parallel bundle [`Solver`]. +#[derive(Debug)] +pub enum Error { + /// An error raised when creating a new master problem solver. + BuildMaster(Box), + /// An error raised by the master problem process. + Master(Box), + /// The iteration limit has been reached. + IterationLimit { limit: usize }, + /// An error raised by a subproblem evaluation. + Evaluation(E), + /// An error raised subproblem update. + Update(E), + /// The dimension of some data is wrong. + Dimension(String), + /// Invalid bounds for a variable. + InvalidBounds { lower: Real, upper: Real }, + /// The value of a variable is outside its bounds. + ViolatedBounds { lower: Real, upper: Real, value: Real }, + /// The variable index is out of bounds. + InvalidVariable { index: usize, nvars: usize }, + /// An error occurred in a subprocess. + Process(Box), + /// A method requiring an initialized solver has been called. + NotInitialized, + /// The problem has not been solved yet. + NotSolved, +} + +impl std::fmt::Display for Error +where + E: std::fmt::Display, +{ + fn fmt(&self, fmt: &mut std::fmt::Formatter) -> std::result::Result<(), std::fmt::Error> { + use Error::*; + match self { + BuildMaster(err) => writeln!(fmt, "Cannot create master problem solver: {}", err), + Master(err) => writeln!(fmt, "Error in master problem: {}", err), + IterationLimit { limit } => writeln!(fmt, "The iteration limit has been reached: {}", limit), + Evaluation(err) => writeln!(fmt, "Error in subproblem evaluation: {}", err), + Update(err) => writeln!(fmt, "Error in subproblem update: {}", err), + Dimension(what) => writeln!(fmt, "Wrong dimension for {}", what), + InvalidBounds { lower, upper } => write!(fmt, "Invalid bounds, lower:{}, upper:{}", lower, upper), + ViolatedBounds { lower, upper, value } => write!( + fmt, + "Violated bounds, lower:{}, upper:{}, value:{}", + lower, upper, value + ), + InvalidVariable { index, nvars } => { + write!(fmt, "Variable index out of bounds, got:{} must be < {}", index, nvars) + } + Process(err) => writeln!(fmt, "Error in subprocess: {}", err), + NotInitialized => writeln!(fmt, "The solver must be initialized (called Solver::init()?)"), + NotSolved => writeln!(fmt, "The problem has not been solved yet"), + } + } +} + +impl std::error::Error for Error +where + E: std::error::Error + 'static, +{ + fn source(&self) -> Option<&(dyn std::error::Error + 'static)> { + use Error::*; + match self { + BuildMaster(err) => Some(err.as_ref()), + Master(err) => Some(err.as_ref()), + Evaluation(err) => Some(err), + Process(err) => Some(err.as_ref()), + _ => None, + } + } +} + +type ClientSender

= + Sender::Primal>,

::Err>>; + +type ClientReceiver

= + Receiver::Primal>,

::Err>>; + +/// Parameters for tuning the solver. +#[derive(Debug, Clone)] +pub struct Parameters { + /// The descent step acceptance factors, must be in (0,1). + /// + /// The default value is 0.1. + acceptance_factor: Real, +} + +impl Default for Parameters { + fn default() -> Self { + Parameters { acceptance_factor: 0.1 } + } +} + +impl Parameters { + /// Change the descent step acceptance factor. + /// + /// The default value is 0.1. + pub fn set_acceptance_factor(&mut self, acceptance_factor: Real) { + assert!( + acceptance_factor > 0.0 && acceptance_factor < 1.0, + "Descent step acceptance factors must be in (0,1), got: {}", + acceptance_factor + ); + self.acceptance_factor = acceptance_factor; + } +} + +/// The step type that has been performed. +#[derive(Clone, Copy, PartialEq, Eq, Debug)] +pub enum Step { + /// A null step has been performed. + Null, + /// A descent step has been performed. + Descent, + /// No step but the algorithm has been terminated. + Term, +} + +pub struct SolverData { + /// Current center of stability. + cur_y: DVector, + + /// Function value in the current point. + cur_val: Real, + + /// Function value at the current candidate. + nxt_val: Real, + + /// Model value at the current candidate. + nxt_mod: Real, + + /// The value of the new minorant in the current center. + new_cutval: Real, + + /// The current expected progress. + /// + /// This value is actually `cur_val - nxt_val`. We store it separately only + /// for debugging purposes because after a descent step `cur_val` will be + /// changed and we could not see the "old" expected progress anymore that + /// led to the descent step. + expected_progress: Real, + + /// Norm of current aggregated subgradient. + sgnorm: Real, + + /// The currently used master problem weight. + cur_weight: Real, +} + +impl SolverData { + /// Reset solver data to initial values. + /// + /// This means that almost everything is set to +infinity so that + /// a null-step is forced after the first evaluation. + fn init(&mut self, y: DVector) { + self.cur_y = y; + self.cur_val = Real::infinity(); + self.nxt_val = Real::infinity(); + self.nxt_mod = -Real::infinity(); + self.new_cutval = -Real::infinity(); + self.expected_progress = Real::infinity(); + self.sgnorm = Real::infinity(); + self.cur_weight = 1.0; + } +} + +impl StandardTerminatable for SolverData { + fn center_value(&self) -> Real { + self.cur_val + } + + fn expected_progress(&self) -> Real { + self.expected_progress + } +} + +impl HKWeightable for SolverData { + fn current_weight(&self) -> Real { + self.cur_weight + } + + fn center(&self) -> &DVector { + &self.cur_y + } + + fn center_value(&self) -> Real { + self.cur_val + } + + fn candidate_value(&self) -> Real { + self.nxt_val + } + + fn candidate_model(&self) -> Real { + self.nxt_mod + } + + fn new_cutvalue(&self) -> Real { + self.new_cutval + } + + fn sgnorm(&self) -> Real { + self.sgnorm + } +} + +/// Internal data used during the main iteration loop. +struct IterData { + /// Maximal number of iterations. + max_iter: usize, + cnt_iter: usize, + cnt_updates: usize, + nxt_ubs: Vec, + cnt_remaining_ubs: usize, + nxt_cutvals: Vec, + cnt_remaining_mins: usize, + nxt_d: Arc, + nxt_y: Arc, + /// True if the problem has been updated after the last evaluation. + updated: bool, +} + +impl IterData { + fn new(num_subproblems: usize, num_variables: usize, max_iter: usize) -> Self { + IterData { + max_iter, + cnt_iter: 0, + cnt_updates: 0, + nxt_ubs: vec![Real::infinity(); num_subproblems], + cnt_remaining_ubs: num_subproblems, + nxt_cutvals: vec![-Real::infinity(); num_subproblems], + cnt_remaining_mins: num_subproblems, + nxt_d: Arc::new(dvec![0.0; num_variables]), + nxt_y: Arc::new(dvec![]), + updated: true, + } + } +} + +/// Data providing access for updating the problem. +struct UpdateData<'a, P, M> +where + P: FirstOrderProblem, + M: MasterProblem, +{ + /// Type of step. + step: Step, + + /// Current center of stability. + cur_y: &'a DVector, + + /// Current candidate. + nxt_y: &'a Arc, + + /// The master process. + master_proc: &'a MasterProcess, +} + +impl<'a, P, M> UpdateState for UpdateData<'a, P, M> +where + P: FirstOrderProblem, + P::Err: Into> + 'static, + M: MasterProblem, + M::MinorantIndex: std::hash::Hash, +{ + fn was_descent(&self) -> bool { + self.step == Step::Descent + } + + fn center(&self) -> Arc { + Arc::new(self.cur_y.clone()) + } + + fn candidate(&self) -> Arc { + self.nxt_y.clone() + } + + fn aggregated_primal(&self, i: usize) -> P::Primal { + self.master_proc + .get_aggregated_primal(i) + .map_err(|_| "get_aggregated_primal".to_string()) + .expect("Cannot get aggregated primal from master process") + } +} + +/// Implementation of a parallel bundle method. +pub struct Solver +where + P: FirstOrderProblem, + M: master::Builder, +{ + /// Parameters for the solver. + pub params: Parameters, + + /// Termination predicate. + pub terminator: T, + + /// Weighter heuristic. + pub weighter: W, + + /// The threadpool of the solver. + pub threadpool: ThreadPool, + + /// The master problem builder. + pub master: M, + + /// The first order problem. + problem: P, + + /// The algorithm data. + data: SolverData, + + /// The master problem process. + master_proc: Option>, + + /// The channel to receive the evaluation results from subproblems. + client_tx: Option>, + + /// The channel to receive the evaluation results from subproblems. + client_rx: Option>, + + /// Number of descent steps. + cnt_descent: usize, + + /// Number of null steps. + cnt_null: usize, + + /// Number of function evaluation. + cnt_evals: usize, + + /// Time when the solution process started. + /// + /// This is actually the time of the last call to `Solver::init`. + start_time: Instant, +} + +impl Solver +where + P: FirstOrderProblem, + P::Err: Into> + 'static, + T: Terminator + Default, + W: Weighter + Default, + M: master::Builder, + M::MasterProblem: MasterProblem, + ::MinorantIndex: std::hash::Hash, +{ + /// Create a new parallel bundle solver. + pub fn new(problem: P) -> Self + where + M: Default, + { + Solver { + params: Parameters::default(), + terminator: Default::default(), + weighter: Default::default(), + problem, + data: SolverData { + cur_y: dvec![], + cur_val: 0.0, + nxt_val: 0.0, + nxt_mod: 0.0, + new_cutval: 0.0, + expected_progress: 0.0, + sgnorm: 0.0, + cur_weight: 1.0, + }, + + threadpool: ThreadPool::with_name("Parallel bundle solver".to_string(), num_cpus::get()), + master: M::default(), + master_proc: None, + client_tx: None, + client_rx: None, + + cnt_descent: 0, + cnt_null: 0, + cnt_evals: 0, + + start_time: Instant::now(), + } + } + + /// Create a new parallel bundle solver. + pub fn with_master(problem: P, master: M) -> Self { + Solver { + params: Parameters::default(), + terminator: Default::default(), + weighter: Default::default(), + problem, + data: SolverData { + cur_y: dvec![], + cur_val: 0.0, + nxt_val: 0.0, + nxt_mod: 0.0, + new_cutval: 0.0, + expected_progress: 0.0, + sgnorm: 0.0, + cur_weight: 1.0, + }, + + threadpool: ThreadPool::with_name("Parallel bundle solver".to_string(), num_cpus::get()), + master, + master_proc: None, + client_tx: None, + client_rx: None, + + cnt_descent: 0, + cnt_null: 0, + cnt_evals: 0, + + start_time: Instant::now(), + } + } + + /// Return the underlying threadpool. + /// + /// In order to use the same threadpool for concurrent processes, + /// just clone the returned `ThreadPool`. + pub fn threadpool(&self) -> &ThreadPool { + &self.threadpool + } + + /// Set the threadpool. + /// + /// This function allows to use a specific threadpool for all processes + /// spawned by the solver. Note that this does not involve any threads + /// used by the problem because the solver is not responsible for executing + /// the evaluation process of the subproblems. However, the problem might + /// use the same threadpool as the solver. + pub fn set_threadpool(&mut self, threadpool: ThreadPool) { + self.threadpool = threadpool; + } + + /// Return the current problem associated with the solver. + pub fn problem(&self) -> &P { + &self.problem + } + + /// Initialize the solver. + /// + /// This will reset the internal data structures so that a new fresh + /// solution process can be started. + /// + /// It will also setup all worker processes. + /// + /// This function is automatically called by [`Solver::solve`]. + pub fn init(&mut self) -> Result<(), Error> { + debug!("Initialize solver"); + + let n = self.problem.num_variables(); + let m = self.problem.num_subproblems(); + + self.data.init(dvec![0.0; n]); + self.cnt_descent = 0; + self.cnt_null = 0; + self.cnt_evals = 0; + + let (tx, rx) = channel(); + self.client_tx = Some(tx); + self.client_rx = Some(rx); + + let master_config = MasterConfig { + num_subproblems: m, + num_vars: n, + lower_bounds: self.problem.lower_bounds().map(DVector), + upper_bounds: self.problem.upper_bounds().map(DVector), + }; + + if master_config + .lower_bounds + .as_ref() + .map(|lb| lb.len() != n) + .unwrap_or(false) + { + return Err(Error::Dimension("lower bounds".to_string())); + } + if master_config + .upper_bounds + .as_ref() + .map(|ub| ub.len() != n) + .unwrap_or(false) + { + return Err(Error::Dimension("upper bounds".to_string())); + } + + debug!("Start master process"); + self.master_proc = Some(MasterProcess::start( + self.master.build().map_err(|err| Error::BuildMaster(err.into()))?, + master_config, + &mut self.threadpool, + )); + + debug!("Initial problem evaluation"); + // We need an initial evaluation of all oracles for the first center. + let y = Arc::new(self.data.cur_y.clone()); + for i in 0..m { + self.problem + .evaluate(i, y.clone(), i, self.client_tx.clone().unwrap()) + .map_err(Error::Evaluation)?; + } + + debug!("Initialization complete"); + + self.start_time = Instant::now(); + + Ok(()) + } + + /// Solve the problem with the default maximal iteration limit [`DEFAULT_ITERATION_LIMIT`]. + pub fn solve(&mut self) -> Result<(), Error> { + self.solve_with_limit(DEFAULT_ITERATION_LIMIT) + } + + /// Solve the problem with a maximal iteration limit. + pub fn solve_with_limit(&mut self, limit: usize) -> Result<(), Error> { + // First initialize the internal data structures. + self.init()?; + + if self.solve_iter(limit)? { + Ok(()) + } else { + Err(Error::IterationLimit { limit }) + } + } + + /// Solve the problem but stop after at most `niter` iterations. + /// + /// The function returns `Ok(true)` if the termination criterion + /// has been satisfied. Otherwise it returns `Ok(false)` or an + /// error code. + /// + /// If this function is called again, the solution process is + /// continued from the previous point. Because of this one *must* + /// call `init()` before the first call to this function. + pub fn solve_iter(&mut self, niter: usize) -> Result> { + debug!("Start solving up to {} iterations", niter); + + let mut itdata = IterData::new(self.problem.num_subproblems(), self.problem.num_variables(), niter); + + loop { + select! { + recv(self.client_rx.as_ref().ok_or(Error::NotInitialized)?) -> msg => { + let msg = msg + .map_err(|err| Error::Process(err.into()))? + .map_err(Error::Evaluation)?; + if self.handle_client_response(msg, &mut itdata)? { + return Ok(false); + } + }, + recv(self.master_proc.as_ref().ok_or(Error::NotInitialized)?.rx) -> msg => { + debug!("Receive master response"); + // Receive result (new candidate) from the master + let master_res = msg + .map_err(|err| Error::Process(err.into()))? + .map_err(|err| Error::Master(err.into()))?; + + if self.handle_master_response(master_res, &mut itdata)? { + return Ok(true); + } + }, + } + } + } + + /// Handle a response from a subproblem evaluation. + /// + /// The function returns `Ok(true)` if the final iteration count has been reached. + fn handle_client_response( + &mut self, + msg: EvalResult::Primal>, + itdata: &mut IterData, + ) -> Result> { + let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; + match msg { + EvalResult::ObjectiveValue { index, value } => { + debug!("Receive objective from subproblem {}: {}", index, value); + if itdata.nxt_ubs[index].is_infinite() { + itdata.cnt_remaining_ubs -= 1; + } + itdata.nxt_ubs[index] = itdata.nxt_ubs[index].min(value); + } + EvalResult::Minorant { + index, + mut minorant, + primal, + } => { + debug!("Receive minorant from subproblem {}", index); + if itdata.nxt_cutvals[index].is_infinite() { + itdata.cnt_remaining_mins -= 1; + } + // move center of minorant to cur_y + minorant.move_center(-1.0, &itdata.nxt_d); + itdata.nxt_cutvals[index] = itdata.nxt_cutvals[index].max(minorant.constant); + // add minorant to master problem + master.add_minorant(index, minorant, primal)?; + } + } + + if itdata.cnt_remaining_ubs > 0 || itdata.cnt_remaining_mins > 0 { + // Haven't received data from all subproblems, yet. + return Ok(false); + } + + // All subproblems have been evaluated, do a step. + let nxt_ub = itdata.nxt_ubs.iter().sum::(); + let descent_bnd = Self::get_descent_bound(self.params.acceptance_factor, &self.data); + + self.data.nxt_val = nxt_ub; + self.data.new_cutval = itdata.nxt_cutvals.iter().sum::(); + + debug!("Step"); + debug!(" cur_val ={}", self.data.cur_val); + debug!(" nxt_mod ={}", self.data.nxt_mod); + debug!(" nxt_ub ={}", nxt_ub); + debug!(" descent_bnd={}", descent_bnd); + + itdata.updated = false; + let step; + if self.data.cur_val.is_infinite() { + // This is the first evaluation. We effectively get + // the function value at the current center but + // we do not have a model estimate yet. Hence, we do not know + // a good guess for the weight. + step = Step::Descent; + self.data.cur_val = nxt_ub; + self.data.cur_weight = Real::infinity(); + master.set_weight(1.0)?; + + itdata.updated = true; + + debug!("First Step"); + debug!(" cur_val={}", self.data.cur_val); + debug!(" cur_y={}", self.data.cur_y); + } else if nxt_ub <= descent_bnd { + step = Step::Descent; + self.cnt_descent += 1; + + // Note that we must update the weight *before* we + // change the internal data, so the old information + // that caused the descent step is still available. + self.data.cur_weight = self.weighter.descent_weight(&self.data); + self.data.cur_y = itdata.nxt_y.as_ref().clone(); + self.data.cur_val = nxt_ub; + + master.move_center(1.0, itdata.nxt_d.clone())?; + master.set_weight(self.data.cur_weight)?; + + debug!("Descent Step"); + debug!(" dir ={}", itdata.nxt_d); + debug!(" newy={}", self.data.cur_y); + } else { + step = Step::Null; + self.cnt_null += 1; + self.data.cur_weight = self.weighter.null_weight(&self.data); + master.set_weight(self.data.cur_weight)?; + } + + Self::show_info( + &self.start_time, + step, + &self.data, + self.cnt_descent, + self.cnt_null, + itdata.cnt_updates, + ); + itdata.cnt_iter += 1; + + // Update problem. + if Self::update_problem(&mut self.problem, step, &mut self.data, itdata, master)? { + itdata.updated = true; + } + + // Compute the new candidate. The main loop will wait for the result of + // this solution process of the master problem. + master.solve(self.data.cur_val)?; + + Ok(itdata.cnt_iter >= itdata.max_iter) + } + + fn handle_master_response( + &mut self, + master_res: MasterResponse, + itdata: &mut IterData, + ) -> Result> { + let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; + + self.data.nxt_mod = master_res.nxt_mod; + self.data.sgnorm = master_res.sgnorm; + self.data.expected_progress = self.data.cur_val - self.data.nxt_mod; + itdata.cnt_updates = master_res.cnt_updates; + + // If this is the very first solution of the model, + // we use its result as to make a good guess for the initial weight + // of the proximal term and resolve. + if self.data.cur_weight.is_infinite() { + self.data.cur_weight = self.weighter.initial_weight(&self.data); + master.set_weight(self.data.cur_weight)?; + master.solve(self.data.cur_val)?; + return Ok(false); + } + + if self.terminator.terminate(&self.data) && !itdata.updated { + Self::show_info( + &self.start_time, + Step::Term, + &self.data, + self.cnt_descent, + self.cnt_null, + itdata.cnt_updates, + ); + info!("Termination criterion satisfied"); + return Ok(true); + } + + // Compress bundle + master.compress()?; + + // Compute new candidate. + let mut next_y = dvec![]; + itdata.nxt_d = Arc::new(master_res.nxt_d); + next_y.add(&self.data.cur_y, &itdata.nxt_d); + itdata.nxt_y = Arc::new(next_y); + + // Reset evaluation data. + itdata.nxt_ubs.clear(); + itdata.nxt_ubs.resize(self.problem.num_subproblems(), Real::infinity()); + itdata.cnt_remaining_ubs = self.problem.num_subproblems(); + itdata.nxt_cutvals.clear(); + itdata + .nxt_cutvals + .resize(self.problem.num_subproblems(), -Real::infinity()); + itdata.cnt_remaining_mins = self.problem.num_subproblems(); + + // Start evaluation of all subproblems at the new candidate. + let client_tx = self.client_tx.as_ref().ok_or(Error::NotInitialized)?; + for i in 0..self.problem.num_subproblems() { + self.problem + .evaluate(i, itdata.nxt_y.clone(), i, client_tx.clone()) + .map_err(Error::Evaluation)?; + } + Ok(false) + } + + fn update_problem( + problem: &mut P, + step: Step, + data: &mut SolverData, + itdata: &mut IterData, + master_proc: &mut MasterProcess, + ) -> Result> { + let (update_tx, update_rx) = channel(); + problem + .update( + &UpdateData { + cur_y: &data.cur_y, + nxt_y: &itdata.nxt_y, + step, + master_proc, + }, + itdata.cnt_iter, + update_tx, + ) + .map_err(Error::Update)?; + + let mut have_update = false; + for update in update_rx { + let update = update.map_err(Error::Update)?; + have_update = true; + match update { + Update::AddVariables { bounds, sgext, .. } => { + let mut newvars = Vec::with_capacity(bounds.len()); + for (lower, upper) in bounds { + if lower > upper { + return Err(Error::InvalidBounds { lower, upper }); + } + let value = if lower > 0.0 { + lower + } else if upper < 0.0 { + upper + } else { + 0.0 + }; + //self.bounds.push((lower, upper)); + newvars.push((None, lower - value, upper - value, value)); + } + if !newvars.is_empty() { + // modify moved variables + for (index, val) in newvars.iter().filter_map(|v| v.0.map(|i| (i, v.3))) { + data.cur_y[index] = val; + } + + // add new variables + data.cur_y.extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3)); + + master_proc.add_vars(newvars.iter().map(|v| (v.0, v.1, v.2)).collect(), sgext)?; + } + } + } + } + + Ok(have_update) + } + + /// Return the bound the function value must be below of to enforce a descent step. + /// + /// If the oracle guarantees that $f(\bar{y}) \le$ this bound, the + /// bundle method will perform a descent step. + /// + /// This value is $f(\hat{y}) + \varrho \cdot \Delta$ where + /// $\Delta = f(\hat{y}) - \hat{f}(\bar{y})$ is the expected + /// progress and $\varrho$ is the `acceptance_factor`. + fn get_descent_bound(acceptance_factor: Real, data: &SolverData) -> Real { + data.cur_val - acceptance_factor * (data.cur_val - data.nxt_mod) + } + + fn show_info( + start_time: &Instant, + step: Step, + data: &SolverData, + cnt_descent: usize, + cnt_null: usize, + cnt_updates: usize, + ) { + let time = start_time.elapsed(); + info!( + "{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4} {:4} {:4}{:1} {:9.4} {:9.4} \ + {:12.6e}({:12.6e}) {:12.6e}", + if step == Step::Term { "_endit" } else { "endit " }, + time.as_secs() / 3600, + (time.as_secs() / 60) % 60, + time.as_secs() % 60, + time.subsec_nanos() / 10_000_000, + cnt_descent, + cnt_descent + cnt_null, + cnt_updates, + if step == Step::Descent { "*" } else { " " }, + data.cur_weight, + data.expected_progress(), + data.nxt_mod, + data.nxt_val, + data.cur_val + ); + } + + /// Return the aggregated primal of the given subproblem. + pub fn aggregated_primal(&self, subproblem: usize) -> Result> { + Ok(self + .master_proc + .as_ref() + .ok_or(Error::NotSolved)? + .get_aggregated_primal(subproblem)?) + } +} DELETED src/vector.rs Index: src/vector.rs ================================================================== --- src/vector.rs +++ /dev/null @@ -1,287 +0,0 @@ -// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer -// -// This program is free software: you can redistribute it and/or -// modify it under the terms of the GNU General Public License as -// published by the Free Software Foundation, either version 3 of the -// License, or (at your option) any later version. -// -// This program is distributed in the hope that it will be useful, but -// WITHOUT ANY WARRANTY; without even the implied warranty of -// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU -// General Public License for more details. -// -// You should have received a copy of the GNU General Public License -// along with this program. If not, see -// - -//! Finite-dimensional sparse and dense vectors. - -use crate::{Aggregatable, Real}; -use std::fmt; -use std::ops::{Deref, DerefMut}; -// use std::cmp::min; -use std::borrow::Borrow; -use std::iter::FromIterator; -use std::vec::IntoIter; - -#[cfg(feature = "blas")] -use {openblas_src as _, rs_blas as blas, std::os::raw::c_int}; - -/// Type of dense vectors. -#[derive(Debug, Clone, PartialEq, Default)] -pub struct DVector(pub Vec); - -impl Deref for DVector { - type Target = Vec; - - fn deref(&self) -> &Vec { - &self.0 - } -} - -impl DerefMut for DVector { - fn deref_mut(&mut self) -> &mut Vec { - &mut self.0 - } -} - -impl fmt::Display for DVector { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - write!(f, "(")?; - for (i, x) in self.iter().enumerate() { - if i > 0 { - write!(f, ", ")?; - } - write!(f, "{}", x)? - } - write!(f, ")")?; - Ok(()) - } -} - -impl FromIterator for DVector { - fn from_iter>(iter: I) -> Self { - DVector(Vec::from_iter(iter)) - } -} - -impl IntoIterator for DVector { - type Item = Real; - type IntoIter = IntoIter; - - fn into_iter(self) -> IntoIter { - self.0.into_iter() - } -} - -/// Type of dense or vectors. -#[derive(Debug, Clone)] -pub enum Vector { - /// A vector with dense storage. - Dense(DVector), - - /** - * A vector with sparse storage. - * - * For each non-zero element this vector stores an index and the - * value of the element in addition to the size of the vector. - */ - Sparse { size: usize, elems: Vec<(usize, Real)> }, -} - -impl fmt::Display for Vector { - fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - match *self { - Vector::Dense(ref v) => write!(f, "{}", v), - Vector::Sparse { size, ref elems } => { - let mut it = elems.iter(); - write!(f, "{}:(", size)?; - if let Some(&(i, x)) = it.next() { - write!(f, "{}:{}", i, x)?; - for &(i, x) in it { - write!(f, ", {}:{}", i, x)?; - } - } - write!(f, ")") - } - } - } -} - -impl DVector { - /// Set all elements to 0. - pub fn init0(&mut self, size: usize) { - self.clear(); - self.extend((0..size).map(|_| 0.0)); - } - - /// Set self = factor * y. - pub fn scal(&mut self, factor: Real, y: &DVector) { - self.clear(); - self.extend(y.iter().map(|y| factor * y)); - } - - /// Return factor * self. - pub fn scaled(&self, factor: Real) -> DVector { - let mut x = DVector::default(); - x.scal(factor, self); - x - } - - /// Return the inner product with another vector. - pub fn dot(&self, other: &DVector) -> Real { - assert_eq!(self.len(), other.len()); - self.dot_begin(other) - } - - /// Return the inner product with another vector. - /// - /// The inner product is computed on the smaller of the two - /// dimensions. All other elements are assumed to be zero. - pub fn dot_begin(&self, other: &DVector) -> Real { - #[cfg(feature = "blas")] - unsafe { - blas::ddot(self.len().min(other.len()) as c_int, &self, 1, &other, 1) - } - #[cfg(not(feature = "blas"))] - { - self.iter().zip(other.iter()).map(|(x, y)| x * y).sum::() - } - } - - /// Add two vectors and store result in this vector. - pub fn add(&mut self, x: &DVector, y: &DVector) { - assert_eq!(x.len(), y.len()); - self.clear(); - self.extend(x.iter().zip(y.iter()).map(|(a, b)| a + b)); - } - - /// Add two vectors and store result in this vector. - pub fn add_scaled(&mut self, alpha: Real, y: &DVector) { - assert_eq!(self.len(), y.len()); - #[cfg(feature = "blas")] - unsafe { - blas::daxpy(self.len() as c_int, alpha, &y, 1, &mut self[..], 1) - } - #[cfg(not(feature = "blas"))] - { - for (x, y) in self.iter_mut().zip(y.iter()) { - *x += alpha * y; - } - } - } - - /// Add two vectors and store result in this vector. - /// - /// In contrast to `add_scaled`, the two vectors might have - /// different sizes. The size of the resulting vector is the - /// larger of the two vector sizes and the remaining entries of - /// the smaller vector are assumed to be 0.0. - pub fn add_scaled_begin(&mut self, alpha: Real, y: &DVector) { - #[cfg(feature = "blas")] - unsafe { - let n = self.len(); - blas::daxpy(n.min(y.len()) as c_int, alpha, &y, 1, &mut self[..], 1); - } - #[cfg(not(feature = "blas"))] - { - for (x, y) in self.iter_mut().zip(y.iter()) { - *x += alpha * y; - } - } - let n = self.len(); - if n < y.len() { - self.extend(y[n..].iter().map(|y| alpha * y)); - } - } - - /// Return the 2-norm of this vector. - pub fn norm2(&self) -> Real { - #[cfg(feature = "blas")] - unsafe { - blas::dnrm2(self.len() as c_int, &self, 1) - } - #[cfg(not(feature = "blas"))] - { - self.iter().map(|x| x * x).sum::().sqrt() - } - } -} - -impl Aggregatable for DVector { - fn new_scaled(alpha: Real, other: A) -> Self - where - A: Borrow, - { - DVector::scaled(&other.borrow(), alpha) - } - - fn add_scaled(&mut self, alpha: Real, other: A) - where - A: Borrow, - { - DVector::add_scaled(self, alpha, &other.borrow()) - } -} - -impl Vector { - /** - * Return a sparse vector with the given non-zeros. - */ - pub fn new_sparse(n: usize, indices: &[usize], values: &[Real]) -> Vector { - assert_eq!(indices.len(), values.len()); - - if indices.is_empty() { - Vector::Sparse { size: n, elems: vec![] } - } else { - let mut ordered: Vec<_> = (0..n).collect(); - ordered.sort_by_key(|&i| indices[i]); - assert!(*indices.last().unwrap() < n); - let mut elems = Vec::with_capacity(indices.len()); - let mut last_idx = n; - for i in ordered { - let val = unsafe { *values.get_unchecked(i) }; - if val != 0.0 { - let idx = unsafe { *indices.get_unchecked(i) }; - if idx != last_idx { - elems.push((idx, val)); - last_idx = idx; - } else { - elems.last_mut().unwrap().1 += val; - if elems.last_mut().unwrap().1 == 0.0 { - elems.pop(); - last_idx = n; - } - } - } - } - Vector::Sparse { size: n, elems } - } - } - - /** - * Convert vector to a dense vector. - * - * This function always returns a copy of the vector. - */ - pub fn to_dense(&self) -> DVector { - match *self { - Vector::Dense(ref x) => x.clone(), - Vector::Sparse { size: n, elems: ref xs } => { - let mut v = vec![0.0; n]; - for &(i, x) in xs { - unsafe { *v.get_unchecked_mut(i) = x }; - } - DVector(v) - } - } - } -} - -#[test] -fn test_add_scaled_begin() { - let mut x = dvec![1.0; 5]; - let y = dvec![2.0; 7]; - x.add_scaled_begin(3.0, &y); - assert_eq!(x, dvec![7.0, 7.0, 7.0, 7.0, 7.0, 6.0, 6.0]); -}