RsBundle  Check-in [1adcbb8b61]

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

//! > Helmberg, C. and Kiwiel, K.C. (2002): A spectral bundle method
//! > with bounds, Math. Programming A 93, 173--194
//!

use Real;
use {BundleState, SolverParams, Step, Weighter};

use std::f64::NEG_INFINITY;
use std::cmp::{max, min};


const FACTOR: Real = 2.0;

/**
 * Weight updating rule according to Helmberg and Kiwiel.
 *
 * The procedure is described in

    }

    /// Create new HKWeighter with weight $m_R$.
    pub fn new_weight(m_r: Real) -> HKWeighter {
        assert!(m_r > 0.0);
        HKWeighter {
            eps_weight: 1e30,
            m_r: m_r,
            iter: 0,
            model_max: NEG_INFINITY,
        }
    }
}

impl Default for HKWeighter {

        let w = 2.0 * state.weight * (1.0 - cur_nxt / cur_mod);

        debug!("  cur_nxt={} cur_mod={} w={}", cur_nxt, cur_mod, w);

        if state.step == Step::Null {
            let sgnorm = state.sgnorm;
            let lin_err = state.cur_val - state.new_cutval;
            self.eps_weight = self.eps_weight
                .min(sgnorm + cur_mod - sgnorm * sgnorm / state.weight);
            let new_weight = if self.iter < -3 && lin_err > self.eps_weight.max(FACTOR * cur_mod) {
                w
            } else {
                state.weight
            }.min(FACTOR * state.weight)
                .min(params.max_weight);
            if new_weight > state.weight {

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

use {DVector, Minorant, Real};
use master::MasterProblem;
use master::UnconstrainedMasterProblem;


use std::result::Result;
use std::f64::{EPSILON, INFINITY, NEG_INFINITY};
use itertools::multizip;

use failure::Error;

/**
 * Turn unconstrained master problem into box-constrained one.
 *
 * This master problem adds box constraints to an unconstrainted
 * master problem implementation. The box constraints are enforced by
 * an additional outer optimization loop.

            primopt: dvec![],
            primoptval: 0.0,
            dualoptnorm2: 0.0,
            model_eps: 0.6,
            max_updates: 100,
            cnt_updates: 0,
            need_new_candidate: true,
            master: master,
        }
    }

    pub fn set_max_updates(&mut self, max_updates: usize) -> Result<(), Error> {
        assert!(max_updates > 0);
        self.max_updates = max_updates;
        Ok(())

    // defined by a fixed $\bar{g}$ while choosing the best possible
    // $\eta$.
    //
    fn compute_candidate(&mut self) {
        self.need_new_candidate = false;

        if self.master.dualopt().len() == self.lb.len() {
            self.primopt
                .scal(-1.0 / self.master.weight(), self.master.dualopt())
        } else {
            self.primopt.init0(self.lb.len());
        }
        self.update_box_multipliers();
    }

    /// Compute $\langle b, \eta \rangle$ with $b$ the bounds of eta.

     * Return $\\|G \alpha - \eta\\|_2\^2$.
     *
     * This is the norm-square of the dual optimal solution including
     * the current box-multipliers $\eta$.
     */
    fn get_norm_subg2(&self) -> Real {
        let dualopt = self.master.dualopt();
        dualopt
            .iter()
            .zip(self.eta.iter())
            .map(|(x, y)| x * y)
            .sum()
    }
}

impl<M: UnconstrainedMasterProblem> MasterProblem for BoxedMasterProblem<M> {
    type MinorantIndex = M::MinorantIndex;

