RsBundle  Changes On Branch minorant-trait

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

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

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];

type Minorant = (Real, DVector, DVector);

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

impl std::fmt::Display for EvalError {

impl CFLProblem {
    fn new() -> CFLProblem {
        CFLProblem { pool: None }
    }
}

fn solve_facility_problem(j: usize, lambda: &[Real]) -> Result<(Real, Minorant<DVector>), EvalError> {
fn solve_facility_problem(j: usize, lambda: &[Real]) -> Result<(Real, Minorant), EvalError> {
    // facility subproblem max\{-F[j] + CAP[j] * lambda[j] * y[j] | y_j \in \{0,1\}\}
    let objective;
    let mut subg = dvec![0.0; lambda.len()];
    let primal;
    if -F[j] + CAP[j] * lambda[j] > 0.0 {
        objective = -F[j] + CAP[j] * lambda[j];
        subg[j] = CAP[j];
        primal = dvec![1.0; 1];
    } else {
        objective = 0.0;
        primal = dvec![0.0; 1];
    }

    Ok((
        objective,
    Ok((objective, (objective, subg, primal)))
        Minorant {
            constant: objective,
            linear: subg,
            primal,
        },
    ))
}

fn solve_customer_problem(i: usize, lambda: &[Real]) -> Result<(Real, Minorant<DVector>), EvalError> {
fn solve_customer_problem(i: usize, lambda: &[Real]) -> Result<(Real, Minorant), EvalError> {
    // customer subproblem
    let opt = (0..Nfac)
        .map(|j| (j, -DEMAND[i] * C[j][i] - DEMAND[i] * lambda[j]))
        .flat_map(|x| NotNan::new(x.1).map(|y| (x.0, y)))
        .max_by_key(|x| x.1)
        .ok_or(EvalError::Customer(i))?;

    let objective = opt.1.into_inner();
    let mut subg = dvec![0.0; lambda.len()];
    subg[opt.0] = -DEMAND[i];
    let mut primal = dvec![0.0; Nfac];
    primal[opt.0] = 1.0;

    Ok((
        objective,
    Ok((objective, (objective, subg, primal)))
        Minorant {
            constant: objective,
            linear: subg,
            primal,
        },
    ))
}

impl ParallelProblem for CFLProblem {
    type Err = EvalError;

    type Primal = DVector;
    type Minorant = Minorant;

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

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

use rustop::opts;
use std::io::Write;

use bundle::master::{Builder, FullMasterBuilder, MasterProblem, MinimalMasterBuilder};
use bundle::mcf::MMCFProblem;
use bundle::solver::asyn::Solver as AsyncSolver;
use bundle::solver::sync::Solver as SyncSolver;
use bundle::DVector;
use bundle::{terminator::StandardTerminator, weighter::HKWeighter};
use bundle::{DVector, Real};

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

type Minorant = (Real, DVector, DVector);

fn solve_parallel<M>(master: M, mmcf: MMCFProblem) -> Result<(), Box<dyn Error>>
where
    M: Builder<DVector>,
    M::MasterProblem: MasterProblem<DVector, MinorantIndex = usize>,
    M: Builder<Minorant>,
    M::MasterProblem: MasterProblem<Minorant, MinorantIndex = usize>,
{
    let mut slv = AsyncSolver::<_, 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();
            slv.problem().get_primal_costs(i, &[aggr_primals])
        })
        .sum();
    info!("Primal costs: {}", costs);
    Ok(())
}

fn solve_sync<M>(master: M, mmcf: MMCFProblem) -> Result<(), Box<dyn Error>>
where
    M: Builder<DVector>,
    M::MasterProblem: MasterProblem<DVector, MinorantIndex = usize>,
    M: Builder<Minorant>,
    M::MasterProblem: MasterProblem<Minorant, MinorantIndex = usize>,
{
    let mut slv = SyncSolver::<_, 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())

use std::error::Error;
use std::io::Write;
use std::sync::Arc;
use std::thread;

use bundle::problem::{FirstOrderProblem as ParallelProblem, ResultSender};
use bundle::solver::sync::{DefaultSolver, NoBundleSolver};
use bundle::{DVector, Minorant, Real};
use bundle::{DVector, 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 ParallelProblem for QuadraticProblem {
    type Err = Box<dyn Error + Send>;
    type Primal = ();

    type Minorant = (Real, DVector);

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

    fn num_subproblems(&self) -> usize {
        1

            }

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

            tx.objective(objective).unwrap();
            tx.minorant(Minorant {
            tx.minorant((objective, g)).unwrap();
                constant: objective,
                linear: g,
                primal: (),
            })
            .unwrap();
        });
        Ok(())
    }
}

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

//! Objects that can be combined linearly.

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

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

    /// 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>;
        A: Borrow<T>;

    /// 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>,
        A: Borrow<T>,
    {
        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();

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

/// Implement for pairs.
impl<T1, T2> Aggregatable for (T1, T2)
where
    T1: Aggregatable,
    T2: Aggregatable,
{
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: Borrow<Self>,
    {
        let x = other.borrow();
        (T1::new_scaled(alpha, &x.0), T2::new_scaled(alpha, &x.1))
    }

    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: Borrow<Self>,
    {
        let x = other.borrow();
        self.0.add_scaled(alpha, &x.0);
        self.1.add_scaled(alpha, &x.1);
    }
}

/// Implement for triples.
impl<T1, T2, T3> Aggregatable for (T1, T2, T3)
where
    T1: Aggregatable,
    T2: Aggregatable,
    T3: Aggregatable,
{
    fn new_scaled<A>(alpha: Real, other: A) -> Self
    where
        A: Borrow<Self>,
    {
        let x = other.borrow();
        (
            T1::new_scaled(alpha, &x.0),
            T2::new_scaled(alpha, &x.1),
            T3::new_scaled(alpha, &x.2),
        )
    }

    fn add_scaled<A>(&mut self, alpha: Real, other: A)
    where
        A: Borrow<Self>,
    {
        let x = other.borrow();
        self.0.add_scaled(alpha, &x.0);
        self.1.add_scaled(alpha, &x.1);
        self.2.add_scaled(alpha, &x.2);
    }
}

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

// 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 num_traits::Zero;
use std::borrow::Borrow;
use std::fmt;

use crate::data::raw::BLAS1;

/// 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$.
///
/// Each minorant may have an additional "primal" data, that is
/// linearly combined along with the minorant itself. This data can
/// be used for separation algorithms.
#[derive(Clone)]
pub struct Minorant<P: Aggregatable> {
    /// The constant term.
    pub constant: Real,
pub trait Minorant: Aggregatable + Default + Send + 'static {
    type Primal: Aggregatable + Send;

    /// Return the constant value.
    fn constant(&self) -> Real;

    /// The linear term.
    /// Return the dimension of the linear term.
    pub linear: DVector,

    fn dim(&self) -> usize;
    /// The associated primal data.
    pub primal: P,
}


    /// Evaluate the linear term at some point.
impl<P: Aggregatable> fmt::Debug for Minorant<P> {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
    fn linear(&self, x: &DVector) -> Real;
        write!(f, "{:?}", (self.constant, &self.linear))?;
        Ok(())
    }
}


    /// Return the associated primal data.
    fn primal(&self) -> &Self::Primal;

    /// Compute the inner product with another minorant.
impl<P: Aggregatable> fmt::Display for Minorant<P> {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
    fn product(&self, other: &Self) -> Real;
        write!(f, "{} + y * {}", self.constant, self.linear)?;
        Ok(())
    }

}

    /// Add a scaled `self.linear` to `target`.
impl<P: Aggregatable> Default for Minorant<P> {
    fn default() -> Minorant<P> {
        Minorant {
            constant: 0.0,
            linear: dvec![],
    ///
    /// This sets `target = target + alpha * self.linear`.
    fn add_scaled_to(&self, alpha: Real, target: &mut [Real]);
            primal: P::default(),
        }

    }
}

    /// Compute linear combination of (linear terms of) the given minorants and
    /// store in a dense vector.
    fn combine_to_vec<I, A>(minorants: I, target: &mut (Real, DVector))
    where
        I: IntoIterator<Item = (Real, A)>,
        A: Borrow<Self>,
    {
impl<P: Aggregatable> Minorant<P> {
    /// Return a new 0 minorant.
        target.0 = Real::zero();
    pub fn new(constant: Real, linear: Vec<Real>, primal: P) -> Minorant<P> {
        Minorant {
            constant,
            linear: DVector(linear),
            primal,
        }
    }

    /**
     * Evaluate minorant at some point.
        let mut it = minorants.into_iter();
        if let Some((alpha, m)) = it.next() {
            let m = m.borrow();
            target.1.init0(m.dim());
            m.add_scaled_to(alpha, &mut target.1);
            target.0 += alpha * m.constant();
            for (alpha, m) in it {
                let m = m.borrow();
                m.add_scaled_to(alpha, &mut target.1);
                target.0 += alpha * m.constant();
            }
        } else {
            target.1.clear();
        }
    }

    /// 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)
    ///
    /// 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$.
    ///
    fn eval(&self, x: &DVector) -> Real {
        self.constant() + self.linear(x)
    }

    /**
     * Move the center of the minorant.
    /// Move the center of the minorant along direction `d`.
     */
    pub fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.constant += alpha * self.linear.dot(d);
    fn move_center(&mut self, alpha: Real, d: &DVector) {
        debug_assert_eq!(d.len(), self.dim(), "Minorant and direction must have same dimension");
        self.move_center_begin(alpha, d)
    }

    /// Move the center of the minorant.
    ///
    /// In contrast to [`move_center`] the vector `d` might be shorter than the
    /// linear term of the minorant. The missing elements are assumed to be
    /// zero.
    ///
    /// [`move_center`]: Minorant::move_center
    pub fn move_center_begin(&mut self, alpha: Real, d: &DVector) {
        debug_assert!(self.linear.len() >= d.len());
        self.constant += alpha * self.linear.dot_begin(d);
    }
}
    fn move_center_begin(&mut self, alpha: Real, d: &DVector);

    fn debug(&self) -> &dyn std::fmt::Debug;
}

