RsBundle  Artifact [0d98bfd044]

Artifact 0d98bfd0447c879000d5495798148eb86c4c95c4:

  • File src/master/boxed/unconstrained/cpx.rs — part of check-in [baa49a9473] at 2020-07-18 14:56:37 on branch minorant-primal — master: remove callback from `compress` method (user: fifr size: 27858) [more...]

// Copyright (c) 2016-2020 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::UnconstrainedMasterProblem;
use crate::problem::SubgradientExtender;
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<P: Aggregatable> {
    minorant: Minorant<P>,
    index: usize,
}

impl<P: Aggregatable> Deref for MinorantInfo<P> {
    type Target = Minorant<P>;
    fn deref(&self) -> &Minorant<P> {
        &self.minorant
    }
}

impl<P: Aggregatable> DerefMut for MinorantInfo<P> {
    fn deref_mut(&mut self) -> &mut Minorant<P> {
        &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<P: Aggregatable> {
    /// 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<P>>,
}

impl<P: Aggregatable> Deref for SubModel<P> {
    type Target = Vec<usize>;
    fn deref(&self) -> &Vec<usize> {
        &self.subproblems
    }
}

impl<P: Aggregatable> DerefMut for SubModel<P> {
    fn deref_mut(&mut self) -> &mut Vec<usize> {
        &mut self.subproblems
    }
}

pub struct CplexMaster<P: Aggregatable> {
    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<P>>>,

    /// 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<P>>,
    /// 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<P>,

    /// Maximal bundle size.
    pub max_bundle_size: usize,
}

unsafe impl<P: Aggregatable> Send for CplexMaster<P> {}

impl<P: Aggregatable> Drop for CplexMaster<P> {
    fn drop(&mut self) {
        unsafe { cpx::freeprob(self.env, &mut self.lp) };
        unsafe { cpx::closeCPLEX(&mut self.env) };
    }
}

impl<P: Aggregatable + 'static> UnconstrainedMasterProblem<P> for CplexMaster<P> {
    type MinorantIndex = usize;

    type Err = CplexMasterError;

    fn new() -> Result<CplexMaster<P>> {
        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) -> Result<()> {
        assert!(self.max_bundle_size >= 2, "Maximal bundle size must be >= 2");
        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<_>>();
                self.aggregate(i, &inds)?;
            }
        }
        Ok(())
    }

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant<P>) -> 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<SErr>(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtender<P, SErr>,
    ) -> std::result::Result<(), Either<CplexMasterError, SErr>> {
        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.primal, &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 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 aggregated_primal(&self, fidx: usize) -> Result<P> {
        Ok(P::combine(
            self.opt_mults[fidx]
                .iter()
                .enumerate()
                .map(|(i, alpha)| (*alpha, &self.minorants[fidx][i].primal)),
        ))
    }

    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<P: Aggregatable + 'static> CplexMaster<P> {
    /// 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(())
    }

    /// 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<(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))
    }
}

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<P: Aggregatable + 'static> super::Builder<P> for Builder {
    type MasterProblem = CplexMaster<P>;

    fn build(&mut self) -> Result<CplexMaster<P>> {
        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)
    }
}