/*
* 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/>
*/
use {Real, DVector, Minorant};
use master::master::{MasterProblem, Error, Result};
use master::UnconstrainedMasterProblem;
use std::f64::{INFINITY, NEG_INFINITY};
/**
* Turn unconstrained master problem into box-constrained one.
*
* This master problem adds box constraints to an unconstrainted
* master problem implementation. The box constraints are enforced by
* an additional outer optimization loop.
*/
pub struct BoxedMasterProblem<M : UnconstrainedMasterProblem> {
lb : DVector,
ub : DVector,
eta : DVector,
/// Primal optimal solution.
primopt : DVector,
/// Primal optimal solution value.
primoptval : Real,
/// Square of norm of dual optimal solution.
dualoptnorm2: Real,
/// Model precision.
model_eps: Real,
need_new_candidate: bool,
/// Maximal number of updates of box multipliers.
max_updates : usize,
/// Current number of updates.
cnt_updates : usize,
/// The unconstrained master problem solver.
master: M,
}
impl<M : UnconstrainedMasterProblem> BoxedMasterProblem<M> {
pub fn new() -> Result<BoxedMasterProblem<M>> {
Ok(BoxedMasterProblem{
lb : dvec![],
ub : dvec![],
eta : dvec![],
primopt : dvec![],
primoptval: 0.0,
dualoptnorm2: 0.0,
model_eps: 0.6,
max_updates : 100,
cnt_updates : 0,
need_new_candidate : true,
master : match M::new() {
Ok(m) => m,
Err(e) => return Err(Error::Solver(Box::new(e))),
},
})
}
pub fn set_max_updates(&mut self, max_updates: usize) {
assert!(max_updates > 0);
self.max_updates = max_updates;
}
/**
* Update box multipliers $\eta$.
*
* This function solves the dual problem for fixed aggregated
* minorant w.r.t. $\eta$. When called, the variable `self.primopt`
* should contain the primal solution (i.e. the new candidate point
* d) without the influence of $\eta$.
*/
fn update_box_multipliers(&mut self) -> bool {
let mut updated_eta = false;
let weight = self.master.weight();
if self.eta.len() != self.lb.len() {
self.eta.resize(self.lb.len(), 0.0);
}
for i in 0..self.lb.len() {
let mut b = self.lb[i];
let x = self.primopt[i];
if x >= b {
b = self.ub[i];
if x <= b {
self.eta[i] = 0.0;
continue;
}
}
self.primopt[i] = b;
let neweta = (b - x) * weight;
if neweta != self.eta[i] { updated_eta = true; }
self.eta[i] = neweta;
}
debug!("Eta update");
debug!(" primopt={}", self.primopt);
debug!(" eta ={}", self.eta);
return updated_eta;
}
/*
* Compute the new candidate point.
*
* This consists of two steps:
*
* 1. the new point is computed as $-\tfrac{1}{u}\bar{g}$, where $\bar{g}$
* is the aggregated minorant
* 2. the multipliers $\eta$ are updated
*
* In other words, this function computes the new candidate
* defined by a fixed $\bar{g}$ while choosing the best possible
* $\eta$.
*/
fn compute_candidate(&mut self) {
self.need_new_candidate = false;
if self.master.dualopt().len() == self.lb.len() {
self.primopt.scal(-1.0/self.master.weight(), self.master.dualopt())
} else {
self.primopt.init0(self.lb.len());
}
self.update_box_multipliers();
}
/// Compute $\langle b, \eta \rangle$ with $b$ the bounds of eta.
fn eta_cutval(&self) -> Real {
let mut val = 0.0;
for i in 0..self.lb.len() {
if self.eta[i] >= 0.0 {
if self.lb[i] > NEG_INFINITY {
val += self.lb[i] * self.eta[i];
}
} else {
if self.ub[i] < INFINITY {
val += self.ub[i] * self.eta[i];
}
}
}
return val;
}
/**
* Return $\\|G \alpha - \eta\\|_2\^2$.
*
* This is the norm-square of the dual optimal solution including
* the current box-multipliers $\eta$.
*/
fn get_norm_subg2(&self) -> Real {
let dualopt = self.master.dualopt();
let mut norm2 = 0.0;
for i in 0..self.eta.len() {
let x = dualopt[i] - self.eta[i];
norm2 += x * x;
}
return norm2;
}
}
impl<M : UnconstrainedMasterProblem> MasterProblem for BoxedMasterProblem<M> {
fn set_vars(&mut self, n: usize, lb : Option<DVector>, ub: Option<DVector>) {
assert!(lb.as_ref().map(|x| x.len()).unwrap_or(n) == n);
assert!(ub.as_ref().map(|x| x.len()).unwrap_or(n) == n);
self.lb = lb.unwrap_or_else(|| dvec![NEG_INFINITY; n]);
self.ub = ub.unwrap_or_else(|| dvec![INFINITY; n]);
}
fn num_minorants(&self) -> usize { self.master.num_minorants() }
fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<()> {
self.master.add_minorant(fidx, minorant).map_err(|err| Error::Solver(Box::new(err)))
}
fn weight(&self) -> Real {
self.master.weight()
}
fn set_weight(&mut self, weight: Real) {
self.master.set_weight(weight);
}
fn solve(&mut self, center_value: Real) -> Result<()> {
debug!("Solve Master");
debug!(" lb ={}", self.lb);
debug!(" ub ={}", self.ub);
if self.need_new_candidate {
self.compute_candidate();
}
let mut cnt_updates = 0;
let mut old_augval = NEG_INFINITY;
loop {
cnt_updates += 1;
self.cnt_updates += 1;
// TODO: relprec is fixed
if let Err(err) = self.master.solve(&self.eta, center_value, old_augval, 1e-3) {
return Err(Error::Solver(Box::new(err)));
}
// compute the primal solution without the influence of eta
self.primopt.scal(-1.0 / self.master.weight(), self.master.dualopt());
// solve w.r.t. eta
let updated_eta = self.update_box_multipliers();
// compute value of the linearized model
self.dualoptnorm2 = self.get_norm_subg2();
let linval = self.master.dualopt().dot(&self.primopt) + self.master.dualopt_cutval();
assert!((self.eta.dot(&self.primopt) - self.eta_cutval()).abs() < 1e-6); // verify complementarity
let augval = linval + 0.5 * self.dualoptnorm2 / self.master.weight();
let mut cutval = linval;
if updated_eta {
cutval = self.master.eval_model(&self.primopt)
}
let curval = center_value;
let model_prec = (cutval - linval) / (curval - linval).max(1e-16);
debug!("Eta Test");
debug!(" dualnorm2={}", self.dualoptnorm2);
debug!(" linval={}", linval);
debug!(" modval={}", self.master.eval_model(&self.primopt));
debug!(" augval={}", augval);
debug!(" cutval={}", cutval);
debug!(" model_prec={}", model_prec);
debug!(" old_augval={}", old_augval);
debug!(" center_value={}", center_value);
debug!(" model_eps={}", self.model_eps);
debug!(" cut-lin={} < eps*(cur-lin)={}", cutval - linval, self.model_eps * (curval - linval));
debug!(" cnt_update={} max_updates={}", cnt_updates, self.max_updates);
self.primoptval = linval;
if augval < old_augval + 1e-10 ||
cutval - linval < self.model_eps * (curval - linval) ||
cnt_updates >= self.max_updates
{
break;
}
old_augval = old_augval.max(augval);
}
debug!("Model");
debug!(" cnt_update={}", cnt_updates);
debug!(" primopt={}", self.primopt);
debug!(" dualopt={}", self.master.dualopt());
debug!(" etaopt={}", self.eta);
debug!(" primoptval={}", self.primoptval);
Ok(())
}
fn aggregate(&mut self, fidx: usize, mins: &[usize]) {
self.master.aggregate(fidx, mins)
}
fn get_primopt(&self) -> DVector { self.primopt.clone() }
fn get_primoptval(&self) -> Real { self.primoptval }
fn get_dualoptnorm2(&self) -> Real { self.dualoptnorm2 }
fn move_center(&mut self, alpha: Real, d: &DVector) {
self.need_new_candidate = true;
self.master.move_center(alpha, d);
for i in 0..self.primopt.len() {
self.lb[i] -= alpha * d[i];
self.ub[i] -= alpha * d[i];
}
}
fn set_max_updates(&mut self, max_updates: usize) {
BoxedMasterProblem::set_max_updates(self, max_updates);
}
fn cnt_updates(&self) -> usize {
self.cnt_updates
}
}