RsBundle  Changes On Branch restructure

 */

#![allow(non_upper_case_globals)]

//! Example implementation for a capacitated facility location problem.

use better_panic;
use crossbeam::channel::unbounded as channel;
use log::{info, Level};
use rustop::opts;
use std::error::Error;
use std::fmt::Write;
use std::io::Write as _;
use std::sync::Arc;

use env_logger::{self, fmt::Color};
use ordered_float::NotNan;
use threadpool::ThreadPool;

use bundle::parallel::{
    DefaultSolver as ParallelSolver, EvalResult, FirstOrderProblem as ParallelProblem,
    NoBundleSolver as NoParallelSolver, ResultSender,
use bundle::problem::{EvalResult, FirstOrderProblem as ParallelProblem, ResultSender};
use bundle::solver::sync::{DefaultSolver, NoBundleSolver};
};
use bundle::{dvec, DVector, Minorant, Real};
use bundle::{DefaultSolver, FirstOrderProblem, NoBundleSolver, SimpleEvaluation};

const Nfac: usize = 3;
const Ncus: usize = 5;
const F: [Real; Nfac] = [1000.0, 1000.0, 1000.0];
const CAP: [Real; Nfac] = [500.0, 500.0, 500.0];
const C: [[Real; Ncus]; Nfac] = [
    [4.0, 5.0, 6.0, 8.0, 10.0], //
    [6.0, 4.0, 3.0, 5.0, 8.0],  //
    [9.0, 7.0, 4.0, 3.0, 4.0],  //
];
const DEMAND: [Real; 5] = [80.0, 270.0, 250.0, 160.0, 180.0];

#[derive(Debug)]
enum EvalError {
    Customer(usize),
    NoObjective(usize),
}

impl std::fmt::Display for EvalError {
    fn fmt(&self, fmt: &mut std::fmt::Formatter) -> Result<(), std::fmt::Error> {
        match self {
            EvalError::Customer(i) => writeln!(fmt, "Customer subproblem {} failed", i),
            EvalError::NoObjective(i) => writeln!(fmt, "No objective value generated for subproblem {}", i),
        }
    }
}

impl Error for EvalError {}

struct CFLProblem {

            constant: objective,
            linear: subg,
        },
        primal,
    ))
}

impl FirstOrderProblem for CFLProblem {
    type Err = EvalError;

    type Primal = DVector;

    type EvalResult = SimpleEvaluation<Self::Primal>;

    fn num_variables(&self) -> usize {
        Nfac
    }

    fn lower_bounds(&self) -> Option<Vec<Real>> {
        Some(vec![0.0; FirstOrderProblem::num_variables(self)])
    }

    fn num_subproblems(&self) -> usize {
        Nfac + Ncus
    }

    fn evaluate(
        &mut self,
        index: usize,
        lambda: &[Real],
        _nullstep_bound: Real,
        _relprec: Real,
    ) -> Result<Self::EvalResult, Self::Err> {
        let (tx, rx) = channel();
        ParallelProblem::evaluate(self, index, Arc::new(lambda.iter().cloned().collect()), index, tx)?;
        let mut objective = None;
        let mut minorants = vec![];

        for r in rx {
            match r {
                Ok(EvalResult::ObjectiveValue { value, .. }) => objective = Some(value),
                Ok(EvalResult::Minorant { minorant, primal, .. }) => minorants.push((minorant, primal)),
                _ => break,
            }
        }

        Ok(SimpleEvaluation {
            objective: objective.ok_or(EvalError::NoObjective(index))?,
            minorants,
        })
    }
}

impl ParallelProblem for CFLProblem {
    type Err = EvalError;

    type Primal = DVector;