    fn set_num_subproblems(&mut self, n: usize) -> Result<(), Error> {

        extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector,
    ) -> Result<(), Error> {
        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.eta.resize(self.lb.len(), 0.0);
            self.need_new_candidate = true;
            let nnew = bounds.iter().filter(|v| v.0.is_none()).count();
            let changed = bounds.iter().filter_map(|v| v.0).collect::<Vec<_>>();
            self.master.add_vars(nnew, &changed, extend_subgradient)
        } else {
            Ok(())

        let mut cnt_updates = 0;
        let mut old_augval = NEG_INFINITY;
        loop {
            cnt_updates += 1;
            self.cnt_updates += 1;

            // TODO: relprec is fixed
            self.master
                .solve(&self.eta, center_value, old_augval, 1e-3)?;

            // compute the primal solution without the influence of eta
            self.primopt
                .scal(-1.0 / self.master.weight(), self.master.dualopt());

            // solve w.r.t. eta
            let updated_eta = self.update_box_multipliers();

            // compute value of the linearized model
            self.dualoptnorm2 = self.get_norm_subg2();
            let linval = self.master.dualopt().dot(&self.primopt) + self.master.dualopt_cutval();

            debug!("  center_value={}", center_value);
            debug!("  model_eps={}", self.model_eps);
            debug!(
                "  cut-lin={} < eps*(cur-lin)={}",
                cutval - linval,
                self.model_eps * (curval - linval)
            );
            debug!(
                "  cnt_update={} max_updates={}",
                cnt_updates, self.max_updates
            );

            self.primoptval = linval;


            if augval < old_augval + 1e-10 || cutval - linval < self.model_eps * (curval - linval)
                || cnt_updates >= self.max_updates
            {
                break;
            }

            old_augval = old_augval.max(augval);
        }

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

use {DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation};
use mcf;



use std::fs::File;
use std::io::Read;
use std::f64::INFINITY;
use std::result::Result;

use failure::Error;

/// A solver error.
#[derive(Debug, Fail)]
#[fail(display = "Format error: {}", msg)]

        let fnod = buffer
            .split_whitespace()
            .map(|x| x.parse::<usize>().unwrap())
            .collect::<Vec<_>>();

        if fnod.len() != 4 {
            return Err(MCFFormatError {
                msg: format!(
                    "Expected 4 numbers in {}.nod, but got {}",
                    basename,
                    fnod.len()
                ),
            }.into());
        }

        let ncom = fnod[0];
        let nnodes = fnod[1];
        let narcs = fnod[2];
        let ncaps = fnod[3];

            let snk = try!(data.next().unwrap().parse::<usize>()) - 1;
            let com = try!(data.next().unwrap().parse::<usize>()) - 1;
            let cost = try!(data.next().unwrap().parse::<Real>());
            let cap = try!(data.next().unwrap().parse::<Real>());
            let mt = try!(data.next().unwrap().parse::<isize>()) - 1;
            assert!(
                arc < narcs,
                format!(
                    "Wrong arc number (got: {}, expected in 1..{})",
                    arc + 1,
                    narcs
                )
            );
            // set internal coeff
            let coeff = arcmap[com].len();
            arcmap[com].push(ArcInfo {
                arc: arc + 1,
                src: src + 1,
                snk: snk + 1,

            rhs[mt] = cap;
        }

        // set lhs
        let mut lhs = vec![vec![vec![]; ncom]; ncaps];
        for i in 0..ncaps {
            for fidx in 0..ncom {
                lhs[i][fidx] = lhsidx[i][fidx]
                    .iter()
                    .map(|&j| Elem { ind: j, val: 1.0 })
                    .collect();
            }
        }

        Ok(MMCFProblem {
            multimodel: false,
            nets: nets,
            lhs: lhs,
            rhs: rhs,
            rhsval: 0.0,
            cbase: cbase,
            c: vec![dvec![]; ncom],
        })
    }

    /// Compute costs for a primal solution.
    pub fn get_primal_costs(&self, fidx: usize, primals: &[DVector]) -> Real {
        if self.multimodel {

            } else {
                subg = dvec![0.0; self.rhs.len()];
                objective = -try!(self.nets[fidx].objective());
            }

            let sol = try!(self.nets[fidx].get_solution());
            for (i, lhs) in self.lhs.iter().enumerate() {
                subg[i] -= lhs[fidx]
                    .iter()
                    .map(|elem| elem.val * sol[elem.ind])
                    .sum::<Real>();
            }

            Ok(SimpleEvaluation {
                objective: objective,
                minorants: vec![
                    (
                        Minorant {
                            constant: objective,
                            linear: subg,
                        },
                        vec![sol],
                    ),
                ],
            })
        } else {
            let mut objective = self.rhsval;
            let mut sols = Vec::with_capacity(self.nets.len());
            for i in 0..self.nets.len() {
                objective -= try!(self.nets[i].objective());
                sols.push(try!(self.nets[i].get_solution()));
            }

            let mut subg = self.rhs.clone();
            for (i, lhs) in self.lhs.iter().enumerate() {
                for (fidx, flhs) in lhs.iter().enumerate() {
                    subg[i] -= flhs.iter()
                        .map(|elem| elem.val * sols[fidx][elem.ind])
                        .sum::<Real>();
                }
            }

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

    fn aggregate_primals(&mut self, primals: Vec<(Real, Vec<DVector>)>) -> Vec<DVector> {
        self.aggregate_primals_ref(&primals
            .iter()
            .map(|&(alpha, ref x)| (alpha, x))
            .collect::<Vec<_>>())
    }
}

#![allow(unused_unsafe)]

use {DVector, Real};

use cplex_sys as cpx;

use std;
use std::ptr;
use std::ffi::CString;
use std::result::Result;

use std::os::raw::{c_char, c_double, c_int};

use failure::Error;

pub struct Solver {

    pub fn new(nnodes: usize) -> Result<Solver, Error> {
        let mut status: c_int;
        let mut net = ptr::null_mut();

        unsafe {
            #[cfg_attr(feature = "cargo-clippy", allow(never_loop))]
            loop {
                status = cpx::setlogfilename(
                    cpx::env(),
                    c_str!("mcf.cpxlog").as_ptr(),
                    c_str!("w").as_ptr(),
                );
                if status != 0 {
                    break;
                }

                net = cpx::NETcreateprob(cpx::env(), &mut status, c_str!("mcf").as_ptr());
                if status != 0 {
                    break;
                }
                status = cpx::NETaddnodes(cpx::env(), net, nnodes as c_int, ptr::null(), ptr::null());
                if status != 0 {
                    break;
                }
                status = cpx::NETchgobjsen(cpx::env(), net, cpx::ObjectiveSense::Minimize.to_c());
                if status != 0 {
                    break;
                }
                break;
            }

            if status != 0 {
                let msg = CString::new(vec![0; cpx::MESSAGE_BUF_SIZE])
                    .unwrap()
                    .into_raw();
                cpx::geterrorstring(cpx::env(), status, msg);
                cpx::NETfreeprob(cpx::env(), &mut net);
                return Err(cpx::CplexError {
                    code: status,
                    msg: CString::from_raw(msg).to_string_lossy().into_owned(),
                }.into());
            }
        }

        Ok(Solver { net: net })
    }

    pub fn num_nodes(&self) -> usize {
        unsafe { cpx::NETgetnumnodes(cpx::env(), self.net) as usize }
    }

    pub fn num_arcs(&self) -> usize {
        unsafe { cpx::NETgetnumarcs(cpx::env(), self.net) as usize }
    }

    pub fn set_balance(&mut self, node: usize, supply: Real) -> Result<(), Error> {
        let n = node as c_int;
        let s = supply as c_double;
        Ok(trycpx!(cpx::NETchgsupply(
            cpx::env(),
            self.net,
            1,
            &n,
            &s as *const c_double
        )))
    }

    pub fn set_objective(&mut self, obj: &DVector) -> Result<(), Error> {
        let inds = (0..obj.len() as c_int).collect::<Vec<_>>();
        Ok(trycpx!(cpx::NETchgobj(
            cpx::env(),
            self.net,
            obj.len() as c_int,
            inds.as_ptr(),
            obj.as_ptr()
        )))

    }

    pub fn add_arc(&mut self, src: usize, snk: usize, cost: Real, cap: Real) -> Result<(), Error> {
        let f = src as c_int;
        let t = snk as c_int;
        let c = cost as c_double;
        let u = cap as c_double;
        let name = CString::new(format!("x{}#{}_{}", self.num_arcs() + 1, f + 1, t + 1)).unwrap();
        let cname = name.as_ptr();
        Ok(trycpx!(cpx::NETaddarcs(
            cpx::env(),
            self.net,
            1,
            &f,
            &t,
            ptr::null(),
            &u,
            &c,
            &cname as *const *const c_char
        )))

    }

    pub fn solve(&mut self) -> Result<(), Error> {
        Ok(trycpx!(cpx::NETprimopt(cpx::env(), self.net)))

    }

    pub fn objective(&self) -> Result<Real, Error> {
        let mut objval: c_double = 0.0;
        trycpx!(cpx::NETgetobjval(
            cpx::env(),
            self.net,
            &mut objval as *mut c_double
        ));
        Ok(objval)
    }

    pub fn get_solution(&self) -> Result<DVector, Error> {
        let mut sol = dvec![0.0; self.num_arcs()];
        let mut stat: c_int = 0;
        let mut objval: c_double = 0.0;

            ptr::null_mut()
        ));
        Ok(sol)
    }