impl<P: Aggregatable + Send + 'static> Minorant for (Real, DVector, P) {
    type Primal = P;

    fn dim(&self) -> usize {
        self.1.len()
    }

    fn constant(&self) -> Real {
        self.0
    }
impl<P: Aggregatable> Aggregatable for Minorant<P> {
    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),
            primal: P::new_scaled(alpha, &m.primal),
        }
    }

    fn linear(&self, x: &DVector) -> Real {
        self.1.dot(x)
    }

    fn primal(&self) -> &P {
        &self.2
    }

    fn product(&self, other: &Self) -> Real {
        self.1.dot(&other.1)
    }

    fn add_scaled_to(&self, alpha: Real, target: &mut [Real]) {
        debug_assert_eq!(self.1.len(), target.len());
        if !alpha.is_zero() {
            BLAS1::add_scaled(target, alpha, &self.1);
        }
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.0 += alpha * self.1.dot(d);
    }

    fn move_center_begin(&mut self, alpha: Real, d: &DVector) {
        debug_assert!(self.1.len() >= d.len());
        self.0 += alpha * self.1.dot_begin(d);
    }

    fn debug(&self) -> &dyn std::fmt::Debug {
        &self.1
    }
}

impl Minorant for (Real, DVector) {
    type Primal = ();

    fn dim(&self) -> usize {
        self.1.len()
    }
    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);
        self.primal.add_scaled(alpha, &m.primal);

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

    fn linear(&self, x: &DVector) -> Real {
        self.1.dot(x)
    }

    fn primal(&self) -> &Self::Primal {
        &()
    }

    fn product(&self, other: &Self) -> Real {
        self.1.dot(&other.1)
    }

    fn add_scaled_to(&self, alpha: Real, target: &mut [Real]) {
        debug_assert_eq!(self.1.len(), target.len());
        if !alpha.is_zero() {
            BLAS1::add_scaled(target, alpha, &self.1);
        }
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.0 += alpha * self.1.dot(d);
    }

    fn move_center_begin(&mut self, alpha: Real, d: &DVector) {
        debug_assert!(self.1.len() >= d.len());
        self.0 += alpha * self.1.dot_begin(d);
    }

    fn debug(&self) -> &dyn std::fmt::Debug {
        self
    }
}

/// The subgradient extender is a callback used to update existing minorants
/// The subgradient extender is a callback used to update an existing minorant.
/// given their associated primal data.
pub type SubgradientExtender<P, E> = dyn FnMut(usize, &P, &[usize]) -> Result<DVector, E> + Send;
pub type SubgradientExtender<M, E> = dyn FnMut(usize, &mut M) -> Result<(), E> + Send;

/*
 * Copyright (c) 2019 Frank Fischer <frank-fischer@shadow-soft.de>
 * Copyright (c) 2019, 2020 Frank Fischer <frank-fischer@shadow-soft.de>
 *
 * This program is free software: you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see  <http://www.gnu.org/licenses/>
 */

//! General types and data structures.

pub(crate) mod raw;

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) 2020 Frank Fischer <frank-fischer@shadow-soft.de>
 *
 * This program is free software: you can redistribute it and/or
 * modify it under the terms of the GNU General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see  <http://www.gnu.org/licenses/>
 */

///! BLAS-like low-level vector routines.

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

pub trait BLAS1<T> {
    /// Compute the inner product.
    fn dot(&self, y: &[T]) -> T;

    /// Compute the inner product.
    ///
    /// The inner product is computed on the smaller of the two
    /// dimensions. All other elements are assumed to be zero.
    fn dot_begin(&self, y: &[T]) -> T;

    /// Compute `self = self + alpha * y`.
    fn add_scaled(&mut self, alpha: T, y: &[T]);

    /// Return the 2-norm of this vector.
    fn norm2(&self) -> T;
}

impl BLAS1<f64> for [f64] {
    fn dot(&self, other: &[f64]) -> f64 {
        debug_assert_eq!(self.len(), other.len(), "Vectors must have the same size");
        Self::dot_begin(self, other)
    }

    fn dot_begin(&self, other: &[f64]) -> f64 {
        #[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::<f64>()
        }
    }

    fn add_scaled(&mut self, alpha: f64, y: &[f64]) {
        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;
            }
        }
    }

    fn norm2(&self) -> f64 {
        #[cfg(feature = "blas")]
        unsafe {
            blas::dnrm2(self.len() as c_int, &self, 1)
        }
        #[cfg(not(feature = "blas"))]
        {
            self.iter().map(|x| x * x).sum::<f64>().sqrt()
        }
    }
}

// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
// Copyright (c) 2016-2020 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but

use num_traits::Zero;
use std::borrow::Borrow;
use std::fmt;
use std::iter::FromIterator;
use std::ops::{Deref, DerefMut};
use std::vec::IntoIter;

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

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

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

    }

    /// 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)
        BLAS1::dot_begin(&self[..], other)
        }
        #[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)
        BLAS1::add_scaled(&mut self[..], alpha, y);
        }
        #[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();
        let n = self.len().min(y.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;
            }
        BLAS1::add_scaled(&mut self[..n], alpha, &y[..n]);

        }
        let n = self.len();
        if n < y.len() {
        if self.len() < 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)
        BLAS1::norm2(&self[..])
        }
        #[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>,

//

pub mod unconstrained;
use self::unconstrained::UnconstrainedMasterProblem;

use super::MasterProblem;
use crate::problem::SubgradientExtender;
use crate::{Aggregatable, DVector, Minorant, Real};
use crate::{DVector, Minorant, Real};

use either::Either;
use log::debug;
use std::f64::{EPSILON, INFINITY, NEG_INFINITY};
use std::marker::PhantomData;

/**
 * Turn unconstrained master problem into box-constrained one.
 *
 * This master problem adds box constraints to an unconstrained
 * master problem implementation. The box constraints are enforced by
 * an additional outer optimization loop.
 */
pub struct BoxedMasterProblem<P, M>
pub struct BoxedMasterProblem<Mn, M>
where
    P: Aggregatable,
    Mn: Minorant,
{
    /// The unconstrained master problem solver.
    pub master: M,

    /// Maximal number of updates of box multipliers.
    pub max_updates: usize,

    model_eps: Real,

    need_new_candidate: bool,

    /// Current number of updates.
    cnt_updates: usize,

    phantom: PhantomData<P>,
    phantom: PhantomData<Mn>,
}