    fn num_variables(&self) -> usize {
        Nfac

    let (args, _) = opts! {
        synopsis "Solver a simple capacitated facility location problem";
        opt minimal:bool, desc:"Use the minimal master model";
    }
    .parse_or_exit();
    if !args.minimal {
        let mut slv = DefaultSolver::new(CFLProblem::new())?;
        let mut slv = DefaultSolver::<_>::new(CFLProblem::new());
        slv.params.max_bundle_size = 5;
        slv.terminator.termination_precision = 1e-9;
        slv.solve()?;
        show_primals(|i| slv.aggregated_primals(i))?;

        let mut slv = ParallelSolver::<_>::new(CFLProblem::new());
        slv.terminator.termination_precision = 1e-9;
        slv.master.max_bundle_size = 5;
        slv.solve()?;

        show_primals(|i| slv.aggregated_primal(i).unwrap())?;
    } else {
        let mut slv = NoBundleSolver::new(CFLProblem::new())?;
        let mut slv = NoBundleSolver::<_>::new(CFLProblem::new());
        slv.params.max_bundle_size = 2;
        slv.terminator.termination_precision = 1e-5;
        slv.solve()?;
        show_primals(|i| slv.aggregated_primals(i))?;

        let mut slv = NoParallelSolver::<_>::new(CFLProblem::new());
        slv.terminator.termination_precision = 1e-5;
        slv.solve()?;

        show_primals(|i| slv.aggregated_primal(i).unwrap())?;
    }

    Ok(())
}

use env_logger;
use env_logger::fmt::Color;
use log::{info, Level};
use rustop::opts;
use std::io::Write;

use bundle::master::{Builder as MasterBuilder, MasterProblem};
use bundle::master::{Builder, FullMasterBuilder, MasterProblem, MinimalMasterBuilder};
use bundle::mcf::MMCFProblem;
use bundle::parallel;
use bundle::solver::sync::Solver;
use bundle::{terminator::StandardTerminator, weighter::HKWeighter};
use bundle::{DefaultSolver, FirstOrderProblem, NoBundleSolver, Solver};
use bundle::{FullMasterBuilder, MinimalMasterBuilder};

use std::error::Error;
use std::result::Result;

fn solve_standard<M>(mut slv: Solver<MMCFProblem, StandardTerminator, HKWeighter, M>) -> Result<(), Box<dyn Error>>
where
    M: MasterBuilder + Default,
    M::MasterProblem: MasterProblem<MinorantIndex = usize>,
{
    slv.weighter.set_weight_bounds(1e-1, 100.0);
    slv.terminator.termination_precision = 1e-6;
    slv.solve()?;

    let costs: f64 = (0..slv.problem().num_subproblems())
        .map(|i| {
            let aggr_primals = slv.aggregated_primals(i);
            slv.problem().get_primal_costs(i, &aggr_primals)
        })
        .sum();
    info!("Primal costs: {}", costs);
    Ok(())
}

fn solve_parallel<M>(master: M, mmcf: MMCFProblem) -> Result<(), Box<dyn Error>>
where
    M: MasterBuilder,
    M: Builder,
    M::MasterProblem: MasterProblem<MinorantIndex = usize>,
{
    let mut slv = parallel::Solver::<_, StandardTerminator, HKWeighter, M>::with_master(mmcf, master);
    let mut slv = Solver::<_, StandardTerminator, HKWeighter, M>::with_master(mmcf, master);
    slv.weighter.set_weight_bounds(1e-1, 100.0);
    slv.terminator.termination_precision = 1e-6;
    slv.solve()?;

    let costs: f64 = (0..slv.problem().num_subproblems())
        .map(|i| {
            let aggr_primals = slv.aggregated_primal(i).unwrap();

    }
    .parse_or_exit();

    let filename = args.file;
    info!("Reading instance: {}", filename);

    if !args.minimal {
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
        let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            if args.aggregated {
                mmcf.multimodel = false;
                mmcf.set_separate_constraints(false);
            } else {
                mmcf.multimodel = true;
                mmcf.set_separate_constraints(args.separate);
        mmcf.set_separate_constraints(args.separate);
            }

            let mut solver = DefaultSolver::new(mmcf)?;
            solver.params.max_bundle_size = if args.bundle_size <= 1 {
                if args.aggregated {
                    50
                } else {
                    5
                }
            } else {
                args.bundle_size
            };
            solve_standard(solver)?;
        }

        println!("---------------------------------");
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            mmcf.set_separate_constraints(args.separate);
            mmcf.multimodel = true;
        mmcf.multimodel = true;

            let mut master = FullMasterBuilder::default();
            if args.aggregated {
                master.max_bundle_size(if args.bundle_size <= 1 { 50 } else { args.bundle_size });
                master.use_full_aggregation();
            } else {
                master.max_bundle_size(if args.bundle_size <= 1 { 5 } else { args.bundle_size });
            }
            solve_parallel(master, mmcf)?;
        let mut master = FullMasterBuilder::default();
        if args.aggregated {
            master.max_bundle_size(if args.bundle_size <= 1 { 50 } else { args.bundle_size });
            master.use_full_aggregation();
        } else {
            master.max_bundle_size(if args.bundle_size <= 1 { 5 } else { args.bundle_size });
        }
        solve_parallel(master, mmcf)?;
        }
    } else {
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
        let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            mmcf.multimodel = false;

            let mut solver = NoBundleSolver::new(mmcf)?;
            solver.params.max_bundle_size = 2;
            solve_standard(solver)?;
        }

        println!("---------------------------------");
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            mmcf.set_separate_constraints(args.separate);
            mmcf.multimodel = true;
        mmcf.set_separate_constraints(args.separate);
        mmcf.multimodel = true;

            let master = MinimalMasterBuilder::default();
            solve_parallel(master, mmcf)?;
        let master = MinimalMasterBuilder::default();
        solve_parallel(master, mmcf)?;
        }
    }

    Ok(())
}

use env_logger::fmt::Color;
use log::{debug, Level};
use rustop::opts;
use std::io::Write;
use std::sync::Arc;
use std::thread;

use bundle::parallel::{
    DefaultSolver as ParallelSolver, EvalResult, FirstOrderProblem as ParallelProblem,
    NoBundleSolver as NoParallelSolver, ResultSender,
use bundle::problem::{EvalResult, FirstOrderProblem as ParallelProblem, ResultSender};
use bundle::solver::sync::{DefaultSolver, NoBundleSolver};
};
use bundle::{DVector, DefaultSolver, FirstOrderProblem, Minorant, NoBundleSolver, Real, SimpleEvaluation};
use bundle::{DVector, Minorant, Real};