    pub fn writelp(&self, filename: &str) -> Result<(), Error> {
        let fname = CString::new(filename).unwrap();
        Ok(trycpx!(cpx::NETwriteprob(
            cpx::env(),
            self.net,
            fname.as_ptr(),
            ptr::null_mut()
        )))
    }
}

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

    }
}

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

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

    /// Combines this minorant with another minorant.
    pub fn combine(&self, self_factor: Real, other_factor: Real, other: &Minorant) -> Minorant {
        Minorant {
            constant: self_factor * self.constant + other_factor * other.constant,
            linear: self.linear
                .combine(self_factor, other_factor, &other.linear),
        }
    }

    /// Combines several minorants storing the result in this minorant.
    pub fn combine_all(&mut self, factors: &[Real], minorants: &[Minorant]) {
        debug_assert_eq!(factors.len(), minorants.len());
        self.constant = 0.0;

// Copyright (c) 2016, 2017 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 {DVector, Real};
use {Evaluation, FirstOrderProblem, HKWeighter, Update};

use master::{BoxedMasterProblem, MasterProblem, UnconstrainedMasterProblem};
use master::{CplexMaster, MinimalMaster};

use std::mem::swap;
use std::f64::{INFINITY, NEG_INFINITY};
use std::time::Instant;
use std::result::Result;


use failure::Error;

/// A solver error.
#[derive(Debug, Fail)]
pub enum SolverError {
    /// An error occured during oracle evaluation.

    #[fail(display = "Parameter error: {}", _0)]
    Parameter(String),
    /// The lower bound of a variable is larger than the upper bound.
    #[fail(display = "Invalid bounds, lower:{} upper:{}", lower, upper)]
    InvalidBounds { lower: Real, upper: Real },
    /// The value of a variable is outside its bounds.
    #[fail(display = "Violated bounds, lower:{} upper:{} value:{}", lower, upper, value)]
    ViolatedBounds {
        lower: Real,
        upper: Real,
        value: Real,
    },
    /// The variable index is out of bounds.
    #[fail(display = "Variable index out of bounds, got:{} must be < {}", index, nvars)]
    InvalidVariable { index: usize, nvars: usize },
    /// Iteration limit has been reached.
    #[fail(display = "The iteration limit of {} has been reached.", limit)]
    IterationLimit { limit: usize },
}

     */
    pub step: Step,
}

