/*
* Copyright (c) 2016 Frank Fischer <frank-fischer@shadow-soft.de>
*
* 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 <http://www.gnu.org/licenses/>
*/
//! The main bundle method solver.
use {Real, DVector};
use {FirstOrderProblem, Evaluation, HKWeighter};
use master::{self, MasterProblem, BoxedMasterProblem, MinimalMaster, CplexMaster};
use std::result;
use std::error;
use std::mem::swap;
use std::f64::{INFINITY, NEG_INFINITY};
use std::time::Instant;
quick_error! {
/// A solver error.
#[derive(Debug)]
pub enum Error {
/// An error occurred during evaluation of the oracle.
Eval(err: Box<error::Error>) {
cause(&**err)
description(err.description())
display("Evaluation error: {}", err)
}
/// Error solving the master problem.
Master(err: Box<error::Error>) {
cause(&**err)
description(err.description())
display("Master problem error: {}", err)
from(err: master::Error) -> (Box::new(err))
}
/// The oracle did not return a minorant.
NoMinorant {
description("No minorant")
display("The oracle did not return a minorant")
}
/// The dimension of some data is wrong.
Dimension(msg: &'static str) {
description("Dimension error")
display("Dimension error: {}", msg)
}
/// Some parameter has an invalid value.
Parameter(msg: String) {
description("Parameter error")
display("Parameter error: {}", msg)
}
}
}
/// Result type for solvers.
pub type Result<T> = result::Result<T, Error>;
/**
* 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,
}
impl<'a> BundleState<'a> {
}
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,
}
};
}
/**
* Termination predicate.
*
* Given the current state of the bundle method, this function returns
* whether the solution process should be stopped.
*/
pub trait Terminator {
/// Return true if the method should stop.
fn terminate(&mut self, state : &BundleState, params: &SolverParams) -> bool;
}
/**
* Terminates if expected progress is small enough.
*/
pub struct StandardTerminator {
pub termination_precision: Real,
}
impl Terminator for StandardTerminator {
#[allow(unused_variables)]
fn terminate(&mut self, state : &BundleState, params: &SolverParams) -> bool {
assert!(self.termination_precision >= 0.0);
state.expected_progress <= self.termination_precision * (state.cur_val.abs() + 1.0)
}
}
/**
* Bundle weight controller.
*
* Given the current state of the bundle method, this function determines the
* weight factor of the quadratic term for the next iteration.
*/
pub trait Weighter {
/// Return the new weight of the quadratic term.
fn weight(&mut self, state : &BundleState, params: &SolverParams) -> Real;
}
/// 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,
/// Minimal allowed bundle weight. Must be > 0 and < max_weight.
pub min_weight: Real,
/// Maximal allowed bundle weight. Must be > min_weight,
pub max_weight: Real,
/**
* Maximal number of updates of box multipliers.
*
* This is the maximal number of iterations for updating the box
* multipliers when solving the master problem with box
* constraints. This is a technical parameter that should probably
* never be changed. If you experience an unexpectedly high number
* of inner iterations, consider removing/fixing the corresponding
* variables.
*/
pub max_updates : usize,
}
impl SolverParams {
/// Verify that all parameters are valid.
fn check(&self) -> Result<()> {
if self.max_bundle_size < 2 {
Err(Error::Parameter(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(Error::Parameter(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(Error::Parameter(format!("nullstep_factor must be in (0,acceptance_factor] (got: {}, acceptance_factor:{})",
self.nullstep_factor,
self.acceptance_factor)))
} else if self.min_weight <= 0.0 {
Err(Error::Parameter(format!("min_weight must be in > 0 (got: {})", self.min_weight)))
} else if self.max_weight < self.min_weight {
Err(Error::Parameter(format!("max_weight must be in >= min_weight (got: {}, min_weight: {})",
self.max_weight, self.min_weight)))
} else if self.max_updates == 0 {
Err(Error::Parameter(format!("max_updates must be in > 0 (got: {})", self.max_updates)))
} else {
Ok(())
}
}
}
impl Default for SolverParams {
fn default() -> SolverParams {
SolverParams {
max_bundle_size: 50,
nullstep_factor: 0.1,
acceptance_factor: 0.1,
min_weight: 0.01,
max_weight: 1000.0,
max_updates: 50,
}
}
}
/// 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, Copy)]
struct MinorantInfo {
/// The minorant's index in the master problem
index: usize,
/// Current multiplier.
multiplier: usize,
}
/**
* Implementation of a bundle method.
*/
pub struct Solver<P>
where P : for <'a> FirstOrderProblem<'a>
{
/// The first order problem description.
problem : P,
/// The solver parameter.
pub params : SolverParams,
/// Termination predicate.
pub terminator: Box<Terminator>,
/// Weighter heuristic.
pub weighter: Box<Weighter>,
/// 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: Box<MasterProblem<MinorantIndex=usize>>,
/// The active minorant indices for each subproblem.
minorants: Vec<Vec<MinorantInfo>>,
}
impl<P> Solver<P>
where P : for <'a> FirstOrderProblem<'a>
{
/**
* 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()`.
*/
pub fn new_params(problem : P, params : SolverParams) -> Result<Solver<P>> {
Ok(Solver{
problem: problem,
params: params,
terminator: Box::new(StandardTerminator { termination_precision: 1e-3 }),
weighter: Box::new(HKWeighter::new()),
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: match BoxedMasterProblem::<MinimalMaster>::new() {
Ok(master) => Box::new(master),
Err(err) => return Err(Error::Master(Box::new(err))),
},
minorants: vec![],
})
}
/// A new solver with default parameter.
pub fn new(problem : P) -> Result<Solver<P>> {
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<()> {
try!(self.params.check());
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().map(|x| x.to_dense());
let ub = self.problem.upper_bounds().map(|x| x.to_dense());
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 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;
}
}
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.
pub fn solve(&mut self) -> Result<()> {
try!(self.init());
for _ in 0..100000 {
let term = try!(self.step());
self.show_info(term);
if term == Step::Term {
break;
}
}
Ok(())
}
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() / 10000000,
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);
}
/**
* 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<()> {
let m = self.problem.num_subproblems();
self.master = if m == 1 && self.params.max_bundle_size == 2 {
debug!("Use minimal master problem");
Box::new(BoxedMasterProblem::<MinimalMaster>::new().unwrap())
} else {
debug!("Use CPLEX master problem");
Box::new(BoxedMasterProblem::<CplexMaster>::new().unwrap())
};
let lb = self.problem.lower_bounds().map(|v| v.to_dense());
let ub = self.problem.upper_bounds().map(|v| v.to_dense());
if let Some(ref x) = lb {
if x.len() != self.problem.num_variables() {
return Err(Error::Dimension("Dimension of lower bounds does not match number of variables"));
}
}
try!(self.master.set_num_subproblems(m));
self.master.set_vars(self.problem.num_variables(), lb, ub);
self.master.set_max_updates(self.params.max_updates);
self.minorants = vec![vec![]; m];
self.cur_val = 0.0;
for i in 0..m {
let result = match self.problem.evaluate(i, &self.cur_y, INFINITY, 0.0) {
Ok(r) => r,
Err(err) => return Err(Error::Eval(Box::new(err))),
};
self.cur_vals[i] = result.objective();
self.cur_val += self.cur_vals[i];
let mut minorants = result.into_iter();
if let Some(minorant) = minorants.next() {
self.cur_mods[i] = minorant.constant;
self.cur_mod += self.cur_mods[i];
self.minorants[i].push(MinorantInfo {
index: try!(self.master.add_minorant(i, minorant)),
multiplier: 0,
});
} else {
return Err(Error::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);
try!(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.weight(&state, &self.params);
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<()> {
try!(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;
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<()> {
for i in 0..self.problem.num_subproblems() {
let n = self.master.num_minorants(i);
if n >= self.params.max_bundle_size {
for m in self.minorants[i].iter_mut() {
m.multiplier = (1e6 * self.master.multiplier(m.index)) as usize;
}
self.minorants[i].sort_by_key(|m| -(m.multiplier as isize));
let remove = self.minorants[i].split_off(self.params.max_bundle_size-2);
self.minorants[i].push(MinorantInfo{
index: try!(self.master.aggregate(i, &remove.iter().map(|m| m.index).collect::<Vec<_>>())),
multiplier: 0,
});
}
}
Ok(())
}
/// Perform a descent step.
fn descent_step(&mut self) {
let new_weight = self.weighter.weight(¤t_state!(self, Step::Descent), &self.params);
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);
}
/// Perform a null step.
fn null_step(&mut self) {
let new_weight = self.weighter.weight(¤t_state!(self, Step::Null), &self.params);
self.master.set_weight(new_weight);
self.cnt_null += 1;
debug!("Null Step");
}
/// Perform one bundle iteration.
pub fn step(&mut self) -> Result<Step> {
if !self.cur_valid {
// current point needs new evaluation
try!(self.init_master());
}
try!(self.solve_model());
if self.terminator.terminate(¤t_state!(self, Step::Term), &self.params) {
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 };
try!(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 = match self.problem.evaluate(fidx, &self.nxt_y, nullstep_bnd, relprec) {
Ok(r) => r,
Err(err) => return Err(Error::Eval(Box::new(err))),
};
let fun_ub = result.objective();
let mut minorants = result.into_iter();
let mut nxt_minorant = match minorants.next() {
Some(m) => m,
None => return Err(Error::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;
self.minorants[fidx].push(MinorantInfo{
index: try!(self.master.add_minorant(fidx, nxt_minorant)),
multiplier: 0,
});
}
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.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.")
}
}
} else {
// lower bound gives already a null step
if self.cur_val - nxt_lb > 0.8 * (self.cur_val - self.nxt_mod) {
// subgradient won't yield much improvement
warn!("Shallow cut (subgradient won't yield much improvement)");
}
}
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();
return Ok(Step::Descent);
} else {
self.null_step();
return 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)
}
}