#[derive(Clone)]
struct QuadraticProblem {
    a: [[Real; 2]; 2],
    b: [Real; 2],
    c: Real,
}

impl QuadraticProblem {
    fn new() -> QuadraticProblem {
        QuadraticProblem {
            a: [[5.0, 1.0], [1.0, 4.0]],
            b: [-12.0, -10.0],
            c: 3.0,
        }
    }
}

impl FirstOrderProblem for QuadraticProblem {
    type Err = Box<dyn Error + Send + Sync>;
    type Primal = ();
    type EvalResult = SimpleEvaluation<()>;

    fn num_variables(&self) -> usize {
        2
    }

    #[allow(unused_variables)]
    fn evaluate(
        &mut self,
        fidx: usize,
        x: &[Real],
        nullstep_bnd: Real,
        relprec: Real,
    ) -> Result<Self::EvalResult, Self::Err> {
        assert_eq!(fidx, 0);
        let mut objective = self.c;
        let mut g = dvec![0.0; 2];

        for i in 0..2 {
            g[i] += (0..2).map(|j| self.a[i][j] * x[j]).sum::<Real>();
            objective += x[i] * (g[i] + self.b[i]);
            g[i] = 2.0 * g[i] + self.b[i];
        }

        debug!("Evaluation at {:?}", x);
        debug!("  objective={}", objective);
        debug!("  subgradient={}", g);

        Ok(SimpleEvaluation {
            objective: objective,
            minorants: vec![(
                Minorant {
                    constant: objective,
                    linear: g,
                },
                (),
            )],
        })
    }
}

impl ParallelProblem for QuadraticProblem {
    type Err = Box<dyn Error + Send + Sync>;
    type Primal = ();

    fn num_variables(&self) -> usize {
        2
    }

    fn num_subproblems(&self) -> usize {
        1
    }

    fn start(&mut self) {}

    fn stop(&mut self) {}

