// Copyright (c) 2016-2022 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::saveable::Saveable;
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 thiserror::Error;

#[cfg(feature = "serialize")]
use serde::{Deserialize, Serialize};

use std::collections::VecDeque;
use std::error::Error;
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, Error)]
#[non_exhaustive]
pub enum CplexMasterError {
    #[error("CPLEX error")]
    Cplex(#[from] cpx::CplexError),
    #[error("master problem contains no minorants")]
    NoMinorants,
}

impl AsRef<dyn Error + 'static> for CplexMasterError {
    fn as_ref(&self) -> &(dyn Error + 'static) {
        self
    }
}

pub type Result<T> = std::result::Result<T, CplexMasterError>;

/// A minorant and its unique index.
#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
#[derive(Clone, Debug)]
struct MinorantInfo<Mn: Minorant> {
    minorant: Mn,
    index: usize,
}

/// Maps a global index to a minorant.
#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
#[derive(Clone, Copy, Debug)]
struct MinorantIdx {
    /// The function (subproblem) index.
    fidx: usize,
    /// The minorant index within the subproblem.
    ///
    /// The index is `None` if the minorant is still pending.
    idx: Option<usize>,
}

/// A submodel.
///
/// A submodel is a list of subproblems being aggregated in that model.
#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
#[derive(Debug, Clone, Default)]
struct SubModel {
    /// The list of subproblems.
    subproblems: Vec<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<(Real, DVector)>,
}

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<Mn: Minorant> {
    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<Mn>>>,

    /// A list of pending minorants.
    ///
    /// Minorants are only added to the submodel if there is a corresponding
    /// minorant for *each* subproblem in this submodel.
    pending_minorants: Vec<VecDeque<MinorantInfo<Mn>>>,

    /// 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: (Real, DVector),

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

#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
#[derive(Debug, Clone)]
pub struct CplexMasterState<Mn: Minorant> {
    num_subproblems: usize,
    force_update: bool,
    freeinds: Vec<usize>,
    updateinds: Vec<usize>,
    index2min: Vec<MinorantIdx>,
    qterm: Vec<DVector>,
    qdiag: Real,
    weight: Real,
    minorants: Vec<Vec<MinorantInfo<Mn>>>,
    pending_minorants: Vec<VecDeque<MinorantInfo<Mn>>>,
    submodels: Vec<SubModel>,
    in_submodel: Vec<usize>,
    opt_mults: Vec<DVector>,
    opt_minorant: (Real, DVector),
    max_bundle_size: usize,
}

unsafe impl<Mn: Minorant> Send for CplexMaster<Mn> {}

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

impl<Mn: Minorant> UnconstrainedMasterProblem<Mn> for CplexMaster<Mn> {
    type MinorantIndex = usize;

    type Err = CplexMasterError;

    fn new() -> Result<CplexMaster<Mn>> {
        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![],
            pending_minorants: vec![],
            select_model: Arc::new(|i| i),
            submodels: vec![],
            in_submodel: vec![],
            opt_mults: vec![],
            opt_minorant: Default::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::Barrier_ConvergeTol.to_c(),
            1e-12
        ));

        self.index2min.clear();
        self.freeinds.clear();
        self.minorants = (0..n).map(|_| Vec::new()).collect();
        self.pending_minorants = (0..n).map(|_| VecDeque::new()).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 {
        // TODO: is this correct? What about the pending minorants?
        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.submodels.len() {
            let n = self.submodels[i].minorants.len();
            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[self.submodels[i][0]][j]) as isize));
                let inds = &inds[self.max_bundle_size - 2..];

                for j in 0..self.submodels[i].len() {
                    let fidx = self.submodels[i][j];
                    // TODO: is this possible without the temporary created by `collect`?
                    let mins = inds.iter().map(|&m| self.minorants[fidx][m].index).collect::<Vec<_>>();
                    self.aggregate(fidx, &mins)?;
                }
            }
        }
        Ok(())
    }

    fn fill_models<F>(&mut self, f: F) -> Result<()>
    where
        F: FnMut(Self::MinorantIndex, &mut dyn Iterator<Item = (Self::MinorantIndex, Real)>),
    {
        let mut f = f;
        for i in 0..self.submodels.len() {
            if self.submodels[i].len() <= 1 {
                continue;
            }
            let mut min = std::usize::MAX;
            let mut max = std::usize::MIN;
            for &fidx in self.submodels[i].iter() {
                let m = self.minorants[fidx].len() + self.pending_minorants[fidx].len();
                min = min.min(m);
                max = max.max(m);
            }
            // if all subproblems have the same size or one contains not a
            // single minorant do nothing
            if 0 < min && min <= max {
                for j in 0..self.submodels[i].len() {
                    let fidx = self.submodels[i][j];
                    let nnew = max - (self.minorants[fidx].len() + self.pending_minorants[fidx].len());
                    let last_idx = self.minorants[fidx].len() - 1;
                    let last_index = self.minorants[fidx][last_idx].index;
                    for _ in 0..nnew {
                        let newindex = CplexMaster::add_minorant(
                            self,
                            fidx,
                            self.minorants[fidx][last_idx].minorant.clone(),
                            None,
                        )?;
                        f(newindex, &mut once((last_index, 1.0)));
                    }
                }

                debug_assert!(
                    self.submodels[i].iter().all(|&fidx| self.minorants[fidx].len() == max),
                    "Wrong lengths: {:?}",
                    self.submodels[i]
                        .iter()
                        .map(|&fidx| (fidx, self.minorants[fidx].len()))
                        .collect::<Vec<_>>()
                );
            }
        }

        Ok(())
    }

    fn add_minorant(&mut self, fidx: usize, minorant: Mn) -> Result<usize> {
        CplexMaster::add_minorant(self, fidx, minorant, None)
    }

    fn add_vars<SErr>(
        &mut self,
        nnew: usize,
        sgext: &dyn SubgradientExtender<Mn, SErr>,
    ) -> std::result::Result<(), Either<CplexMasterError, SErr>> {
        debug_assert!(!self.minorants[0].is_empty());
        if nnew == 0 {
            return Ok(());
        }

        for (fidx, mins) in self
            .minorants
            .iter_mut()
            .zip(self.pending_minorants.iter_mut())
            .map(|(mins, pmins)| mins.iter_mut().chain(pmins.iter_mut()))
            .enumerate()
        {
            for m in mins {
                sgext.extend_subgradient(fidx, &mut m.minorant).map_err(Either::Right)?;
            }
        }

        // update qterm because all minorants have changed
        self.force_update = true;
        // TODO: this might be inefficient, but we need to update all submodels
        // (can we do this lazily?)
        self.update_submodels();

        Ok(())
    }

    fn solve(&mut self, eta: &DVector, _fbound: Real, _augbound: Real, _relprec: Real) -> Result<()> {
        // all (pending) minorants should have been added to the model
        debug_assert!(self
            .submodels
            .iter()
            .all(|submodel| submodel.iter().all(|&fidx| self.pending_minorants[fidx].is_empty())));

        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 m in &submodel.minorants {
                    let cost = -m.constant() * self.weight - m.linear(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);

        for submodel in &self.submodels {
            debug_assert!(submodel
                .iter()
                .all(|&fidx| self.opt_mults[fidx].len() == submodel.minorants.len()));
            for i in 0..submodel.minorants.len() {
                for &fidx in submodel.iter() {
                    self.opt_mults[fidx][i] = sol[idx];
                    mults.push((sol[idx], &self.minorants[fidx][i].minorant));
                }
                idx += 1;
            }
        }

        Mn::combine_to_vec(mults, &mut self.opt_minorant);

        Ok(())
    }

    fn dualopt(&self) -> &DVector {
        &self.opt_minorant.1
    }

    fn dualopt_cutval(&self) -> Real {
        self.opt_minorant.constant()
    }

    fn multiplier(&self, min: usize) -> Real {
        let MinorantIdx { fidx, idx } = self.index2min[min];
        idx.map(|i| self.opt_mults[fidx][i]).unwrap_or(0.0)
    }

    fn eval_model(&self, y: &DVector) -> Real {
        let result = self
            .submodels
            .iter()
            .map(|s| {
                let mut value = NEG_INFINITY;
                for m in &s.minorants {
                    value = value.max(m.eval(y));
                }
                value
            })
            .sum::<Real>();

        // Verify that the value of the single models and the aggregated
        // submodels are the same.

        // debug_assert!(
        //     result
        //         <= (0..self.submodels.len())
        //             .map(|i| self.eval_model_subproblem(i, y))
        //             .sum::<Real>(),
        //     "Value of single models: {}, aggregated: {}",
        //     result,
        //     (0..self.submodels.len())
        //         .map(|i| self.eval_model_subproblem(i, y))
        //         .sum::<Real>()
        // );

        result
    }

    fn eval_model_subproblem(&self, fidx: usize, y: &DVector) -> Real {
        let mut result = NEG_INFINITY;
        for minfo in &self.minorants[fidx] {
            result = result.max(minfo.minorant.eval(y));
        }
        result
    }

    fn aggregated_primal(&self, fidx: usize) -> Result<Mn::Primal> {
        Ok(Aggregatable::combine(
            self.opt_mults[fidx]
                .iter()
                .enumerate()
                .map(|(i, alpha)| (*alpha, self.minorants[fidx][i].minorant.primal())),
        ))
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        for mins in &mut self.minorants {
            for m in mins.iter_mut() {
                m.minorant.move_center_begin(alpha, d);
            }
        }
        for mins in &mut self.pending_minorants {
            for m in mins.iter_mut() {
                m.minorant.move_center_begin(alpha, d);
            }
        }
        for submod in &mut self.submodels {
            for m in &mut submod.minorants {
                m.move_center_begin(alpha, d);
            }
        }
    }
}

