// 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 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 either::Either;
use log::{debug, warn};
use std;
use std::f64::NEG_INFINITY;
use std::iter::{once, 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),
NoMinorants,
}
impl std::fmt::Display for CplexMasterError {
fn fmt(&self, fmt: &mut std::fmt::Formatter) -> std::fmt::Result {
use CplexMasterError::*;
match self {
Cplex(err) => err.fmt(fmt),
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),
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>;
/// 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<usize>,
/// 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<Minorant>,
}
impl Deref for SubModel {
type Target = Vec<usize>;
fn deref(&self) -> &Vec<usize> {
&self.subproblems
}
}
impl DerefMut for SubModel {
fn deref_mut(&mut self) -> &mut Vec<usize> {
&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<usize>,
/// List of minorant indices to be updated.
updateinds: Vec<usize>,
/// Mapping index to minorant.
index2min: Vec<MinorantIdx>,
/// The quadratic term.
qterm: Vec<DVector>,
/// 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<Vec<MinorantInfo>>,
/// The callback for the submodel for each subproblem.
select_model: Arc<dyn Fn(usize) -> usize>,
/// For each submodel the list of subproblems contained in that model.
submodels: Vec<SubModel>,
/// For each subproblem the submodel it is contained in.
in_submodel: Vec<usize>,
/// 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![],
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<F>(&mut self, f: F) -> Result<()>
where
F: FnMut(Self::MinorantIndex, &mut dyn 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_unstable_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::<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);
debug_assert!(
self.minorants[fidx]
.last()
.map_or(true, |m| m.linear.len() == minorant.linear.len()),
"Newly added minorant has wrong size (got: {}, expected: {})",
minorant.linear.len(),
self.minorants[fidx].last().map_or(0, |m| m.linear.len())
);
// 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<S>(
&mut self,
nnew: usize,
changed: &[usize],
mut extend_subgradient: S,
) -> std::result::Result<(), Either<CplexMasterError, S::Err>>
where
S: SubgradientExtension<Self::MinorantIndex>,
{
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
.subgradient_extension(fidx, m.index, &changedvars)
.map_err(Either::Right)?;
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::<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 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 solstat = unsafe { cpx::getstat(self.env, self.lp) };
if solstat != cpx::Stat::Optimal.to_c() && solstat != cpx::Stat::NumBest.to_c() {
warn!("Problem solving the master QP with Cplex (status: {:?})", solstat)
}
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<dyn Iterator<Item = (Self::MinorantIndex, Real)> + '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::<DVector>()
} else {
// All coefficients are zero, we can simply remove all but one of the minorants.
once(1.0).chain(repeat(0.0)).take(mins.len()).collect()
};
// 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_begin(alpha, d);
}
}
for submod in &mut self.submodels {
for m in &mut submod.minorants {
m.move_center_begin(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<F>(&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::<Vec<_>>();
// 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 <g_i, g_j> 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<c_char> = 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<dyn Fn(usize) -> 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<CplexMaster> {
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<F>(&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)
}
}