    fn evaluate<I>(
        &mut self,
        i: usize,
        y: Arc<DVector>,
        fidx: usize,
        x: Arc<DVector>,
        index: I,
        tx: ResultSender<I, Self::Primal, Self::Err>,
    ) -> Result<(), Self::Err>
    where
        I: Send + Copy + 'static,
    {
        let y = y.clone();
        let mut p = self.clone();
        thread::spawn(move || match FirstOrderProblem::evaluate(&mut p, i, &y, 0.0, 0.0) {
            Ok(res) => {
                tx.send(Ok(EvalResult::ObjectiveValue {
                    index,
                    value: res.objective,
                }))
                .unwrap();
        let x = x.clone();
        let p = self.clone();
        thread::spawn(move || {
            assert_eq!(fidx, 0);
            let mut objective = p.c;
            let mut g = dvec![0.0; 2];

            for i in 0..2 {
                g[i] += (0..2).map(|j| p.a[i][j] * x[j]).sum::<Real>();
                objective += x[i] * (g[i] + p.b[i]);
                g[i] = 2.0 * g[i] + p.b[i];
            }

            debug!("Evaluation at {:?}", x);
            debug!("  objective={}", objective);
            debug!("  subgradient={}", g);

            tx.send(Ok(EvalResult::ObjectiveValue {
                index,
                value: objective,
            }))
            .unwrap();
                for (minorant, primal) in res.minorants {
                    tx.send(Ok(EvalResult::Minorant {
                        index,
                        minorant,
                        primal,
                    }))
                    .unwrap();
            tx.send(Ok(EvalResult::Minorant {
                index,
                minorant: Minorant {
                    constant: objective,
                    linear: g,
                },
                primal: (),
            }))
            .unwrap();
                }
            }
            Err(err) => tx.send(Err(err)).unwrap(),
        });
        Ok(())
    }
}

fn main() -> Result<(), Box<dyn Error>> {
    better_panic::install();

    let (args, _) = opts! {
        synopsis "Solver a simple quadratic optimization problem";
        opt minimal:bool, desc:"Use the minimal master model";
    }
    .parse_or_exit();

    {
        let f = QuadraticProblem::new();
        if !args.minimal {
            let mut solver = DefaultSolver::new(f).map_err(|e| format!("{}", e))?;
    let f = QuadraticProblem::new();
    if !args.minimal {
        let mut solver = DefaultSolver::<_>::new(f);
            solver.weighter.set_weight_bounds(1.0, 1.0);
            solver.solve().map_err(|e| format!("{}", e))?;
        } else {
            let mut solver = NoBundleSolver::new(f).map_err(|e| format!("{}", e))?;
        solver.solve().map_err(|e| format!("{}", e))?;
    } else {
        let mut solver = NoBundleSolver::<_>::new(f);
            solver.params.max_bundle_size = 2;
            solver.weighter.set_weight_bounds(1.0, 1.0);
            solver.solve().map_err(|e| format!("{}", e))?;
        }
    }
        solver.solve().map_err(|e| format!("{}", e))?;
    }


    println!("-------------------------");

    {
        let f = QuadraticProblem::new();
        if !args.minimal {
            let mut solver = ParallelSolver::<_>::new(f);
            solver.solve().map_err(|e| format!("{}", e))?;
        } else {
            let mut solver = NoParallelSolver::<_>::new(f);
            solver.solve().map_err(|e| format!("{}", e))?;
        }
    }

    Ok(())
}

/*
 * Copyright (c) 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/>
 */

//! Objects that can be combined linearly.

use super::Real;
use std::borrow::Borrow;

/// An aggregatable object.
pub trait Aggregatable: Default {
    /// Return a scaled version of `other`, i.e. `alpha * other`.
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: Borrow<Self>;

    /// Add a scaled version of `other` to `self`.
    ///
    /// This sets `self = self + alpha * other`.
    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: Borrow<Self>;

    /// Return $\sum\_{i=1}\^n alpha_i m_i$.
    ///
    /// If `aggregates` is empty return the default value.
    fn combine<I, A>(aggregates: I) -> Self
    where
        I: IntoIterator<Item = (Real, A)>,
        A: Borrow<Self>,
    {
        let mut it = aggregates.into_iter();
        let mut x;
        if let Some((alpha, y)) = it.next() {
            x = Self::new_scaled(alpha, y);
        } else {
            return Self::default();
        }

        for (alpha, y) in it {
            x.add_scaled(alpha, y);
        }

        x
    }
}

/// Implement for empty tuples.
impl Aggregatable for () {
    fn new_scaled<A>(_alpha: Real, _other: A) -> Self
    where
        A: Borrow<Self>,
    {
    }

    fn add_scaled<A>(&mut self, _alpha: Real, _other: A)
    where
        A: Borrow<Self>,
    {
    }
}

/// Implement for scalar values.
impl Aggregatable for Real {
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: Borrow<Self>,
    {
        alpha * other.borrow()
    }

    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: Borrow<Self>,
    {
        *self += alpha * other.borrow()
    }
}

/// Implement for vectors of aggregatable objects.
impl<T> Aggregatable for Vec<T>
where
    T: Aggregatable,
{
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: std::borrow::Borrow<Self>,
    {
        other
            .borrow()
            .iter()
            .map(|y| Aggregatable::new_scaled(alpha, y))
            .collect()
    }

    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: std::borrow::Borrow<Self>,
    {
        debug_assert_eq!(self.len(), other.borrow().len(), "Vectors must have the same size");
        for (ref mut x, y) in self.iter_mut().zip(other.borrow()) {
            x.add_scaled(alpha, y)
        }
    }
}

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

//! A linear minorant.

use super::{Aggregatable, DVector, Real};

use std::borrow::Borrow;
use std::fmt;

/// A linear minorant of a convex function.
///
/// A linear minorant of a convex function $f \colon \mathbb{R}\^n \to
/// \mathbb{R}$ is a linear function of the form
///
///   \\[ l \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto \langle g, x
///   \rangle + c \\]
///
/// such that $l(x) \le f(x)$ for all $x \in \mathbb{R}\^n$.
#[derive(Clone, Debug)]
pub struct Minorant {
    /// The constant term.
    pub constant: Real,

    /// The linear term.
    pub linear: DVector,
}

impl fmt::Display for Minorant {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "{} + y * {}", self.constant, self.linear)?;
        Ok(())
    }
}

impl Default for Minorant {
    fn default() -> Minorant {
        Minorant {
            constant: 0.0,
            linear: dvec![],
        }
    }
}

impl Minorant {
    /// Return a new 0 minorant.
    pub fn new(constant: Real, linear: Vec<Real>) -> Minorant {
        Minorant {
            constant,
            linear: DVector(linear),
        }
    }

    /**
     * Evaluate minorant at some point.
     *
     * This function computes $c + \langle g, x \rangle$ for this minorant
     *   \\[\ell \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto c + \langle g, x \rangle\\]
     * and the given point $x \in \mathbb{R}\^n$.
     */
    pub fn eval(&self, x: &DVector) -> Real {
        self.constant + self.linear.dot(x)
    }

    /**
     * Move the center of the minorant.
     */
    pub fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.constant += alpha * self.linear.dot(d);
    }
}

impl Aggregatable for Minorant {
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: Borrow<Self>,
    {
        let m = other.borrow();
        Minorant {
            constant: alpha * m.constant,
            linear: DVector::scaled(&m.linear, alpha),
        }
    }

    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: Borrow<Self>,
    {
        let m = other.borrow();
        self.constant += alpha * m.constant;
        self.linear.add_scaled(alpha, &m.linear);
    }
}

/*
 * Copyright (c) 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/>
 */

//! General types and data structures.

pub mod vector;
pub use vector::{DVector, Vector};

pub mod aggregatable;
pub use aggregatable::Aggregatable;

pub mod minorant;
pub use minorant::Minorant;

