// Copyright (c) 2016, 2017, 2018, 2019 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/>
//
//! 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::{self, NEG_INFINITY};
use std::os::raw::{c_char, c_int};
use std::ptr;
#[derive(Debug)]
pub enum CplexMasterError {
Cplex(cpx::CplexError),
SubgradientExtension(Box<dyn std::error::Error + Send + Sync>),
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<cpx::CplexError> for CplexMasterError {
fn from(err: cpx::CplexError) -> Self {
CplexMasterError::Cplex(err)
}
}
pub type Result<T> = std::result::Result<T, CplexMasterError>;
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<usize>,
/// List of minorant indices to be updated.
updateinds: Vec<usize>,
/// Mapping minorant to index.
min2index: Vec<Vec<usize>>,
/// Mapping index to minorant.
index2min: Vec<(usize, usize)>,
/// The quadratic term.
qterm: Vec<DVector>,
/// The weight of the quadratic term.
weight: Real,
/// The minorants for each subproblem in the model.
minorants: Vec<Vec<Minorant>>,
/// Optimal multipliers for each subproblem in the model.
opt_mults: Vec<DVector>,
/// 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<CplexMaster> {
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![],
min2index: vec![],
index2min: vec![],
qterm: vec![],
weight: 1.0,
minorants: 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.min2index = vec![vec![]; n];
self.index2min.clear();
self.freeinds.clear();
self.minorants = vec![vec![]; n];
self.opt_mults = vec![dvec![]; n];
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<F>(&mut self, f: F) -> Result<()>
where
F: FnMut(Self::MinorantIndex, &mut Iterator<Item = (Self::MinorantIndex, Real)>),
{
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::<Vec<_>>();
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.min2index[i][j])
.collect::<Vec<_>>();
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<usize> {
debug!("Add minorant");
debug!(" fidx={} index={}: {}", fidx, self.minorants[fidx].len(), minorant);
let min_idx = self.minorants[fidx].len();
self.minorants[fidx].push(minorant);
self.opt_mults[fidx].push(0.0);
self.force_update = true;
if let Some(idx) = self.freeinds.pop() {
self.min2index[fidx].push(idx);
self.index2min[idx] = (fidx, min_idx);
self.updateinds.push(idx);
Ok(idx)
} else {
let idx = self.index2min.len();
self.min2index[fidx].push(idx);
self.index2min.push((fidx, min_idx));
self.updateinds.push(idx);
Ok(idx)
}
}
fn add_vars(
&mut self,
nnew: usize,
changed: &[usize],
extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
) -> Result<()> {
debug_assert!(!self.minorants[0].is_empty());
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() {
if !mins.is_empty() {
for (i, m) in mins.iter_mut().enumerate() {
let new_subg = extend_subgradient(fidx, self.min2index[fidx][i], &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
if self.force_update {
return Ok(());
}
for (fidx_i, mins_i) in self.minorants.iter().enumerate() {
for (i, m_i) in mins_i.iter().enumerate() {
let idx_i = self.min2index[fidx_i][i];
for (fidx_j, mins_j) in self.minorants.iter().enumerate() {
for (j, m_j) in mins_j.iter().enumerate() {
let idx_j = self.min2index[fidx_j][j];
if idx_i <= idx_j {
let x: Real = (nnewvars..noldvars).map(|k| m_i.linear[k] * m_j.linear[k]).sum();
self.qterm[idx_i][idx_j] += x;
self.qterm[idx_j][idx_i] = self.qterm[idx_i][idx_j];
}
}
}
}
}
// WORST CASE: DO THIS
// 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 };
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 mins in &self.minorants {
for m in mins {
inds.push(c.len() as c_int);
c.push(-m.constant * self.weight - m.linear.dot(eta));
}
}
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 fidx in 0..self.minorants.len() {
for i in 0..self.minorants[fidx].len() {
self.opt_mults[fidx][i] = sol[idx];
mults.push(sol[idx]);
mins.push(&self.minorants[fidx][i]);
idx += 1;
}
}
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 (fidx, idx) = self.index2min[min];
self.opt_mults[fidx][idx]
}
fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box<Iterator<Item = (Self::MinorantIndex, Real)> + 'a> {
Box::new(
self.opt_mults[fidx]
.iter()
.cloned()
.enumerate()
.map(move |(i, alpha)| (self.min2index[fidx][i], alpha)),
)
}
fn eval_model(&self, y: &DVector) -> Real {
let mut result = 0.0;
for mins in &self.minorants {
let mut this_val = NEG_INFINITY;
for m in 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 {
sum_coeffs += self.opt_mults[fidx][self.index2min[i].1];
}
let aggr_coeffs = if sum_coeffs != 0.0 {
mins.iter()
.map(|&i| self.opt_mults[fidx][self.index2min[i].1] / sum_coeffs)
.collect::<DVector>()
} else {
dvec![0.0; mins.len()]
};
// compute aggregated diagonal term
let mut aggr_diag = 0.0;
for (idx_i, &i) in mins.iter().enumerate() {
for (idx_j, &j) in mins.iter().enumerate() {
aggr_diag += aggr_coeffs[idx_i] * aggr_coeffs[idx_j] * self.qterm[i][j];
}
}
// compute aggregated off-diagonal terms
let mut aggr_qterm = dvec![0.0; self.qterm.len()];
for fidx_i in 0..self.minorants.len() {
for idx_i in 0..self.minorants[fidx_i].len() {
let i = self.min2index[fidx_i][idx_i];
for (idx_j, &j) in mins.iter().enumerate() {
aggr_qterm[i] += aggr_coeffs[idx_j] * self.qterm[i][j];
}
}
}
// aggregate the minorants
let aggr = {
let mut aggr_mins = Vec::with_capacity(mins.len());
for &i in mins {
let (min_fidx, min_idx) = self.index2min[i];
debug_assert_eq!(min_fidx, fidx);
let m = self.minorants[fidx].swap_remove(min_idx);
let idx = self.min2index[fidx].swap_remove(min_idx);
self.opt_mults[fidx].swap_remove(min_idx);
self.freeinds.push(idx);
debug_assert_eq!(idx, i);
aggr_mins.push(m);
// update index2min table for moved minorant
if min_idx < self.minorants[fidx].len() {
self.index2min[self.min2index[fidx][min_idx]].1 = min_idx;
}
}
Aggregatable::combine(aggr_coeffs.iter().cloned().zip(&aggr_mins))
};
// save aggregated minorant
let aggr_idx = self.freeinds.pop().unwrap();
self.minorants[fidx].push(aggr);
self.opt_mults[fidx].push(sum_coeffs);
self.min2index[fidx].push(aggr_idx);
self.index2min[aggr_idx] = (fidx, self.minorants[fidx].len() - 1);
// update qterm
for fidx_i in 0..self.minorants.len() {
for idx_i in 0..self.minorants[fidx_i].len() {
let i = self.min2index[fidx_i][idx_i];
self.qterm[i][aggr_idx] = aggr_qterm[i];
self.qterm[aggr_idx][i] = aggr_qterm[i];
}
}
self.qterm[aggr_idx][aggr_idx] = aggr_diag;
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);
}
}
}
}
impl CplexMaster {
fn init_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
});
if self.force_update {
self.updateinds.clear();
for inds in &self.min2index {
self.updateinds.extend(inds.iter());
}
}
let nfun = self.minorants.len();
let mut nvars;
// add constraints
{
let sense: Vec<c_char> = vec!['E' as c_char; nfun];
let rhs = dvec![1.0; nfun];
let mut rmatbeg = Vec::with_capacity(nfun);
let mut rmatind = Vec::with_capacity(self.index2min.len());
let mut rmatval = Vec::with_capacity(self.index2min.len());
nvars = 0;
for i in 0..nfun {
rmatbeg.push(nvars as c_int);
for _ in 0..self.minorants[i].len() {
rmatind.push(nvars as c_int);
rmatval.push(1.0);
nvars += 1;
}
}
trycpx!(cpx::addrows(
self.env,
self.lp,
nvars as c_int,
nfun as c_int,
nvars as c_int,
rhs.as_ptr(),
sense.as_ptr(),
rmatbeg.as_ptr(),
rmatind.as_ptr(),
rmatval.as_ptr(),
ptr::null(),
ptr::null()
));
}
// build quadratic term
{
self.qterm.resize(self.index2min.len(), dvec![]);
for i in 0..self.qterm.len() {
self.qterm[i].resize(self.index2min.len(), 0.0);
}
// the global indices for each minorant in order
let mut activeinds = vec![];
for (fidx, mins_i) in self.minorants.iter().enumerate() {
for (i, m_i) in mins_i.iter().enumerate() {
let idx_i = self.min2index[fidx][i];
activeinds.push(idx_i);
for &idx_j in &self.updateinds {
let (fidx_j, j) = self.index2min[idx_j];
let x = m_i.linear.dot(&self.minorants[fidx_j][j].linear);
self.qterm[idx_i][idx_j] = x;
self.qterm[idx_j][idx_i] = x;
}
}
}
if cfg!(debug_assertions) {
for &idx_i in &activeinds {
for &idx_j in &activeinds {
let (fidx_i, i) = self.index2min[idx_i];
let (fidx_j, j) = self.index2min[idx_j];
let m_i = &self.minorants[fidx_i][i];
let m_j = &self.minorants[fidx_j][j];
debug_assert!((self.qterm[idx_i][idx_j] - m_i.linear.dot(&m_j.linear)).abs() < 1e-6);
}
}
}
// main diagonal plus small identity to ensure Q being semi-definite
let mut maxq = 0.0;
for &i in &activeinds {
maxq = f64::max(maxq, self.qterm[i][i])
}
maxq *= 1e-8;
// update coefficients
for (i, &idx_i) in activeinds.iter().enumerate() {
for (j, &idx_j) in activeinds.iter().enumerate() {
if i != j {
trycpx!(cpx::chgqpcoef(
self.env,
self.lp,
i as c_int,
j as c_int,
self.qterm[idx_i][idx_j]
));
} else {
trycpx!(cpx::chgqpcoef(
self.env,
self.lp,
i as c_int,
j as c_int,
self.qterm[idx_i][idx_j] + maxq
));
}
}
}
}
self.updateinds.clear();
self.force_update = false;
Ok(())
}
}
pub struct Builder {
/// The maximal bundle size used in the master problem.
pub max_bundle_size: usize,
}
impl Default for Builder {
fn default() -> Self {
Builder { max_bundle_size: 50 }
}
}
impl unconstrained::Builder for Builder {
type MasterProblem = CplexMaster;
fn build(&mut self) -> Result<CplexMaster> {
let mut cpx = CplexMaster::new()?;
cpx.max_bundle_size = self.max_bundle_size;
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
}
}