impl<Mn: Minorant> CplexMaster<Mn> {
    /// 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
                .minorants
                .resize_with(minorants[submodel[0]].len(), Default::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 {
                    debug_assert!(idx.is_some());
                    debug_assert!(idx.map(|idx| idx < self.minorants[fidx].len()).unwrap_or(false));
                    if let Some(idx) = idx {
                        if idx
                            < submodel
                                .iter()
                                .map(|&fidx| self.minorants[fidx].len())
                                .min()
                                .unwrap_or(0)
                        {
                            return Some((index, mod_i, idx));
                        }
                    }
                }
                None
            })
            .collect::<Vec<_>>();

        // Compute the aggregated minorants.
        for &(_, mod_i, i) in &updateinds {
            let submodel = &mut self.submodels[mod_i];
            Mn::combine_to_vec(
                submodel
                    .subproblems
                    .iter()
                    .map(|&fidx| (1.0, &minorants[fidx][i].minorant)),
                &mut submodel.minorants[i],
            );
        }

        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.minorants.len() {
                // 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].1;
                // 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(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.minorants.len() {
                    let idx_i = self.minorants[submod_i[0]][i].index;
                    for submod_j in submodels.iter() {
                        for j in 0..submod_j.minorants.len() {
                            let idx_j = self.minorants[submod_j[0]][j].index;
                            let x = submod_i.minorants[i].linear(&submod_j.minorants[j].1);
                            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.minorants.len() {
                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.minorants.len()));
            rmatval.extend(repeat(1.0).take(submodel.minorants.len()));
            nvars += submodel.minorants.len();
        }

        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.minorants.len() {
                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.minorants.len() {
                        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");
        debug_assert!(mins.iter().all(|&i| self.index2min[i].idx.is_some()));

        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.index2min[i]
                .idx
                .map(|idx| self.opt_mults[fidx][idx])
                .unwrap_or(0.0);
        }
        let aggr_coeffs = if sum_coeffs != 0.0 {
            mins.iter()
                .map(|&i| {
                    self.index2min[i]
                        .idx
                        .map(|idx| self.opt_mults[fidx][idx])
                        .unwrap_or(0.0)
                        / 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())
                .filter_map(|(alpha, &index)| {
                    self.index2min[index]
                        .idx
                        .map(|idx| (alpha, &self.minorants[fidx][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 min_idx = min_idx.unwrap();
            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);

            if let Some(i) = self.updateinds.iter().position(|idx| *idx == m.index) {
                self.updateinds.remove(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() {
                debug_assert!(self.index2min[self.minorants[fidx][min_idx].index].idx.is_some());
                self.index2min[self.minorants[fidx][min_idx].index].idx = Some(min_idx);
                self.updateinds.push(self.minorants[fidx][min_idx].index);
            }
        }

        // Finally add the aggregated minorant.
        let aggr_idx = CplexMaster::add_minorant(self, fidx, aggr, Some(sum_coeffs))?;
        Ok((aggr_idx, aggr_coeffs))
    }

    fn add_minorant(&mut self, fidx: usize, minorant: Mn, alpha: Option<Real>) -> Result<usize> {
        debug!(
            "Add minorant fidx={} index={}: {:?}",
            fidx,
            self.minorants[fidx].len(),
            minorant.debug()
        );

        debug_assert!(
            self.minorants[fidx]
                .last()
                .map_or(true, |m| m.minorant.dim() == minorant.dim()),
            "Newly added minorant has wrong size (got: {}, expected: {})",
            minorant.dim(),
            self.minorants[fidx].last().map_or(0, |m| m.minorant.dim())
        );

        // find new unique minorant index
        //
        // the new minorant is not added immediately but put on the queue
        // until there is a corresponding minorant for all subproblems of the same
        // submodel
        let index = if let Some(index) = self.freeinds.pop() {
            self.index2min[index] = MinorantIdx { fidx, idx: None };
            index
        } else {
            self.index2min.push(MinorantIdx { fidx, idx: None });
            self.index2min.len() - 1
        };

        // enqueue minorant
        if let Some(alpha) = alpha {
            let min_idx = self.minorants[fidx].len();
            self.updateinds.push(index);
            self.index2min[index].idx = Some(min_idx);
            self.minorants[fidx].push(MinorantInfo { minorant, index });
            self.opt_mults[fidx].push(alpha);
            return Ok(index);
        }
        self.pending_minorants[fidx].push_back(MinorantInfo { minorant, index });

        // check if there are enough minorants to form a new one.
        let submodel = &self.submodels[self.in_submodel[fidx]];
        if submodel.iter().all(|&f| !self.pending_minorants[f].is_empty()) {
            for &f in submodel.iter() {
                let min_idx = self.minorants[f].len();
                let min = self.pending_minorants[f].pop_front().unwrap();
                debug!(
                    "Use pending ({}) minorant fidx:{} idx:{}, index:{}",
                    self.pending_minorants[f].len(),
                    f,
                    min_idx,
                    min.index
                );
                self.updateinds.push(min.index);
                self.index2min[min.index].idx = Some(min_idx);
                self.minorants[f].push(min);
                self.opt_mults[f].push(0.0);
            }
        }

        Ok(index)
    }
}

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<Mn: Minorant> super::Builder<Mn> for Builder {
    type MasterProblem = CplexMaster<Mn>;

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

impl<Mn: Minorant> Saveable for CplexMaster<Mn> {
    type State = CplexMasterState<Mn>;

    type Err = CplexMasterError;

    fn get_state(&self) -> std::result::Result<Self::State, Self::Err> {
        Ok(CplexMasterState {
            num_subproblems: self.num_subproblems(),
            force_update: self.force_update,
            freeinds: self.freeinds.clone(),
            updateinds: self.updateinds.clone(),
            index2min: self.index2min.clone(),
            qterm: self.qterm.clone(),
            qdiag: self.qdiag,
            weight: self.weight,
            minorants: self.minorants.clone(),
            pending_minorants: self.pending_minorants.clone(),
            submodels: self.submodels.clone(),
            in_submodel: self.in_submodel.clone(),
            opt_mults: self.opt_mults.clone(),
            opt_minorant: self.opt_minorant.clone(),
            max_bundle_size: self.max_bundle_size,
        })
    }

    fn set_state(&mut self, state: Self::State) -> std::result::Result<(), Self::Err> {
        self.set_num_subproblems(state.num_subproblems)?;
        self.force_update = state.force_update;
        self.freeinds = state.freeinds;
        self.updateinds = state.updateinds;
        self.index2min = state.index2min;
        self.qterm = state.qterm;
        self.qdiag = state.qdiag;
        self.weight = state.weight;
        self.minorants = state.minorants;
        self.pending_minorants = state.pending_minorants;
        self.submodels = state.submodels;
        self.in_submodel = state.in_submodel;
        self.opt_mults = state.opt_mults;
        self.opt_minorant = state.opt_minorant;
        self.max_bundle_size = state.max_bundle_size;
        Ok(())
    }
}