/// Type used for real numbers throughout the library.
pub type Real = f64;

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

//! Finite-dimensional sparse and dense vectors.

use crate::{Aggregatable, Real};
use std::fmt;
use std::ops::{Deref, DerefMut};
// use std::cmp::min;
use std::borrow::Borrow;
use std::iter::FromIterator;
use std::vec::IntoIter;

#[cfg(feature = "blas")]
use {openblas_src as _, rs_blas as blas, std::os::raw::c_int};

/// Type of dense vectors.
#[derive(Debug, Clone, PartialEq, Default)]
pub struct DVector(pub Vec<Real>);

impl Deref for DVector {
    type Target = Vec<Real>;

    fn deref(&self) -> &Vec<Real> {
        &self.0
    }
}

impl DerefMut for DVector {
    fn deref_mut(&mut self) -> &mut Vec<Real> {
        &mut self.0
    }
}

impl fmt::Display for DVector {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "(")?;
        for (i, x) in self.iter().enumerate() {
            if i > 0 {
                write!(f, ", ")?;
            }
            write!(f, "{}", x)?
        }
        write!(f, ")")?;
        Ok(())
    }
}

impl FromIterator<Real> for DVector {
    fn from_iter<I: IntoIterator<Item = Real>>(iter: I) -> Self {
        DVector(Vec::from_iter(iter))
    }
}

impl IntoIterator for DVector {
    type Item = Real;
    type IntoIter = IntoIter<Real>;

    fn into_iter(self) -> IntoIter<Real> {
        self.0.into_iter()
    }
}

/// Type of dense or vectors.
#[derive(Debug, Clone)]
pub enum Vector {
    /// A vector with dense storage.
    Dense(DVector),

    /**
     * A vector with sparse storage.
     *
     * For each non-zero element this vector stores an index and the
     * value of the element in addition to the size of the vector.
     */
    Sparse { size: usize, elems: Vec<(usize, Real)> },
}

impl fmt::Display for Vector {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        match *self {
            Vector::Dense(ref v) => write!(f, "{}", v),
            Vector::Sparse { size, ref elems } => {
                let mut it = elems.iter();
                write!(f, "{}:(", size)?;
                if let Some(&(i, x)) = it.next() {
                    write!(f, "{}:{}", i, x)?;
                    for &(i, x) in it {
                        write!(f, ", {}:{}", i, x)?;
                    }
                }
                write!(f, ")")
            }
        }
    }
}

impl DVector {
    /// Set all elements to 0.
    pub fn init0(&mut self, size: usize) {
        self.clear();
        self.extend((0..size).map(|_| 0.0));
    }

    /// Set self = factor * y.
    pub fn scal(&mut self, factor: Real, y: &DVector) {
        self.clear();
        self.extend(y.iter().map(|y| factor * y));
    }

    /// Return factor * self.
    pub fn scaled(&self, factor: Real) -> DVector {
        let mut x = DVector::default();
        x.scal(factor, self);
        x
    }

    /// Return the inner product with another vector.
    pub fn dot(&self, other: &DVector) -> Real {
        assert_eq!(self.len(), other.len());
        self.dot_begin(other)
    }

    /// Return the inner product with another vector.
    ///
    /// The inner product is computed on the smaller of the two
    /// dimensions. All other elements are assumed to be zero.
    pub fn dot_begin(&self, other: &DVector) -> Real {
        #[cfg(feature = "blas")]
        unsafe {
            blas::ddot(self.len().min(other.len()) as c_int, &self, 1, &other, 1)
        }
        #[cfg(not(feature = "blas"))]
        {
            self.iter().zip(other.iter()).map(|(x, y)| x * y).sum::<Real>()
        }
    }

    /// Add two vectors and store result in this vector.
    pub fn add(&mut self, x: &DVector, y: &DVector) {
        assert_eq!(x.len(), y.len());
        self.clear();
        self.extend(x.iter().zip(y.iter()).map(|(a, b)| a + b));
    }

    /// Add two vectors and store result in this vector.
    pub fn add_scaled(&mut self, alpha: Real, y: &DVector) {
        assert_eq!(self.len(), y.len());
        #[cfg(feature = "blas")]
        unsafe {
            blas::daxpy(self.len() as c_int, alpha, &y, 1, &mut self[..], 1)
        }
        #[cfg(not(feature = "blas"))]
        {
            for (x, y) in self.iter_mut().zip(y.iter()) {
                *x += alpha * y;
            }
        }
    }

    /// Add two vectors and store result in this vector.
    ///
    /// In contrast to `add_scaled`, the two vectors might have
    /// different sizes. The size of the resulting vector is the
    /// larger of the two vector sizes and the remaining entries of
    /// the smaller vector are assumed to be 0.0.
    pub fn add_scaled_begin(&mut self, alpha: Real, y: &DVector) {
        #[cfg(feature = "blas")]
        unsafe {
            let n = self.len();
            blas::daxpy(n.min(y.len()) as c_int, alpha, &y, 1, &mut self[..], 1);
        }
        #[cfg(not(feature = "blas"))]
        {
            for (x, y) in self.iter_mut().zip(y.iter()) {
                *x += alpha * y;
            }
        }
        let n = self.len();
        if n < y.len() {
            self.extend(y[n..].iter().map(|y| alpha * y));
        }
    }