impl<P, M> BoxedMasterProblem<P, M>
impl<Mn, M> BoxedMasterProblem<Mn, M>
where
    P: Aggregatable,
    M: UnconstrainedMasterProblem<P>,
    Mn: Minorant,
    M: UnconstrainedMasterProblem<Mn>,
{
    pub fn with_master(master: M) -> BoxedMasterProblem<P, M> {
    pub fn with_master(master: M) -> BoxedMasterProblem<Mn, M> {
        BoxedMasterProblem {
            master,
            max_updates: 50,
            lb: dvec![],
            ub: dvec![],
            eta: dvec![],
            primopt: dvec![],

     * the current box-multipliers $\eta$.
     */
    fn get_norm_subg2(&self) -> Real {
        self.eta.dot(self.master.dualopt())
    }
}

impl<P, M> MasterProblem<P> for BoxedMasterProblem<P, M>
impl<Mn, M> MasterProblem<Mn> for BoxedMasterProblem<Mn, M>
where
    P: Aggregatable + Send + 'static,
    M: UnconstrainedMasterProblem<P>,
    Mn: Minorant,
    M: UnconstrainedMasterProblem<Mn>,
{
    type MinorantIndex = M::MinorantIndex;

    type Err = M::Err;

    fn new() -> Result<Self, Self::Err> {
        M::new().map(BoxedMasterProblem::with_master)

        Ok(())
    }

    fn num_minorants(&self, fidx: usize) -> usize {
        self.master.num_minorants(fidx)
    }

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant<P>) -> Result<Self::MinorantIndex, Self::Err> {
    fn add_minorant(&mut self, fidx: usize, minorant: Mn) -> Result<Self::MinorantIndex, Self::Err> {
        self.master.add_minorant(fidx, minorant)
    }

    fn weight(&self) -> Real {
        self.master.weight()
    }

    fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err> {
        self.master.set_weight(weight)
    }

    fn add_vars<SErr>(
        &mut self,
        bounds: &[(Option<usize>, Real, Real)],
        extend_subgradient: &mut SubgradientExtender<P, SErr>,
        bounds: &[(Real, Real)],
        extend_subgradient: &mut SubgradientExtender<Mn, SErr>,
    ) -> Result<(), Either<Self::Err, SErr>> {
        debug_assert!(bounds.iter().all(|&(l, u)| l <= u));
        if !bounds.is_empty() {
            for (index, l, u) in bounds.iter().filter_map(|v| v.0.map(|i| (i, v.1, v.2))) {
                self.lb[index] = l;
                self.ub[index] = u;
            }
            self.lb.extend(bounds.iter().filter(|v| v.0.is_none()).map(|x| x.1));
            self.ub.extend(bounds.iter().filter(|v| v.0.is_none()).map(|x| x.2));
            self.lb.extend(bounds.iter().map(|x| x.0));
            self.ub.extend(bounds.iter().map(|x| x.1));
            self.eta.resize(self.lb.len(), 0.0);
            self.need_new_candidate = true;
            let nnew = bounds.iter().filter(|v| v.0.is_none()).count();
            let nnew = bounds.len();
            let changed = bounds.iter().filter_map(|v| v.0).collect::<Vec<_>>();
            self.master.add_vars(nnew, &changed, extend_subgradient)
            self.master.add_vars(nnew, extend_subgradient)
        } else {
            Ok(())
        }
    }

    #[allow(clippy::cognitive_complexity)]
    fn solve(&mut self, center_value: Real) -> Result<(), Self::Err> {

        self.dualoptnorm2
    }

    fn multiplier(&self, min: Self::MinorantIndex) -> Real {
        self.master.multiplier(min)
    }

    fn aggregated_primal(&self, fidx: usize) -> Result<P, Self::Err> {
    fn aggregated_primal(&self, fidx: usize) -> Result<Mn::Primal, Self::Err> {
        self.master.aggregated_primal(fidx)
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.need_new_candidate = true;
        self.master.move_center(alpha, d);
        self.lb.add_scaled_begin(-alpha, d);

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

impl<P, B> super::Builder<P> for Builder<B>
impl<Mn, B> super::Builder<Mn> for Builder<B>
where
    P: Aggregatable + Send + 'static,
    B: unconstrained::Builder<P>,
    B::MasterProblem: UnconstrainedMasterProblem<P>,
    Mn: Minorant,
    B: unconstrained::Builder<Mn>,
    B::MasterProblem: UnconstrainedMasterProblem<Mn>,
{
    type MasterProblem = BoxedMasterProblem<P, B::MasterProblem>;
    type MasterProblem = BoxedMasterProblem<Mn, B::MasterProblem>;

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

impl<B> std::ops::Deref for Builder<B> {
    type Target = B;

// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

pub mod cpx;
pub mod minimal;

use crate::problem::SubgradientExtender;
use crate::{Aggregatable, DVector, Minorant, Real};
use crate::{DVector, Minorant, Real};

use either::Either;
use std::error::Error;

/**
 * Trait for master problems without box constraints.
 *

 * \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<P: Aggregatable>: Send + 'static {
pub trait UnconstrainedMasterProblem<M: Minorant>: 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>

    ///
    /// 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(&mut self) -> Result<(), Self::Err>;

    /// Add a new minorant to the model.
    fn add_minorant(&mut self, fidx: usize, minorant: Minorant<P>) -> Result<Self::MinorantIndex, Self::Err>;
    fn add_minorant(&mut self, fidx: usize, minorant: M) -> 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
    /// Add `nnew` new variables.
    /// are added.
    fn add_vars<SErr>(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtender<P, SErr>,
        extend_subgradient: &mut SubgradientExtender<M, SErr>,
    ) -> Result<(), Either<Self::Err, SErr>>;

    /// 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 value of the current model at the given point.
    fn eval_model(&self, y: &DVector) -> Real;

    /// Return the aggregated primal.
    fn aggregated_primal(&self, fidx: usize) -> Result<P, Self::Err>;
    fn aggregated_primal(&self, fidx: usize) -> Result<M::Primal, 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<P: Aggregatable> {
pub trait Builder<Mn: Minorant> {
    /// The master problem to be build.
    type MasterProblem: UnconstrainedMasterProblem<P>;
    type MasterProblem: UnconstrainedMasterProblem<Mn>;

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

        CplexMasterError::Cplex(err)
    }
}

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

/// A minorant and its unique index.
struct MinorantInfo<P: Aggregatable> {
    minorant: Minorant<P>,
struct MinorantInfo<Mn: Minorant> {
    minorant: Mn,
    index: usize,
}

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

    fn new(minorant: Mn, index: usize) -> Self {
impl<P: Aggregatable> DerefMut for MinorantInfo<P> {
    fn deref_mut(&mut self) -> &mut Minorant<P> {
        &mut self.minorant
        MinorantInfo { minorant, index }
    }
}

/// Maps a global index to a minorant.
#[derive(Clone, Copy)]
struct MinorantIdx {
    /// The function (subproblem) index.
    fidx: usize,
    /// The minorant index within the subproblem.
    idx: usize,
}

/// A submodel.
///
/// A submodel is a list of subproblems being aggregated in that model.
#[derive(Debug, Clone, Default)]
struct SubModel<P: Aggregatable> {
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<P>>,
    minorants: Vec<(Real, DVector)>,
}

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

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

pub struct CplexMaster<P: Aggregatable> {
pub struct CplexMaster<Mn: Minorant> {
    env: *mut cpx::Env,

    lp: *mut cpx::Lp,

    /// True if the QP must be updated.
    force_update: bool,

    /// The additional diagonal term to ensure positive definiteness.
    qdiag: Real,

    /// The weight of the quadratic term.
    weight: Real,

    /// The minorants for each subproblem in the model.
    minorants: Vec<Vec<MinorantInfo<P>>>,
    minorants: Vec<Vec<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<P>>,
    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<P>,
    opt_minorant: (Real, DVector),

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

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

impl<P: Aggregatable> Drop for CplexMaster<P> {
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<P: Aggregatable + 'static> UnconstrainedMasterProblem<P> for CplexMaster<P> {
impl<Mn: Minorant> UnconstrainedMasterProblem<Mn> for CplexMaster<Mn> {
    type MinorantIndex = usize;

    type Err = CplexMasterError;

    fn new() -> Result<CplexMaster<P>> {
    fn new() -> Result<CplexMaster<Mn>> {
        let env;

        trycpx!({
            let mut status = 0;
            env = cpx::openCPLEX(&mut status);
            status
        });

            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(),
            opt_minorant: Default::default(),
            max_bundle_size: 50,
        })
    }

    fn num_subproblems(&self) -> usize {
        self.minorants.len()
    }

                    .collect::<Vec<_>>();
                self.aggregate(i, &inds)?;
            }
        }
        Ok(())
    }

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

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

        // 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.minorants[fidx].push(MinorantInfo::new(minorant, index));
        self.opt_mults[fidx].push(0.0);

        Ok(index)
    }

    fn add_vars<SErr>(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtender<P, SErr>,
        extend_subgradient: &mut SubgradientExtender<Mn, SErr>,
    ) -> std::result::Result<(), Either<CplexMasterError, SErr>> {
        debug_assert!(!self.minorants[0].is_empty());
        if changed.is_empty() && nnew == 0 {
        if nnew == 0 {
            return Ok(());
        }
        let noldvars = self.minorants[0][0].linear.len();
        let nnewvars = noldvars + nnew;

        let mut changedvars = vec![];
        changedvars.extend_from_slice(changed);
        changedvars.extend(noldvars..nnewvars);
        for (fidx, mins) in self.minorants.iter_mut().enumerate() {
            for m in &mut mins[..] {
                let new_subg = (extend_subgradient)(fidx, &m.primal, &changedvars).map_err(Either::Right)?;
                (extend_subgradient)(fidx, &mut m.minorant).map_err(Either::Right)?;
                for (&j, &g) in changed.iter().zip(new_subg.iter()) {
                    m.linear[j] = g;
                }
                m.linear.extend_from_slice(&new_subg[changed.len()..]);
            }
        }

        // update qterm because all minorants have changed
        self.force_update = true;

        Ok(())

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

            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));
        Mn::combine_to_vec(mults.into_iter().zip(mins), &mut self.opt_minorant);

        Ok(())
    }

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

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

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

    fn eval_model(&self, y: &DVector) -> Real {
        let mut result = 0.0;
        for submodel in &self.submodels {
            let mut this_val = NEG_INFINITY;
            for m in &submodel.minorants[0..submodel.num_mins] {
                this_val = this_val.max(m.eval(y));
            }
            result += this_val;
        }
        result
    }

    fn aggregated_primal(&self, fidx: usize) -> Result<P> {
        Ok(P::combine(
    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].primal)),
                .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.move_center_begin(alpha, d);
                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<P: Aggregatable + 'static> CplexMaster<P> {
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

        }

        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);
            submodel.minorants.resize_with(submodel.num_mins, 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

                }
            })
            .collect::<Vec<_>>();

        // Compute the aggregated minorants.
        for &(_, mod_i, i) in &updateinds {
            let submodel = &mut self.submodels[mod_i];
            Mn::combine_to_vec(
            submodel.minorants[i] =
                Aggregatable::combine(submodel.iter().map(|&fidx| (1.0, &minorants[fidx][i].minorant)));
                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.
        //

        // - 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;
                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.dot(g_i);
                    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.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);
                            let x = submod_i.minorants[i].linear(&submod_j.minorants[j].1);
                            debug_assert!((x - self.qterm[idx_i][idx_j]) < 1e-6);
                        }
                    }
                }
            }
        }

        Builder {
            max_bundle_size: 50,
            select_model: Arc::new(|i| i),
        }
    }
}

impl<P: Aggregatable + 'static> super::Builder<P> for Builder {
    type MasterProblem = CplexMaster<P>;
impl<Mn: Minorant> super::Builder<Mn> for Builder {
    type MasterProblem = CplexMaster<Mn>;

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

 * is that this model can be solved explicitely and very quickly, but
 * it is only a very loose approximation of the objective function.
 *
 * Because of its properties, it can only be used if the problem to be
 * solved has a maximal number of minorants of two and only one
 * subproblem.
 */
pub struct MinimalMaster<P: Aggregatable> {
pub struct MinimalMaster<Mn: Minorant> {
    /// The weight of the quadratic term.
    weight: Real,

    /// The minorants in the model.
    ///
    /// There are up to two minorants, each minorant consists of one part for
    /// each subproblem.
    ///
    /// The minorants for the i-th subproblem have the indices `2*i` and
    /// `2*i+1`.
    minorants: [Vec<Minorant<P>>; 2],
    minorants: [Vec<Mn>; 2],
    /// The number of minorants. Only the minorants with index less than this
    /// number are valid.
    num_minorants: usize,
    /// The number of minorants for each subproblem.
    num_minorants_of: Vec<usize>,
    /// The number of subproblems.
    num_subproblems: usize,
    /// The number of subproblems with at least 1 minorant.
    num_subproblems_with_1: usize,
    /// The number of subproblems with at least 2 minorants.
    num_subproblems_with_2: usize,
    /// Optimal multipliers.
    opt_mult: [Real; 2],
    /// Optimal aggregated minorant.
    opt_minorant: Minorant<P>,
    opt_minorant: (Real, DVector),
}

impl<P: Aggregatable + Send + 'static> UnconstrainedMasterProblem<P> for MinimalMaster<P> {
impl<Mn: Minorant> UnconstrainedMasterProblem<Mn> for MinimalMaster<Mn> {
    type MinorantIndex = usize;

    type Err = MinimalMasterError;

    fn new() -> Result<MinimalMaster<P>, Self::Err> {
    fn new() -> Result<MinimalMaster<Mn>, Self::Err> {
        Ok(MinimalMaster {
            weight: 1.0,
            num_minorants: 0,
            num_minorants_of: vec![],
            num_subproblems: 0,
            num_subproblems_with_1: 0,
            num_subproblems_with_2: 0,
            minorants: [vec![], vec![]],
            opt_mult: [0.0, 0.0],
            opt_minorant: Minorant::default(),
            opt_minorant: Default::default(),
        })
    }

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

    fn set_num_subproblems(&mut self, n: usize) -> Result<(), Self::Err> {
        self.num_subproblems = n;
        self.num_minorants = 0;
        self.num_minorants_of = vec![0; n];
        self.num_subproblems_with_1 = 0;
        self.num_subproblems_with_2 = 0;
        self.minorants = [
            (0..n).map(|_| Minorant::default()).collect(),
            (0..n).map(|_| Minorant::default()).collect(),
            (0..n).map(|_| Default::default()).collect(),
            (0..n).map(|_| Default::default()).collect(),
        ];
        Ok(())
    }

    fn compress(&mut self) -> Result<(), Self::Err> {
        if self.num_minorants == 2 {
            debug!("Aggregate");
            debug!(
                "  {} * {:?}",
            debug!("  {} * {:?}", self.opt_mult[0], self.minorants[0]);
            debug!("  {} * {:?}", self.opt_mult[1], self.minorants[1]);
                self.opt_mult[0],
                self.minorants[0].iter().map(|m| m.debug()).collect::<Vec<_>>()
            );
            debug!(
                "  {} * {:?}",
                self.opt_mult[1],
                self.minorants[1].iter().map(|m| m.debug()).collect::<Vec<_>>()
            );

            self.minorants[0] = Aggregatable::combine(self.opt_mult.iter().cloned().zip(&self.minorants));
            self.opt_mult[0] = 1.0;
            self.num_minorants = 1;
            self.num_minorants_of.clear();
            self.num_minorants_of.resize(self.num_subproblems, 1);
            self.num_subproblems_with_2 = 0;

            debug!("  {:?}", self.minorants[0]);
            debug!(
                "  {:?}",
                self.minorants[0].iter().map(|m| m.debug()).collect::<Vec<_>>()
            );
        }
        Ok(())
    }

    fn weight(&self) -> Real {
        self.weight
    }

    fn set_weight(&mut self, weight: Real) -> Result<(), Self::Err> {
        assert!(weight > 0.0);
        self.weight = weight;
        Ok(())
    }

    fn num_minorants(&self, fidx: usize) -> usize {
        self.num_minorants_of[fidx]
    }

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant<P>) -> Result<usize, Self::Err> {
    fn add_minorant(&mut self, fidx: usize, minorant: Mn) -> Result<usize, Self::Err> {
        if self.num_minorants_of[fidx] >= 2 {
            return Err(MinimalMasterError::MaxMinorants { subproblem: fidx });
        }

        let minidx = self.num_minorants_of[fidx];
        self.num_minorants_of[fidx] += 1;
        self.minorants[minidx][fidx] = minorant;

            _ => unreachable!("Invalid number of minorants in subproblem {}", fidx),
        }
    }

    fn add_vars<SErr>(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtender<P, SErr>,
        extend_subgradient: &mut SubgradientExtender<Mn, SErr>,
    ) -> Result<(), Either<Self::Err, SErr>> {
        if self.num_subproblems_with_1 == 0 {
        if nnew == 0 || self.num_subproblems_with_1 == 0 {
            return Ok(());
        }

        let noldvars = self.minorants[0][self
            .num_minorants_of
            .iter()
            .position(|&n| n > 0)
            .ok_or(MinimalMasterError::NoMinorants)
            .map_err(Either::Left)?]
        .linear
        .len();
        let mut changedvars = vec![];
        changedvars.extend_from_slice(changed);
        changedvars.extend(noldvars..noldvars + nnew);

        for fidx in 0..self.num_subproblems {
            for i in 0..self.num_minorants_of[fidx] {
                let new_subg =
                    (extend_subgradient)(fidx, &self.minorants[i][fidx].primal, &changedvars).map_err(Either::Right)?;
                let m = &mut self.minorants[i][fidx];
                (extend_subgradient)(fidx, &mut self.minorants[i][fidx]).map_err(Either::Right)?;
                for (&j, &g) in changed.iter().zip(new_subg.iter()) {
                    m.linear[j] = g;
                }
                m.linear.extend_from_slice(&new_subg[changed.len()..]);
            }
        }

        Ok(())
    }

    #[allow(unused_variables)]
    fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<(), Self::Err> {
        for fidx in 0..self.num_subproblems {
            for i in 0..self.num_minorants_of[fidx] {
                debug!("  min(fidx:{}, i:{}) = {}", fidx, i, self.minorants[i][fidx]);
                debug!("  min(fidx:{}, i:{}) = {:?}", fidx, i, self.minorants[i][fidx].debug());
            }
        }

        if self.num_minorants == 2 {
            let min0 = Minorant::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx])));
            let min1 = Minorant::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[1][fidx])));
            let xx = min0.linear.dot(&min0.linear);
            let yy = min1.linear.dot(&min1.linear);
            let xy = min0.linear.dot(&min1.linear);
            let xeta = min0.linear.dot(eta);
            let yeta = min1.linear.dot(eta);
            let min0: Mn = Aggregatable::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx])));
            let min1: Mn = Aggregatable::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[1][fidx])));
            let xx = min0.product(&min0);
            let yy = min1.product(&min1);
            let xy = min0.product(&min1);
            let xeta = min0.linear(eta);
            let yeta = min1.linear(eta);
            let a = yy - 2.0 * xy + xx;
            let b = xy - xx - yeta + xeta;

            let mut alpha2 = 0.0;
            if a > 0.0 {
                alpha2 = ((min1.constant - min0.constant) * self.weight - b) / a;
                alpha2 = ((min1.constant() - min0.constant()) * self.weight - b) / a;
                alpha2 = alpha2.max(0.0).min(1.0);
            }
            self.opt_mult[0] = 1.0 - alpha2;
            self.opt_mult[1] = alpha2;
            Mn::combine_to_vec(
            self.opt_minorant = Aggregatable::combine(self.opt_mult.iter().cloned().zip([min0, min1].iter()));
                self.opt_mult.iter().cloned().zip([min0, min1].iter()),
                &mut self.opt_minorant,
            );
        } else if self.num_minorants == 1 {
            let min0 = Aggregatable::combine((0..self.num_subproblems).map(|fidx| (1.0, &self.minorants[0][fidx])));
            self.opt_minorant = min0;
            let mins0 = &self.minorants[0];
            Mn::combine_to_vec(
                std::iter::repeat(1.0).zip(mins0).take(self.num_subproblems),
                &mut self.opt_minorant,
            );
            self.opt_mult[0] = 1.0;
        } else {
            return Err(MinimalMasterError::NoMinorants);
        }

        debug!("Unrestricted");
        debug!("  opt_minorant={}", self.opt_minorant);
        debug!("  opt_minorant={:?}", self.opt_minorant.debug());
        debug!("  opt_mult={:?}", &self.opt_mult[0..self.num_minorants]);

        Ok(())
    }

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

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

    fn multiplier(&self, min: usize) -> Real {
        self.opt_mult[min % 2]
    }

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

    fn eval_model(&self, y: &DVector) -> Real {
        let mut result = NEG_INFINITY;
        for mins in &self.minorants[0..self.num_minorants] {
            result = result.max(mins.iter().map(|m| m.eval(y)).sum());

impl Default for Builder {
    fn default() -> Self {
        Builder
    }
}

impl<P: Aggregatable + Send + 'static> super::Builder<P> for Builder {
    type MasterProblem = MinimalMaster<P>;
impl<Mn: Minorant> super::Builder<Mn> for Builder {
    type MasterProblem = MinimalMaster<Mn>;

    fn build(&mut self) -> Result<MinimalMaster<P>, MinimalMasterError> {
    fn build(&mut self) -> Result<MinimalMaster<Mn>, MinimalMasterError> {
        MinimalMaster::new()
    }
}

//!   bounds), i.e. replacing $\hat{f}$ by $d \mapsto \hat{f}(d -
//!   \hat{d})$ for some given $\hat{d} \in \mathbb{R}\^n$.

pub mod boxed;
pub use self::boxed::BoxedMasterProblem;

use crate::problem::SubgradientExtender;
use crate::{Aggregatable, DVector, Minorant, Real};
use crate::{DVector, Minorant, Real};
use either::Either;
use std::error::Error;
use std::result::Result;

/// Trait for master problems of a proximal bundle methods.
///
/// Note that solvers are allowed to create multiple master problems potentially
/// running them in different threads. Because of this they must implement `Send
/// + 'static`, i.e. they must be quite self-contained and independent from
/// everything else.
pub trait MasterProblem<P: Aggregatable>: Send + 'static {
pub trait MasterProblem<M: Minorant>: Send + 'static {
    /// Unique index for a minorant.
    type MinorantIndex: Copy + Eq;

    /// The error returned by the master problem.
    type Err: Error + Send + Sync;

    /// Create a new master problem.

    /// Set the maximal number of inner iterations.
    fn set_max_updates(&mut self, max_updates: usize) -> Result<(), Self::Err>;

    /// Return the current number of inner iterations.
    fn cnt_updates(&self) -> usize;

    /// Add or move some variables with bounds.
    /// Add some variables with bounds.
    ///
    /// If an index is specified, existing variables are moved,
    /// otherwise new variables are generated.
    fn add_vars<SErr>(
        &mut self,
        bounds: &[(Option<usize>, Real, Real)],
        extend_subgradient: &mut SubgradientExtender<P, SErr>,
        bounds: &[(Real, Real)],
        extend_subgradient: &mut SubgradientExtender<M, SErr>,
    ) -> Result<(), Either<Self::Err, SErr>>;

    /// Add a new minorant to the model.
    ///
    /// The function returns a unique (among all minorants of all
    /// subproblems) index of the minorant. This index must remain
    /// valid until the minorant is aggregated.
    fn add_minorant(&mut self, fidx: usize, minorant: Minorant<P>) -> Result<Self::MinorantIndex, Self::Err>;
    fn add_minorant(&mut self, fidx: usize, minorant: M) -> Result<Self::MinorantIndex, Self::Err>;

    /// Solve the master problem.
    fn solve(&mut self, cur_value: Real) -> Result<(), Self::Err>;

    /// Compress the bundle.
    ///
    /// When some minorants are compressed, the callback is called with the

    /// $g\^*$ is the optimal aggregated subgradient.
    fn get_dualoptnorm2(&self) -> Real;

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

    /// Return the aggregated primal.
    fn aggregated_primal(&self, fidx: usize) -> Result<P, Self::Err>;
    fn aggregated_primal(&self, fidx: usize) -> Result<M::Primal, Self::Err>;

    /// Move the center of the master problem to $\alpha \cdot d$.
    ///
    /// The vector `d` is allowed to be shorter than the current number of
    /// variables in the master. The missing elements are assumed to be zero.
    fn move_center(&mut self, alpha: Real, d: &DVector);
}

/// A builder for creating a master problem solvers.
pub trait Builder<P: Aggregatable> {
pub trait Builder<M: Minorant> {
    /// The master problem to be build.
    type MasterProblem: MasterProblem<P>;
    type MasterProblem: MasterProblem<M>;

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

/// The full bundle builder.
pub type FullMasterBuilder = boxed::Builder<boxed::unconstrained::cpx::Builder>;

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

        aggr
    }
}

impl ParallelProblem for MMCFProblem {
    type Err = Error;

    type Primal = DVector;
    type Minorant = (Real, DVector, DVector);

    fn num_variables(&self) -> usize {
        self.active_constraints.read().unwrap().len()
    }

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

        let active_constraints = self.active_constraints.read().unwrap()[0..y.len()].to_vec();
        self.pool
            .as_ref()
            .unwrap()
            .execute(move || match sub.write().unwrap().evaluate(&y, active_constraints) {
                Ok((objective, subg, primal)) => {
                    tx.objective(objective).unwrap();
                    tx.minorant(Minorant {
                    tx.minorant((objective, subg, primal)).unwrap();
                        constant: objective,
                        linear: subg,
                        primal,
                    })
                    .unwrap();
                }
                Err(err) => tx.error(err).unwrap(),
            });
        Ok(())
    }

    fn update<U, S>(&mut self, state: U, tx: S) -> Result<()>
    where
        U: ParallelUpdateState<Self::Primal>,
        U: ParallelUpdateState<<Self::Minorant as Minorant>::Primal>,
        S: UpdateSender<Self> + 'static,
    {
        if self.inactive_constraints.read().unwrap().is_empty() {
            return Ok(());
        }

        let inactive_constraints = self.inactive_constraints.clone();

            //     }));

            let nnew = active.len() - nold;
            if nnew > 0 {
                let actives = active.clone();
                if let Err(err) = tx.add_variables(
                    vec![(0.0, Real::infinity()); nnew],
                    Box::new(move |fidx: usize, primal: &DVector, inds: &[usize]| {
                    Box::new(move |fidx: usize, m: &mut Self::Minorant| {
                        let sub = subs[fidx].read().unwrap();
                        Ok(inds
                            .iter()
                            .map(|&i| sub.evaluate_constraint(primal, actives[i]))
                            .collect())
                        let n = m.dim();
                        let primal = &m.2;
                        m.1.extend((n..n + nnew).map(|i| sub.evaluate_constraint(primal, actives[i])));
                        Ok(())
                    }),
                ) {
                    warn!("Error sending problem update: {}", err);
                }
            }
        });

        Ok(())
    }
}

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

//! An asynchronous first-order oracle.

pub use crate::data::minorant::SubgradientExtender;
use crate::{Aggregatable, DVector, Minorant, Real};
use crate::{DVector, Minorant, Real};
use std::sync::Arc;

/// Channel to send evaluation results to.
pub trait ResultSender<P>: Send
where
    P: FirstOrderProblem,
{
    type Err: std::error::Error + Send;

    /// Send a new objective `value`.
    fn objective(&self, value: Real) -> Result<(), Self::Err>;

    /// Send a new `minorant` with associated `primal`.
    fn minorant(&self, minorant: Minorant<P::Primal>) -> Result<(), Self::Err>;
    fn minorant(&self, minorant: P::Minorant) -> Result<(), Self::Err>;

    /// Send an error message.
    fn error(&self, err: P::Err) -> Result<(), Self::Err>;
}

/// Channel to send problem updates to.
pub trait UpdateSender<P>: Send
where
    P: FirstOrderProblem,
{
    type Err: std::error::Error + Send;

    /// Add new variables with given bounds.
    ///
    /// The oracle must provide a "subgradient extension callback", that
    /// computes subgradients for the given variables based on the associated
    /// primals.
    fn add_variables(
        &self,
        bounds: Vec<(Real, Real)>,
        sgext: Box<SubgradientExtender<P::Primal, P::Err>>,
        sgext: Box<SubgradientExtender<P::Minorant, P::Err>>,
    ) -> Result<(), Self::Err>;

    /// Send an error message.
    fn error(&self, err: P::Err) -> Result<(), Self::Err>;
}

/// This trait provides information available in the

///
/// Note that all computations made by an implementation are supposed to be
/// asynchronous.
pub trait FirstOrderProblem {
    /// Error raised by this oracle.
    type Err;

    /// The primal information associated with a minorant.
    type Primal: Aggregatable + Send + 'static;
    /// The minorant information.
    type Minorant: Minorant;

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

    /// Return the lower bounds on the variables.
    ///
    /// If no lower bounds are specified, $-\infty$ is assumed.

    /// information (e.g. adding new variables) to the provided channel.
    ///
    /// The updates might be generated asynchronously.
    ///
    /// The default implementation does nothing.
    fn update<U, S>(&mut self, _state: U, _tx: S) -> Result<(), Self::Err>
    where
        U: UpdateState<Self::Primal>,
        U: UpdateState<<Self::Minorant as Minorant>::Primal>,
        S: UpdateSender<Self> + 'static,
        Self: Sized,
    {
        Ok(())
    }
}

use num_cpus;
use num_traits::Float;
use std::iter::repeat;
use std::sync::Arc;
use std::time::Instant;
use threadpool::ThreadPool;

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

use super::channels::{
    ChannelResultSender, ChannelUpdateSender, ClientReceiver, ClientSender, EvalResult, Message, Update,
};
use super::masterprocess::{self, MasterConfig, MasterProcess, MasterResponse, Response};
use crate::master::{Builder as MasterBuilder, MasterProblem};
use crate::problem::{FirstOrderProblem, UpdateState};

///
/// - `T` is the value type,
/// - `P` is the `FirstOrderProblem` associated with the solver,
/// - `M` is the `MasterBuilder` associated with the solver.
pub type Result<T, P, M> = std::result::Result<
    T,
    Error<
        <<M as MasterBuilder<<P as FirstOrderProblem>::Primal>>::MasterProblem as MasterProblem<
            <P as FirstOrderProblem>::Primal,
        <<M as MasterBuilder<<P as FirstOrderProblem>::Minorant>>::MasterProblem as MasterProblem<
            <P as FirstOrderProblem>::Minorant,
        >>::Err,
        <P as FirstOrderProblem>::Err,
    >,
>;

impl<MErr, PErr> std::fmt::Display for Error<MErr, PErr>
where

}