impl<'a> BundleState<'a> {}

macro_rules! current_state {
    ($slf: ident, $step: expr) => {
        BundleState{
            cur_y : &$slf.cur_y,
            cur_val : $slf.cur_val,
            nxt_y : &$slf.nxt_y,
            nxt_mod : $slf.nxt_mod,
            nxt_val : $slf.nxt_val,
            new_cutval : $slf.new_cutval,
            sgnorm : $slf.sgnorm,
            weight: $slf.master.weight(),
            step: $step,
            expected_progress: $slf.expected_progress,
        }
    };
}

     * Note that the solver owns the problem, so you cannot use the
     * same problem description elsewhere as long as it is assigned to
     * the solver. However, it is possible to get a reference to the
     * internally stored problem using `Solver::problem()`.
     */
    pub fn new_params(problem: P, params: SolverParams) -> Result<Solver<P, Pr, E>, SolverError> {
        Ok(Solver {
            problem: problem,
            params: params,
            terminator: Box::new(StandardTerminator {
                termination_precision: 1e-3,
            }),
            weighter: Box::new(HKWeighter::new()),
            bounds: vec![],
            cur_y: dvec![],
            cur_val: 0.0,

            nxt_mods: dvec![],
            new_cutval: 0.0,
            sgnorm: 0.0,
            expected_progress: 0.0,
            cnt_descent: 0,
            cnt_null: 0,
            start_time: Instant::now(),
            master: Box::new(BoxedMasterProblem::new(MinimalMaster::new()
                .map_err(SolverError::Master)?)),

            minorants: vec![],
            iterinfos: vec![],
        })
    }

    /// A new solver with default parameter.
    pub fn new(problem: P) -> Result<Solver<P, Pr, E>, SolverError> {

                        upper
                    } else {
                        0.0
                    };
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, value));
                }
                Update::AddVariableValue {
                    lower,
                    upper,
                    value,
                } => {
                    if lower > upper {
                        return Err(SolverError::InvalidBounds { lower, upper });
                    }
                    if value < lower || value > upper {
                        return Err(SolverError::ViolatedBounds {
                            lower,
                            upper,
                            value,
                        });
                    }
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, value));
                }
                Update::MoveVariable { index, value } => {
                    if index >= self.bounds.len() {
                        return Err(SolverError::InvalidVariable {
                            index,
                            nvars: self.bounds.len(),
                        });
                    }
                    let (lower, upper) = self.bounds[index];
                    if value < lower || value > upper {
                        return Err(SolverError::ViolatedBounds {
                            lower,
                            upper,
                            value,
                        });
                    }
                    newvars.push((Some(index), 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.cur_y[index] = val;
                self.nxt_y[index] = val;
                self.nxt_d[index] = 0.0;
            }
            // add new variables
            self.cur_y
                .extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3));
            self.nxt_y
                .extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3));
            self.nxt_d.resize(self.nxt_y.len(), 0.0);
            Ok(true)
        } else {
            Ok(false)
        }
    }

    }

    fn show_info(&self, step: Step) {
        let time = self.start_time.elapsed();
        info!(
            "{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4} {:4} {:4}{:1}  {:9.4} {:9.4} \
             {:12.6e}({:12.6e}) {:12.6e}",
            if step == Step::Term {
                "_endit"
            } else {
                "endit "
            },
            time.as_secs() / 3600,
            (time.as_secs() / 60) % 60,
            time.as_secs() % 60,
            time.subsec_nanos() / 10_000_000,
            self.cnt_descent,
            self.cnt_descent + self.cnt_null,
            self.master.cnt_updates(),

     * information.
     */
    fn init_master(&mut self) -> Result<(), SolverError> {
        let m = self.problem.num_subproblems();

        self.master = if m == 1 && self.params.max_bundle_size == 2 {
            debug!("Use minimal master problem");
            Box::new(BoxedMasterProblem::new(MinimalMaster::new()
                .map_err(SolverError::Master)?))

        } else {
            debug!("Use CPLEX master problem");
            Box::new(BoxedMasterProblem::new(CplexMaster::new()
                .map_err(SolverError::Master)?))

        };

        let lb = self.problem.lower_bounds().map(DVector);
        let ub = self.problem.upper_bounds().map(DVector);


        if lb.as_ref()
            .map(|lb| lb.len() != self.problem.num_variables())
            .unwrap_or(false)
        {
            return Err(SolverError::Dimension);
        }

        if ub.as_ref()
            .map(|ub| ub.len() != self.problem.num_variables())
            .unwrap_or(false)
        {
            return Err(SolverError::Dimension);
        }

        self.master
            .set_num_subproblems(m)
            .map_err(SolverError::Master)?;
        self.master
            .set_vars(self.problem.num_variables(), lb, ub)
            .map_err(SolverError::Master)?;
        self.master
            .set_max_updates(self.params.max_updates)
            .map_err(SolverError::Master)?;

        self.minorants = (0..m).map(|_| vec![]).collect();

        self.cur_val = 0.0;
        for i in 0..m {
            let result = self.problem

                .evaluate(i, &self.cur_y, INFINITY, 0.0)
                .map_err(SolverError::Evaluation)?;
            self.cur_vals[i] = result.objective();
            self.cur_val += self.cur_vals[i];

            let mut minorants = result.into_iter();
            if let Some((minorant, primal)) = minorants.next() {
                self.cur_mods[i] = minorant.constant;
                self.cur_mod += self.cur_mods[i];
                self.minorants[i].push(MinorantInfo {
                    index: self.master
                        .add_minorant(i, minorant)
                        .map_err(SolverError::Master)?,
                    multiplier: 0.0,
                    primal: Some(primal),
                });
            } else {
                return Err(SolverError::NoMinorant);
            }
        }

        self.cur_valid = true;

        // Solve the master problem once to compute the initial
        // subgradient.
        //
        // We could compute that subgradient directly by
        // adding up the initial minorants, but this would not include
        // the eta terms. However, this is a heuristic anyway because
        // we assume an initial weight of 1.0, which, in general, will
        // *not* be the initial weight for the first iteration.
        self.master.set_weight(1.0).map_err(SolverError::Master)?;
        self.master
            .solve(self.cur_val)
            .map_err(SolverError::Master)?;
        self.sgnorm = self.master.get_dualoptnorm2().sqrt();

        // Compute the real initial weight.
        let state = current_state!(self, Step::Term);
        let new_weight = self.weighter.weight(&state, &self.params);
        self.master
            .set_weight(new_weight)
            .map_err(SolverError::Master)?;

        debug!("Init master completed");

        Ok(())
    }

    /// Solve the model (i.e. master problem) to compute the next candidate.
    fn solve_model(&mut self) -> Result<(), SolverError> {
        self.master
            .solve(self.cur_val)
            .map_err(SolverError::Master)?;
        self.nxt_d = self.master.get_primopt();
        self.nxt_y.add(&self.cur_y, &self.nxt_d);
        self.nxt_mod = self.master.get_primoptval();
        self.sgnorm = self.master.get_dualoptnorm2().sqrt();
        self.expected_progress = self.cur_val - self.nxt_mod;

        // update multiplier from master solution

        for i in 0..self.problem.num_subproblems() {
            let n = self.master.num_minorants(i);
            if n >= self.params.max_bundle_size {
                // aggregate minorants with smallest coefficients
                self.minorants[i].sort_by_key(|m| -((1e6 * m.multiplier) as isize));
                let aggr = self.minorants[i].split_off(self.params.max_bundle_size - 2);
                let aggr_sum = aggr.iter().map(|m| m.multiplier).sum();
                let (aggr_mins, aggr_primals): (Vec<_>, Vec<_>) = aggr.into_iter()
                    .map(|m| (m.index, m.primal.unwrap()))
                    .unzip();
                let (aggr_min, aggr_coeffs) = self.master
                    .aggregate(i, &aggr_mins)
                    .map_err(SolverError::Master)?;
                // append aggregated minorant
                self.minorants[i].push(MinorantInfo {
                    index: aggr_min,
                    multiplier: aggr_sum,
                    primal: Some(
                        self.problem.aggregate_primals(
                            aggr_coeffs
                                .into_iter()
                                .zip(aggr_primals.into_iter())
                                .collect(),
                        ),
                    ),
                });
            }
        }
        Ok(())
    }

    /// Perform a descent step.
    fn descent_step(&mut self) -> Result<(), SolverError> {
        let new_weight = self.weighter
            .weight(&current_state!(self, Step::Descent), &self.params);
        self.master
            .set_weight(new_weight)
            .map_err(SolverError::Master)?;
        self.cnt_descent += 1;
        swap(&mut self.cur_y, &mut self.nxt_y);
        swap(&mut self.cur_val, &mut self.nxt_val);
        swap(&mut self.cur_mod, &mut self.nxt_mod);
        swap(&mut self.cur_vals, &mut self.nxt_vals);
        swap(&mut self.cur_mods, &mut self.nxt_mods);
        self.master.move_center(1.0, &self.nxt_d);
        debug!("Descent Step");
        debug!("  dir ={}", self.nxt_d);
        debug!("  newy={}", self.cur_y);
        Ok(())
    }

    /// Perform a null step.
    fn null_step(&mut self) -> Result<(), SolverError> {
        let new_weight = self.weighter
            .weight(&current_state!(self, Step::Null), &self.params);
        self.master
            .set_weight(new_weight)
            .map_err(SolverError::Master)?;
        self.cnt_null += 1;
        debug!("Null Step");
        Ok(())
    }

    /// Perform one bundle iteration.
    #[cfg_attr(feature = "cargo-clippy", allow(collapsible_if))]
    pub fn step(&mut self) -> Result<Step, SolverError> {
        self.iterinfos.clear();

        if !self.cur_valid {
            // current point needs new evaluation
            self.init_master()?;
        }

        self.solve_model()?;
        if self.terminator

            .terminate(&current_state!(self, Step::Term), &self.params)
        {
            return Ok(Step::Term);
        }

        let m = self.problem.num_subproblems();
        let descent_bnd = self.get_descent_bound();
        let nullstep_bnd = if m == 1 {
            self.get_nullstep_bound()
        } else {
            INFINITY
        };
        let relprec = if m == 1 {
            self.get_relative_precision()
        } else {
            0.0
        };

        self.compress_bundle()?;

        let mut nxt_lb = 0.0;
        let mut nxt_ub = 0.0;
        self.new_cutval = 0.0;
        for fidx in 0..self.problem.num_subproblems() {
            let result = self.problem

                .evaluate(fidx, &self.nxt_y, nullstep_bnd, relprec)
                .map_err(SolverError::Evaluation)?;
            let fun_ub = result.objective();

            let mut minorants = result.into_iter();
            let mut nxt_minorant;
            let nxt_primal;

            nxt_ub += fun_ub;
            self.nxt_vals[fidx] = fun_ub;

            // move center of minorant to cur_y
            nxt_minorant.move_center(-1.0, &self.nxt_d);
            self.new_cutval += nxt_minorant.constant;
            self.minorants[fidx].push(MinorantInfo {
                index: self.master

                    .add_minorant(fidx, nxt_minorant)
                    .map_err(SolverError::Master)?,
                multiplier: 0.0,
                primal: Some(nxt_primal),
            });
        }

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

//! Finite-dimensional sparse and dense vectors.

use Real;
use std::fmt;
use std::ops::{Deref, DerefMut};

// use std::cmp::min;
use std::iter::FromIterator;
use std::vec::IntoIter;

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

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

impl fmt::Display for Vector {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        match *self {
            Vector::Dense(ref v) => write!(f, "{}", v),
            Vector::Sparse { size, ref elems } => {

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

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

                        if elems.last_mut().unwrap().1 == 0.0 {
                            elems.pop();
                            last_idx = n;
                        }
                    }
                }
            }
            Vector::Sparse {
                size: n,
                elems: elems,
            }
        }
    }

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