    /// Return the 2-norm of this vector.
    pub fn norm2(&self) -> Real {
        #[cfg(feature = "blas")]
        unsafe {
            blas::dnrm2(self.len() as c_int, &self, 1)
        }
        #[cfg(not(feature = "blas"))]
        {
            self.iter().map(|x| x * x).sum::<Real>().sqrt()
        }
    }
}

impl Aggregatable for DVector {
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: Borrow<Self>,
    {
        DVector::scaled(&other.borrow(), alpha)
    }

    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: Borrow<Self>,
    {
        DVector::add_scaled(self, alpha, &other.borrow())
    }
}

impl Vector {
    /**
     * Return a sparse vector with the given non-zeros.
     */
    pub fn new_sparse(n: usize, indices: &[usize], values: &[Real]) -> Vector {
        assert_eq!(indices.len(), values.len());

        if indices.is_empty() {
            Vector::Sparse { size: n, elems: vec![] }
        } else {
            let mut ordered: Vec<_> = (0..n).collect();
            ordered.sort_by_key(|&i| indices[i]);
            assert!(*indices.last().unwrap() < n);
            let mut elems = Vec::with_capacity(indices.len());
            let mut last_idx = n;
            for i in ordered {
                let val = unsafe { *values.get_unchecked(i) };
                if val != 0.0 {
                    let idx = unsafe { *indices.get_unchecked(i) };
                    if idx != last_idx {
                        elems.push((idx, val));
                        last_idx = idx;
                    } else {
                        elems.last_mut().unwrap().1 += val;
                        if elems.last_mut().unwrap().1 == 0.0 {
                            elems.pop();
                            last_idx = n;
                        }
                    }
                }
            }
            Vector::Sparse { size: n, elems }
        }
    }

    /**
     * Convert vector to a dense vector.
     *
     * This function always returns a copy of the vector.
     */
    pub fn to_dense(&self) -> DVector {
        match *self {
            Vector::Dense(ref x) => x.clone(),
            Vector::Sparse { size: n, elems: ref xs } => {
                let mut v = vec![0.0; n];
                for &(i, x) in xs {
                    unsafe { *v.get_unchecked_mut(i) = x };
                }
                DVector(v)
            }
        }
    }
}

#[test]
fn test_add_scaled_begin() {
    let mut x = dvec![1.0; 5];
    let y = dvec![2.0; 7];
    x.add_scaled_begin(3.0, &y);
    assert_eq!(x, dvec![7.0, 7.0, 7.0, 7.0, 7.0, 6.0, 6.0]);
}

// Copyright (c) 2016, 2017, 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/>
//

//! Problem description of a first-order convex optimization problem.

use crate::solver::UpdateState;
use crate::{Aggregatable, Minorant, Real};

use std::result::Result;
use std::vec::IntoIter;

/**
 * Trait for results of an evaluation.
 *
 * An evaluation returns the function value at the point of evaluation
 * and one or more subgradients.
 *
 * The subgradients (linear minorants) can be obtained by iterating over the result. The
 * subgradients are centered around the point of evaluation.
 */
pub trait Evaluation<P>: IntoIterator<Item = (Minorant, P)> {
    /// Return the function value at the point of evaluation.
    fn objective(&self) -> Real;
}

/**
 * Simple standard evaluation result.
 *
 * This result consists of the function value and a list of one or
 * more minorants and associated primal information.
 */
pub struct SimpleEvaluation<P> {
    pub objective: Real,
    pub minorants: Vec<(Minorant, P)>,
}

impl<P> IntoIterator for SimpleEvaluation<P> {
    type Item = (Minorant, P);
    type IntoIter = IntoIter<(Minorant, P)>;

    fn into_iter(self) -> Self::IntoIter {
        self.minorants.into_iter()
    }
}

impl<P> Evaluation<P> for SimpleEvaluation<P> {
    fn objective(&self) -> Real {
        self.objective
    }
}

/// Problem update information.
///
/// The solver calls the `update` method of the problem regularly.
/// This method can modify the problem by adding (or removing)
/// variables. The possible updates are encoded in this type.
#[derive(Debug, Clone, Copy)]
pub enum Update {
    /// Add a variable with bounds.
    ///
    /// The initial value of the variable will be the feasible value
    /// closest to 0.
    AddVariable { lower: Real, upper: Real },
    /// Add a variable with bounds and initial value.
    AddVariableValue { lower: Real, upper: Real, value: Real },
    /// Change the current value of a variable. The bounds remain
    /// unchanged.
    MoveVariable { index: usize, value: Real },
}

/**
 * Trait for implementing a first-order problem description.
 *
 */
pub trait FirstOrderProblem {
    /// Error raised by this oracle.
    type Err;

    /// The primal information associated with a minorant.
    type Primal: Aggregatable;