/// Implementation of a parallel bundle method.
pub struct Solver<P, T = StandardTerminator, W = HKWeighter, M = crate::master::FullMasterBuilder>
where
    P: FirstOrderProblem,
    P::Err: 'static,
    M: MasterBuilder<P::Primal>,
    M: MasterBuilder<P::Minorant>,
{
    /// Parameters for the solver.
    pub params: Parameters,

    /// Termination predicate.
    pub terminator: T,

impl<P, T, W, M> Solver<P, T, W, M>
where
    P: FirstOrderProblem,
    P::Err: Send + 'static,
    T: Terminator<SolverData> + Default,
    W: Weighter<SolverData> + Default,
    M: MasterBuilder<P::Primal>,
    M: MasterBuilder<P::Minorant>,
{
    /// Create a new parallel bundle solver.
    pub fn new(problem: P) -> Self
    where
        M: Default,
    {
        Self::with_master(problem, M::default())

    /// variables), and the master problem is started to compute a new
    /// candidate.
    ///
    /// The function returns
    ///   - `Ok(true)` if the final iteration count has been reached,
    ///   - `Ok(false)` if the final iteration count has not been reached,
    ///   - `Err(_)` on error.
    fn handle_client_response(&mut self, msg: EvalResult<usize, P::Primal, P::Err>) -> Result<bool, P, M> {
    fn handle_client_response(&mut self, msg: EvalResult<usize, P::Minorant, P::Err>) -> Result<bool, P, M> {
        let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?;
        match msg {
            EvalResult::ObjectiveValue { index, value } => {
                debug!("Receive objective from subproblem {}: {}", index, value);
                if self.data.nxt_ubs[index].is_infinite() {
                    self.data.cnt_remaining_ubs -= 1;
                }
                self.data.nxt_ubs[index] = self.data.nxt_ubs[index].min(value);
            }
            EvalResult::Minorant { index, mut minorant } => {
                debug!("Receive minorant from subproblem {}", index);
                if self.data.nxt_cutvals[index].is_infinite() {
                    self.data.cnt_remaining_mins -= 1;
                }
                // move center of minorant to cur_y
                minorant.move_center(-1.0, &self.data.nxt_d);
                self.data.nxt_cutvals[index] = self.data.nxt_cutvals[index].max(minorant.constant);
                self.data.nxt_cutvals[index] = self.data.nxt_cutvals[index].max(minorant.constant());
                // add minorant to master problem
                master.add_minorant(index, minorant)?;
            }
            EvalResult::Done { .. } => return Ok(false), // currently we do nothing ...
            EvalResult::Error { err, .. } => return Err(Error::Evaluation(err)),
        }

    }

    /// Handles an update response `update` of the problem.
    ///
    /// The method is called if the problem informs the solver about a change of
    /// the problem, e.g. adding a new variable. This method updates the other
    /// parts of the solver, e.g. the master problem, about the modification.
    fn handle_update_response(&mut self, update: Update<usize, P::Primal, P::Err>) -> Result<(), P, M> {
    fn handle_update_response(&mut self, update: Update<usize, P::Minorant, P::Err>) -> Result<(), P, M> {
        match update {
            Update::AddVariables { bounds, sgext, .. } => {
                if !bounds.is_empty() {
                    // add new variables
                    self.master_proc
                        .as_mut()
                        .unwrap()
                        .add_vars(bounds.into_iter().map(|(l, u)| (None, l, u)).collect(), sgext)?;
                        .add_vars(bounds.into_iter().map(|(l, u)| (l, u)).collect(), sgext)?;
                }
            }
            Update::Done { .. } => self.data.update_in_progress = false, // currently we do nothing ...
            Update::Error { err, .. } => {
                self.data.update_in_progress = false;
                return Err(Error::Update(err));
            }

            self.data.nxt_mod,
            self.data.nxt_val,
            self.data.cur_val
        );
    }

    /// Return the aggregated primal of the given subproblem.
    pub fn aggregated_primal(&self, subproblem: usize) -> Result<P::Primal, P, M> {
    pub fn aggregated_primal(&self, subproblem: usize) -> Result<<P::Minorant as Minorant>::Primal, P, M> {
        Ok(self
            .master_proc
            .as_ref()
            .ok_or(Error::NotSolved)?
            .get_aggregated_primal(subproblem)?)
    }
}

 */

//! Implementation for `ResultSender` using channels.

use super::masterprocess::Response as MasterResponse;
use crate::master::MasterProblem;
use crate::problem::{FirstOrderProblem, ResultSender, SubgradientExtender, UpdateSender};
use crate::{Aggregatable, Minorant, Real};
use crate::{Minorant, Real};

#[cfg(feature = "crossbeam")]
use rs_crossbeam::channel::{Receiver, SendError, Sender};
#[cfg(not(feature = "crossbeam"))]
use std::sync::mpsc::{Receiver, SendError, Sender};

/// A message received from a client.
///
/// The client might be either a subproblem or a master problem process.
///
/// A client may send one of two possible messages: either an evaluation result
/// or a problem update.
pub enum Message<I, P, E, MErr>
pub enum Message<I, Mn, E, MErr>
where
    P: Aggregatable,
    Mn: Minorant,
{
    /// An evaluation result.
    Eval(EvalResult<I, P, E>),
    Eval(EvalResult<I, Mn, E>),
    /// A problem update.
    Update(Update<I, P, E>),
    Update(Update<I, Mn, E>),
    /// A master problem result.
    Master(MasterResponse<MErr, E>),
}

/// The sender for client messages.
pub type ClientSender<P, M> = Sender<
    Message<
        usize,
        <P as FirstOrderProblem>::Primal,
        <P as FirstOrderProblem>::Minorant,
        <P as FirstOrderProblem>::Err,
        <M as MasterProblem<<P as FirstOrderProblem>::Primal>>::Err,
        <M as MasterProblem<<P as FirstOrderProblem>::Minorant>>::Err,
    >,
>;

/// The receiver for client messages.
pub type ClientReceiver<P, M> = Receiver<
    Message<
        usize,
        <P as FirstOrderProblem>::Primal,
        <P as FirstOrderProblem>::Minorant,
        <P as FirstOrderProblem>::Err,
        <M as MasterProblem<<P as FirstOrderProblem>::Primal>>::Err,
        <M as MasterProblem<<P as FirstOrderProblem>::Minorant>>::Err,
    >,
>;

/// Parameters for tuning the solver.

/// Evaluation result.
///
/// The result of an evaluation is new information to be made
/// available to the solver and the master problem. There are
/// essentially two types of information:
///
///    1. The (exact) function value of a sub-function at some point.
///    2. A minorant of some sub-function.
#[derive(Debug)]
pub enum EvalResult<I, P, E>
pub enum EvalResult<I, Mn, E>
where
    P: Aggregatable,
    Mn: Minorant,
{
    /// The objective value at some point.
    ObjectiveValue { index: I, value: Real },
    /// A minorant with an associated primal.
    Minorant { index: I, minorant: Minorant<P> },
    Minorant { index: I, minorant: Mn },
    /// An error occurred.
    Error { index: I, err: E },
    /// Done, no further message will be received from this evaluation.
    Done { index: I },
}

/// Problem update information.
///
/// The solver calls the `update` method of the problem regularly.
/// This method can modify the problem by adding (or moving)
/// variables. The possible updates are encoded in this type.
pub enum Update<I, P, E> {
pub enum Update<I, Mn, E>
where
    Mn: Minorant,
{
    /// Add new variables with bounds.
    ///
    /// The initial value of each variable will be the feasible value
    /// closest to 0.
    AddVariables {
        index: I,
        bounds: Vec<(Real, Real)>,
        sgext: Box<SubgradientExtender<P, E>>,
        sgext: Box<SubgradientExtender<Mn, E>>,
    },
    /// An error occurred.
    Error { index: I, err: E },
    /// Done, no further message will be received from this evaluation.
    Done { index: I },
}

pub struct ChannelResultSender<I, P, E, M>
pub struct ChannelResultSender<I, Mn, E, M>
where
    I: Clone,
    P: Aggregatable,
    Mn: Minorant,
{
    index: I,
    tx: Sender<Message<I, P, E, M>>,
    tx: Sender<Message<I, Mn, E, M>>,
}

impl<I, P, E, M> ChannelResultSender<I, P, E, M>
impl<I, Mn, E, M> ChannelResultSender<I, Mn, E, M>
where
    I: Clone,
    P: Aggregatable,
    Mn: Minorant,
{
    pub fn new(index: I, tx: Sender<Message<I, P, E, M>>) -> Self {
    pub fn new(index: I, tx: Sender<Message<I, Mn, E, M>>) -> Self {
        ChannelResultSender { index, tx }
    }
}

impl<I, P, M> ResultSender<P> for ChannelResultSender<I, P::Primal, P::Err, M>
impl<I, P, M> ResultSender<P> for ChannelResultSender<I, P::Minorant, P::Err, M>
where
    I: Clone + Send,
    P: FirstOrderProblem,
    P::Err: Send,
    M: Send,
{
    type Err = SendError<Message<I, P::Primal, P::Err, M>>;
    type Err = SendError<Message<I, P::Minorant, P::Err, M>>;

    fn objective(&self, value: Real) -> Result<(), Self::Err> {
        self.tx.send(Message::Eval(EvalResult::ObjectiveValue {
            index: self.index.clone(),
            value,
        }))
    }

    fn minorant(&self, minorant: Minorant<P::Primal>) -> Result<(), Self::Err> {
    fn minorant(&self, minorant: P::Minorant) -> Result<(), Self::Err> {
        self.tx.send(Message::Eval(EvalResult::Minorant {
            index: self.index.clone(),
            minorant,
        }))
    }

    fn error(&self, err: P::Err) -> Result<(), Self::Err> {
        self.tx.send(Message::Eval(EvalResult::Error {
            index: self.index.clone(),
            err,
        }))
    }
}

impl<I, P, E, M> Drop for ChannelResultSender<I, P, E, M>
impl<I, Mn, E, M> Drop for ChannelResultSender<I, Mn, E, M>
where
    I: Clone,
    P: Aggregatable,
    Mn: Minorant,
{
    fn drop(&mut self) {
        self.tx
            .send(Message::Eval(EvalResult::Done {
                index: self.index.clone(),
            }))
            .unwrap()
    }
}

pub struct ChannelUpdateSender<I, P, E, M>
pub struct ChannelUpdateSender<I, Mn, E, M>
where
    I: Clone,
    P: Aggregatable,
    Mn: Minorant,
{
    index: I,
    tx: Sender<Message<I, P, E, M>>,
    tx: Sender<Message<I, Mn, E, M>>,
}

impl<I, P, E, M> ChannelUpdateSender<I, P, E, M>
impl<I, Mn, E, M> ChannelUpdateSender<I, Mn, E, M>
where
    I: Clone,
    P: Aggregatable,
    Mn: Minorant,
{
    pub fn new(index: I, tx: Sender<Message<I, P, E, M>>) -> Self {
    pub fn new(index: I, tx: Sender<Message<I, Mn, E, M>>) -> Self {
        ChannelUpdateSender { index, tx }
    }
}

impl<I, P, M> UpdateSender<P> for ChannelUpdateSender<I, P::Primal, P::Err, M>
impl<I, P, M> UpdateSender<P> for ChannelUpdateSender<I, P::Minorant, P::Err, M>
where
    I: Clone + Send,
    P: FirstOrderProblem,
    P::Err: Send,
    M: Send,
{
    type Err = SendError<Message<I, P::Primal, P::Err, M>>;
    type Err = SendError<Message<I, P::Minorant, P::Err, M>>;

    fn add_variables(
        &self,
        bounds: Vec<(Real, Real)>,
        sgext: Box<SubgradientExtender<P::Primal, P::Err>>,
        sgext: Box<SubgradientExtender<P::Minorant, P::Err>>,
    ) -> Result<(), Self::Err> {
        self.tx.send(Message::Update(Update::AddVariables {
            index: self.index.clone(),
            bounds,
            sgext,
        }))
    }

    fn error(&self, err: P::Err) -> Result<(), Self::Err> {
        self.tx.send(Message::Update(Update::Error {
            index: self.index.clone(),
            err,
        }))
    }
}

impl<I, P, E, M> Drop for ChannelUpdateSender<I, P, E, M>
impl<I, Mn, E, M> Drop for ChannelUpdateSender<I, Mn, E, M>
where
    I: Clone,
    P: Aggregatable,
    Mn: Minorant,
{
    fn drop(&mut self) {
        self.tx
            .send(Message::Update(Update::Done {
                index: self.index.clone(),
            }))
            .unwrap()
    }
}

use log::{debug, warn};
use std::sync::Arc;
use threadpool::ThreadPool;

use super::channels::{ClientSender, Message};
use crate::master::MasterProblem;
use crate::problem::{FirstOrderProblem, SubgradientExtender};
use crate::{Aggregatable, DVector, Minorant, Real};
use crate::{DVector, Minorant, Real};

#[derive(Debug)]
pub enum Error<MErr, PErr> {
    /// The communication channel for sending requests has been disconnected.
    DisconnectedSender,
    /// The communication channel for sending responds has been disconnected.
    DisconnectedReceiver,

    /// The lower bounds on the variables.
    pub lower_bounds: Option<DVector>,
    /// The lower bounds on the variables.
    pub upper_bounds: Option<DVector>,
}

/// A task for the master problem.
enum MasterTask<Pr, PErr, M>
enum MasterTask<Mn, PErr, M>
where
    Mn: Minorant,
    M: MasterProblem<Pr>,
    M: MasterProblem<Mn>,
    Pr: Aggregatable,
{
    /// Add new variables to the master problem.
    AddVariables(Vec<(Option<usize>, Real, Real)>, Box<SubgradientExtender<Pr, PErr>>),
    /// Add new variables with bounds to the master problem.
    AddVariables(Vec<(Real, Real)>, Box<SubgradientExtender<Mn, PErr>>),

    /// Add a new minorant for a subfunction to the master problem.
    AddMinorant(usize, Minorant<Pr>),
    AddMinorant(usize, Mn),

    /// Move the center of the master problem in the given direction.
    MoveCenter(Real, Arc<DVector>),

    /// Start a new computation of the master problem.
    Solve { center_value: Real },

    /// Compress the bundle.
    Compress,

    /// Set the weight parameter of the master problem.
    SetWeight { weight: Real },

    /// Return the current aggregated primal.
    GetAggregatedPrimal {
        subproblem: usize,
        tx: Sender<Result<Pr, M::Err>>,
        tx: Sender<Result<Mn::Primal, M::Err>>,
    },
}

/// The response send from a master process.
///
/// The response contains the evaluation results of the latest call to `solve`.
pub enum Response<MErr, PErr> {

    },
    /// An error occurred.
    Error(Error<MErr, PErr>),
}

/// Convenient type for `Response`.
pub type MasterResponse<P, M> =
    Response<<M as MasterProblem<<P as FirstOrderProblem>::Primal>>::Err, <P as FirstOrderProblem>::Err>;
    Response<<M as MasterProblem<<P as FirstOrderProblem>::Minorant>>::Err, <P as FirstOrderProblem>::Err>;

type ToMasterSender<P, M> = Sender<MasterTask<<P as FirstOrderProblem>::Primal, <P as FirstOrderProblem>::Err, M>>;
type ToMasterSender<P, M> = Sender<MasterTask<<P as FirstOrderProblem>::Minorant, <P as FirstOrderProblem>::Err, M>>;

type ToMasterReceiver<P, M> = Receiver<MasterTask<<P as FirstOrderProblem>::Primal, <P as FirstOrderProblem>::Err, M>>;
type ToMasterReceiver<P, M> =
    Receiver<MasterTask<<P as FirstOrderProblem>::Minorant, <P as FirstOrderProblem>::Err, M>>;

pub struct MasterProcess<P, M>
where
    P: FirstOrderProblem,
    M: MasterProblem<P::Primal>,
    M: MasterProblem<P::Minorant>,
{
    /// The channel to transmit new tasks to the master problem.
    tx: ToMasterSender<P, M>,

    phantom: std::marker::PhantomData<M>,
}

impl<P, M> MasterProcess<P, M>
where
    P: FirstOrderProblem,
    P::Primal: Send + 'static,
    P::Minorant: Send + 'static,
    P::Err: Send + 'static,
    M: MasterProblem<P::Primal> + Send + 'static,
    M: MasterProblem<P::Minorant> + Send + 'static,
    M::Err: Send + 'static,
{
    pub fn start(master: M, master_config: MasterConfig, tx: ClientSender<P, M>, threadpool: &mut ThreadPool) -> Self {
        // Create a pair of communication channels.
        let (to_master_tx, to_master_rx) = channel();

        // The the master process thread.

            phantom: std::marker::PhantomData,
        }
    }

    /// Add new variables to the master problem.
    pub fn add_vars(
        &mut self,
        vars: Vec<(Option<usize>, Real, Real)>,
        sgext: Box<SubgradientExtender<P::Primal, P::Err>>,
        vars: Vec<(Real, Real)>,
        sgext: Box<SubgradientExtender<P::Minorant, P::Err>>,
    ) -> Result<(), Error<M::Err, P::Err>>
    where
        P::Err: 'static,
    {
        Ok(self.tx.send(MasterTask::AddVariables(vars, sgext))?)
    }

    /// Add a new minorant to the master problem model.
    ///
    /// This adds the specified `minorant` with associated `primal` data to the
    /// model of subproblem `i`.
    pub fn add_minorant(&mut self, i: usize, minorant: Minorant<P::Primal>) -> Result<(), Error<M::Err, P::Err>> {
    pub fn add_minorant(&mut self, i: usize, minorant: P::Minorant) -> Result<(), Error<M::Err, P::Err>> {
        Ok(self.tx.send(MasterTask::AddMinorant(i, minorant))?)
    }

    /// Move the center of the master problem.
    ///
    /// This moves the master problem's center in direction $\\alpha \\cdot d$.
    pub fn move_center(&mut self, alpha: Real, d: Arc<DVector>) -> Result<(), Error<M::Err, P::Err>> {

    /// Sets the new weight of the proximal term in the master problem.
    pub fn set_weight(&mut self, weight: Real) -> Result<(), Error<M::Err, P::Err>> {
        Ok(self.tx.send(MasterTask::SetWeight { weight })?)
    }

    /// Get the current aggregated primal for a certain subproblem.
    pub fn get_aggregated_primal(&self, subproblem: usize) -> Result<P::Primal, Error<M::Err, P::Err>> {
    pub fn get_aggregated_primal(
        &self,
        subproblem: usize,
    ) -> Result<<P::Minorant as Minorant>::Primal, Error<M::Err, P::Err>> {
        let (tx, rx) = channel();
        self.tx.send(MasterTask::GetAggregatedPrimal { subproblem, tx })?;
        rx.recv()?.map_err(Error::Aggregation)
    }

    /// The main loop of the master process.
    fn master_main(
        master: M,
        master_config: MasterConfig,
        tx: &mut ClientSender<P, M>,
        rx: ToMasterReceiver<P, M>,
        subgradient_extender: &mut Option<Box<SubgradientExtender<P::Primal, P::Err>>>,
        subgradient_extender: &mut Option<Box<SubgradientExtender<P::Minorant, P::Err>>>,
    ) -> Result<(), Error<M::Err, P::Err>> {
        let mut master = master;

        // Initialize the master problem.
        master
            .set_num_subproblems(master_config.num_subproblems)
            .map_err(Error::Master)?;

                MasterTask::AddMinorant(i, mut m) => {
                    debug!("master: add minorant to subproblem {}", i);
                    // It may happen the number new minorant belongs to an earlier evaluation
                    // with less variables (i.e. new variables have been added
                    // after the start of the evaluation but before the new
                    // minorant is added, i.e. now). In this case we must add
                    // extend the minorant accordingly.
                    if m.linear.len() < master.num_variables() {
                    if m.dim() < master.num_variables() {
                        if let Some(ref mut sgext) = subgradient_extender {
                            let newinds = (m.linear.len()..master.num_variables()).collect::<Vec<_>>();
                            let new_subg =
                                sgext(i, &m.primal, &newinds).map_err(Error::<M::Err, P::Err>::SubgradientExtension)?;
                            sgext(i, &mut m).map_err(Error::<M::Err, P::Err>::SubgradientExtension)?;
                            m.linear.extend(new_subg);
                        }
                    }
                    master.add_minorant(i, m).map_err(Error::Master)?;
                }
                MasterTask::MoveCenter(alpha, d) => {
                    debug!("master: move center");
                    master.move_center(alpha, &d);

#[cfg(feature = "crossbeam")]
use rs_crossbeam::channel::{unbounded as channel, RecvError};
#[cfg(not(feature = "crossbeam"))]
use std::sync::mpsc::{channel, RecvError};

use log::{debug, info, warn};
use num_cpus;
use num_traits::Float;
use num_traits::{Float, Zero};
use std::sync::Arc;
use std::time::Instant;
use threadpool::ThreadPool;

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

use super::channels::{
    ChannelResultSender, ChannelUpdateSender, ClientReceiver, ClientSender, EvalResult, Message, Update,
};
use super::masterprocess::{self, MasterConfig, MasterProcess, MasterResponse, Response};
use crate::master::{Builder as MasterBuilder, MasterProblem};
use crate::problem::{FirstOrderProblem, UpdateState};

///
/// - `T` is the value type,
/// - `P` is the `FirstOrderProblem` associated with the solver,
/// - `M` is the `MasterBuilder` associated with the solver.
pub type Result<T, P, M> = std::result::Result<
    T,
    Error<
        <<M as MasterBuilder<<P as FirstOrderProblem>::Primal>>::MasterProblem as MasterProblem<
            <P as FirstOrderProblem>::Primal,
        <<M as MasterBuilder<<P as FirstOrderProblem>::Minorant>>::MasterProblem as MasterProblem<
            <P as FirstOrderProblem>::Minorant,
        >>::Err,
        <P as FirstOrderProblem>::Err,
    >,
>;

impl<MErr, PErr> std::fmt::Display for Error<MErr, PErr>
where

    }
}

/// Implementation of a parallel bundle method.
pub struct Solver<P, T = StandardTerminator, W = HKWeighter, M = crate::master::FullMasterBuilder>
where
    P: FirstOrderProblem,
    M: MasterBuilder<P::Primal>,
    M: MasterBuilder<P::Minorant>,
    P::Err: 'static,
{
    /// Parameters for the solver.
    pub params: Parameters,

    /// Termination predicate.
    pub terminator: T,

impl<P, T, W, M> Solver<P, T, W, M>
where
    P: FirstOrderProblem,
    P::Err: Send + 'static,
    T: Terminator<SolverData> + Default,
    W: Weighter<SolverData> + Default,
    M: MasterBuilder<P::Primal>,
    M: MasterBuilder<P::Minorant>,
{
    /// Create a new parallel bundle solver.
    pub fn new(problem: P) -> Self
    where
        M: Default,
    {
        Self::with_master(problem, M::default())

            }
        }
    }

    /// Handle a response from a subproblem evaluation.
    ///
    /// The function returns `Ok(true)` if the final iteration count has been reached.
    fn handle_client_response(&mut self, msg: EvalResult<usize, P::Primal, P::Err>) -> Result<bool, P, M> {
    fn handle_client_response(&mut self, msg: EvalResult<usize, P::Minorant, P::Err>) -> Result<bool, P, M> {
        let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?;
        match msg {
            EvalResult::ObjectiveValue { index, value } => {
                debug!("Receive objective from subproblem {}: {}", index, value);
                if self.data.nxt_ubs[index].is_infinite() {
                    self.data.cnt_remaining_ubs -= 1;
                }
                self.data.nxt_ubs[index] = self.data.nxt_ubs[index].min(value);
            }
            EvalResult::Minorant { index, mut minorant } => {
                debug!("Receive minorant from subproblem {}", index);
                if self.data.nxt_cutvals[index].is_infinite() {
                    self.data.cnt_remaining_mins -= 1;
                }
                // move center of minorant to cur_y
                minorant.move_center(-1.0, &self.data.nxt_d);
                self.data.nxt_cutvals[index] = self.data.nxt_cutvals[index].max(minorant.constant);
                self.data.nxt_cutvals[index] = self.data.nxt_cutvals[index].max(minorant.constant());
                // add minorant to master problem
                master.add_minorant(index, minorant)?;
            }
            EvalResult::Done { .. } => return Ok(false), // nothing to do here
            EvalResult::Error { err, .. } => return Err(Error::Evaluation(err)),
        }

        let mut have_update = false;
        for msg in update_rx {
            if let Message::Update(update) = msg {
                match update {
                    Update::AddVariables { bounds, sgext, .. } => {
                        have_update = true;
                        let mut newvars = Vec::with_capacity(bounds.len());
                        for (lower, upper) in bounds {
                            if lower > upper {
                        let newvars = bounds
                            .into_iter()
                            .map(|(lower, upper)| {
                                if lower <= upper {
                                return Err(Error::InvalidBounds { lower, upper });
                            }
                            let value = if lower > 0.0 {
                                lower
                                    let value = lower.max(Real::zero()).min(upper);
                                    Ok((lower - value, upper - value, value))
                            } else if upper < 0.0 {
                                upper
                            } else {
                                } else {
                                0.0
                            };
                            //self.bounds.push((lower, upper));
                                    Err(Error::InvalidBounds { lower, upper })
                            newvars.push((None, lower - value, upper - value, value));
                        }
                                }
                        if !newvars.is_empty() {
                            // modify moved variables
                            for (index, val) in newvars.iter().filter_map(|v| v.0.map(|i| (i, v.3))) {
                                self.data.cur_y[index] = val;
                            }

                            })
                            .collect::<Result<Vec<_>, P, M>>()?;
                        if !newvars.is_empty() {
                            // add new variables
                            self.data
                            self.data.cur_y.extend(newvars.iter().map(|v| v.2));
                                .cur_y
                                .extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3));

                            master_proc.add_vars(newvars.iter().map(|v| (v.0, v.1, v.2)).collect(), sgext)?;
                            master_proc.add_vars(newvars.iter().map(|v| (v.0, v.1)).collect(), sgext)?;
                        }
                    }
                    Update::Done { .. } => (), // there's nothing to do
                    Update::Error { err, .. } => return Err(Error::Update(err)),
                }
            } else {
                unreachable!("Only update results allowed during update");

            self.data.nxt_mod,
            self.data.nxt_val,
            self.data.cur_val
        );
    }

    /// Return the aggregated primal of the given subproblem.
    pub fn aggregated_primal(&self, subproblem: usize) -> Result<P::Primal, P, M> {
    pub fn aggregated_primal(&self, subproblem: usize) -> Result<<P::Minorant as Minorant>::Primal, P, M> {
        Ok(self
            .master_proc
            .as_ref()
            .ok_or(Error::NotSolved)?
            .get_aggregated_primal(subproblem)?)
    }
}