    /// Custom evaluation result value.
    type EvalResult: Evaluation<Self::Primal>;

    /// Return the number of variables.
    fn num_variables(&self) -> usize;

    /**
     * Return the lower bounds on the variables.
     *
     * If no lower bounds a specified, $-\infty$ is assumed.
     *
     * The lower bounds must be less then or equal the upper bounds.
     */
    fn lower_bounds(&self) -> Option<Vec<Real>> {
        None
    }

    /**
     * Return the upper bounds on the variables.
     *
     * If no lower bounds a specified, $+\infty$ is assumed.
     *
     * The upper bounds must be greater than or equal the upper bounds.
     */
    fn upper_bounds(&self) -> Option<Vec<Real>> {
        None
    }

    /// Return the number of subproblems.
    fn num_subproblems(&self) -> usize {
        1
    }

    /**
     * Evaluate the i^th subproblem at the given point.
     *
     * The returned evaluation result must contain (an upper bound on)
     * the objective value at $y$ as well as at least one subgradient
     * centered at $y$.
     *
     * If the evaluation process reaches a lower bound on the function
     * value at $y$ and this bound is larger than $nullstep_bound$,
     * the evaluation may stop and return the lower bound and a
     * minorant. In this case the function value is guaranteed to be
     * large enough so that the new point is rejected as candidate.
     *
     * The returned objective value should be an upper bound on the
     * true function value within $relprec \cdot (\\|f(y)\\| + 1.0)$,
     * otherwise the returned objective should be the maximum of all
     * linear minorants at $y$.
     *
     * Note that `nullstep_bound` and `relprec` are usually only
     * useful if there is only a `single` subproblem.
     */
    fn evaluate(
        &mut self,
        i: usize,
        y: &[Real],
        nullstep_bound: Real,
        relprec: Real,
    ) -> Result<Self::EvalResult, Self::Err>;

    /// Return updates of the problem.
    ///
    /// The default implementation returns no updates.
    fn update(&mut self, _state: &UpdateState<Self::Primal>) -> Result<Vec<Update>, Self::Err> {
        Ok(vec![])
    }

    /// Return new components for a subgradient.
    ///
    /// The components are typically generated by some primal information. The
    /// corresponding primal along with its subproblem index is passed as a
    /// parameter.
    ///
    /// The default implementation fails because it should never be
    /// called.
    fn extend_subgradient(
        &mut self,
        _i: usize,
        _primal: &Self::Primal,
        _vars: &[usize],
    ) -> Result<Vec<Real>, Self::Err> {
        unimplemented!()
    }
}

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! Proximal bundle method implementation.

/// Type used for real numbers throughout the library.
pub type Real = f64;

#[macro_export]
macro_rules! dvec {
    ( $ elem : expr ; $ n : expr ) => { DVector(vec![$elem; $n]) };
    ( $ ( $ x : expr ) , * ) => { DVector(vec![$($x),*]) };
    ( $ ( $ x : expr , ) * ) => { DVector(vec![$($x,)*]) };
}

pub mod vector;
mod data;
pub use crate::vector::{DVector, Vector};

pub use data::{Aggregatable, DVector, Minorant, Real, Vector};
pub mod minorant;
pub use crate::minorant::{Aggregatable, Minorant};

pub mod firstorderproblem;
pub use crate::firstorderproblem::{Evaluation, FirstOrderProblem, SimpleEvaluation, Update};

pub mod solver;
pub mod problem;
pub use crate::solver::{BundleState, FullMasterBuilder, IterationInfo, Solver, SolverParams, Step, UpdateState};

pub mod parallel;
pub mod solver;

pub mod weighter;

pub mod terminator;

pub mod master;

pub mod mcf;

/// The minimal bundle builder.
pub type MinimalMasterBuilder = master::boxed::Builder<master::minimal::Builder>;

/// The default bundle solver with general master problem.
pub type DefaultSolver<P> = Solver<P, terminator::StandardTerminator, weighter::HKWeighter, FullMasterBuilder>;

/// A bundle solver with a minimal cutting plane model.
pub type NoBundleSolver<P> = Solver<P, terminator::StandardTerminator, weighter::HKWeighter, MinimalMasterBuilder>;

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

pub mod unconstrained;
use crate::master::{self, MasterProblem, SubgradientExtension, UnconstrainedMasterProblem};
use self::unconstrained::UnconstrainedMasterProblem;

use super::MasterProblem;
pub use super::SubgradientExtension;
use crate::{DVector, Minorant, Real};

use itertools::multizip;
use log::debug;
use std::f64::{EPSILON, INFINITY, NEG_INFINITY};

/**

/// Builder for `BoxedMasterProblem`.
///
/// `B` is a builder of the underlying `UnconstrainedMasterProblem`.
#[derive(Default)]
pub struct Builder<B>(B);

impl<B> master::Builder for Builder<B>
impl<B> super::Builder for Builder<B>
where
    B: master::unconstrained::Builder,
    B: unconstrained::Builder,
    B::MasterProblem: UnconstrainedMasterProblem,
{
    type MasterProblem = BoxedMasterProblem<B::MasterProblem>;

    fn build(&mut self) -> Result<Self::MasterProblem, <Self::MasterProblem as MasterProblem>::Err> {
        self.0.build().map(BoxedMasterProblem::with_master)
    }

// Copyright (c) 2016, 2017, 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/>
//

pub mod cpx;
pub mod minimal;

use crate::{DVector, Minorant, Real};

pub use super::SubgradientExtension;

use std::error::Error;

/**
 * Trait for master problems without box constraints.
 *
 * Implementors of this trait are supposed to solve quadratic
 * optimization problems of the form
 *
 * \\[ \min \left\\{ \hat{f}(d) + \frac{u}{2} \\| d \\|\^2 \colon
 *     d \in \mathbb{R}\^n \right\\}. \\]
 *
 * where $\hat{f}$ is a piecewise linear model, i.e.
 *
 * \\[ \hat{f}(d) = \max \\{ \ell_i(d) = c_i + \langle g_i, d \rangle \colon
 *                           i=1,\dotsc,k \\}
 *                = \max \left\\{ \sum_{i=1}\^k \alpha_i \ell_i(d) \colon
 *                                \alpha \in \Delta \right\\}, \\]
 *
 * where $\Delta := \left\\{ \alpha \in \mathbb{R}\^k \colon \sum_{i=1}\^k
 * \alpha_i = 1 \right\\}$. Note, the unconstrained solver is expected
 * to compute *dual* optimal solutions, i.e. the solver must compute
 * optimal coefficients $\bar{\alpha}$ for the dual problem
 *
 * \\[ \max_{\alpha \in \Delta} \min_{d \in \mathbb{R}\^n}
 *     \sum_{i=1}\^k \alpha_i \ell_i(d) + \frac{u}{2} \\| d \\|\^2. \\]
 */
pub trait UnconstrainedMasterProblem: Send + 'static {
    /// Unique index for a minorant.
    type MinorantIndex: Copy + Eq;

    /// Error type for this master problem.
    type Err: Error + Send + Sync;

    /// Return a new instance of the unconstrained master problem.
    fn new() -> Result<Self, Self::Err>
    where
        Self: Sized;

    /// Return the number of subproblems.
    fn num_subproblems(&self) -> usize;

    /// Set the number of subproblems (different function models.)
    fn set_num_subproblems(&mut self, n: usize) -> Result<(), Self::Err>;

    /// Return the current weight.
    fn weight(&self) -> Real;

    /// Set the weight of the quadratic term, must be > 0.
    fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err>;

    /// Return the number of minorants of subproblem `fidx`.
    fn num_minorants(&self, fidx: usize) -> usize;

    /// Compress the bundle.
    ///
    /// When some minorants are compressed, the callback is called with the
    /// coefficients and indices of the compressed minorants and the index of
    /// the new minorant. The callback may be called several times.
    fn compress<F>(&mut self, f: F) -> Result<(), Self::Err>
    where
        F: FnMut(Self::MinorantIndex, &mut dyn Iterator<Item = (Self::MinorantIndex, Real)>);

    /// Add a new minorant to the model.
    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<Self::MinorantIndex, Self::Err>;

    /// Add or move some variables.
    ///
    /// The variables in `changed` have been changed, so the subgradient
    /// information must be updated. Furthermore, `nnew` new variables
    /// are added.
    fn add_vars(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<(), Self::Err>;

    /// Solve the master problem.
    fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<(), Self::Err>;

    /// Return the current dual optimal solution.
    fn dualopt(&self) -> &DVector;

    /// Return the current dual optimal solution value.
    fn dualopt_cutval(&self) -> Real;

    /// Return the multiplier associated with a minorant.
    fn multiplier(&self, min: Self::MinorantIndex) -> Real;

    /// Return the multipliers associated with a subproblem.
    fn opt_multipliers<'a>(&'a self, fidx: usize) -> Box<dyn Iterator<Item = (Self::MinorantIndex, Real)> + 'a>;

    /// Return the value of the current model at the given point.
    fn eval_model(&self, y: &DVector) -> Real;

    /// 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<(Self::MinorantIndex, DVector), Self::Err>;

    /// Move the center of the master problem along $\alpha \cdot d$.
    fn move_center(&mut self, alpha: Real, d: &DVector);
}

/// A builder for creating unconstrained master problem solvers.
pub trait Builder {
    /// The master problem to be build.
    type MasterProblem: UnconstrainedMasterProblem;

    /// Create a new master problem.
    fn build(&mut self) -> Result<Self::MasterProblem, <Self::MasterProblem as UnconstrainedMasterProblem>::Err>;
}

// 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 log::debug;

use std;
use std::f64::NEG_INFINITY;
use std::iter::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),
    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>;

/// 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_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);

        // 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(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<()> {
        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.index, &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 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 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 {
            dvec![0.0; mins.len()]
        };

        // 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(alpha, d);
            }
        }
        for submod in &mut self.submodels {
            for m in &mut submod.minorants {
                m.move_center(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)
    }
}