RsBundle  Check-in [998cbd7227]

indent_style="Block"
max_width=120
trailing_comma="Vertical"
attributes_on_same_line_as_field=false

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

extern crate bundle;
extern crate env_logger;
#[macro_use]
extern crate log;


use bundle::{FirstOrderProblem, Solver, SolverParams, StandardTerminator};
use bundle::mcf;
use std::env;

fn main() {
    env_logger::init().unwrap();

    let mut args = env::args();
    let program = args.next().unwrap();

    if let Some(filename) = args.next() {
        info!("Reading instance: {}", filename);
        let mut mmcf = mcf::MMCFProblem::read_mnetgen(&filename).unwrap();
        mmcf.multimodel = false;

        let mut solver = Solver::new_params(
            mmcf,
            SolverParams {
                max_bundle_size: 25,
                min_weight: 1e-3,
                max_weight: 100.0,
                ..Default::default()
            },
        ).unwrap();
        solver.terminator = Box::new(StandardTerminator {
            termination_precision: 1e-6,
        });
        solver.solve().unwrap();

        let costs: f64 = (0..solver.problem().num_subproblems())
            .map(|i| {
                let primals = solver.aggregated_primals(i);
                let aggr_primals = solver.problem().aggregate_primals_ref(&primals);
                solver.problem().get_primal_costs(i, &aggr_primals)
            })
            .sum();
        info!("Primal costs: {}", costs);
    } else {
        panic!("Usage: {} FILENAME", program);
    }
}

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

#[macro_use]
extern crate bundle;

extern crate env_logger;
extern crate failure;
#[macro_use]
extern crate log;

use bundle::{DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation, Solver, SolverParams};
use failure::Error;


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<'a> FirstOrderProblem<'a> for QuadraticProblem {
    type Primal = ();
    type EvalResult = SimpleEvaluation<()>;

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

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

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

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

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

fn main() {
    env_logger::init().unwrap();

    let f = QuadraticProblem::new();
    let mut solver = Solver::new_params(
        f,
        SolverParams {
            min_weight: 1.0,
            max_weight: 1.0,
            ..Default::default()
        },
    ).unwrap();
    solver.solve().unwrap();
}

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

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

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

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

/**

    AddVariableValue {
        lower: Real,
        upper: Real,
        value: Real,
    },
    /// Change the current value of a variable. The bounds remain
    /// unchanged.
    MoveVariable { index: usize, value: Real },



}

/**
 * Trait for implementing a first-order problem description.
 *
 */
pub trait FirstOrderProblem<'a> {

     * true function value within $relprec \cdot (\\|f(y)\\| + 1.0)$,
     * otherwise the returned objective should be the maximum of all
     * linear minorants at $y$.
     *
     * Note that `nullstep_bound` and `relprec` are usually only
     * useful if there is only a `single` subproblem.
     */
    fn evaluate(
        &'a mut self,
        i: usize,
        y: &[Real],
        nullstep_bound: Real,
        relprec: Real,
    ) -> Result<Self::EvalResult, Error>;

    /// Aggregate primal information.
    ///
    /// This function is called from the solver when minorants are
    /// aggregated. The problem can use this information to aggregate
    /// the corresponding primal information.
    ///

//! The procedure is described in
//!
//! > 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

        debug!("HKWeighter {:?} iter:{}", state.step, self.iter);

        if state.step == Step::Term {
            self.eps_weight = 1e30;
            self.iter = 0;
            return if state.cur_y.len() == 0 || state.sgnorm < state.cur_y.len() as Real * 1e-10 {
                1.0
            } else {
                state.sgnorm.max(1e-4)

            }.max(params.min_weight)
                .min(params.max_weight);
        }

        let cur_nxt = state.cur_val - state.nxt_val;
        let cur_mod = state.cur_val - state.nxt_mod;
        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 {
                self.iter = -1
            } else {
                self.iter = min(self.iter - 1, -1);
            }

            debug!(
                "  sgnorm={} cur_val={} new_cutval={} lin_err={} eps_weight={}",
                sgnorm,
                state.cur_val,
                state.new_cutval,
                lin_err,
                self.eps_weight
            );
            debug!("  new_weight={}", new_weight);

            new_weight
        } else {
            self.model_max = self.model_max.max(state.nxt_mod);
            let new_weight = if self.iter > 0 && cur_nxt > self.m_r * cur_mod {
                w
            } else if self.iter > 3 || state.nxt_val < self.model_max {
                state.weight / 2.0
            } else {
                state.weight

            }.max(state.weight / FACTOR)
                .max(params.min_weight);
            self.eps_weight = self.eps_weight.max(2.0 * cur_mod);
            if new_weight < state.weight {
                self.iter = 1;
                self.model_max = NEG_INFINITY;
            } else {
                self.iter = max(self.iter + 1, 1);

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

//! Proximal bundle method implementation.

#[macro_use]
extern crate const_cstr;
extern crate failure;
#[macro_use]
extern crate failure_derive;

extern crate itertools;

#[macro_use]
extern crate log;

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

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

pub mod minorant;
pub use minorant::Minorant;

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

pub mod solver;
pub use solver::{BundleState, IterationInfo, Solver, SolverParams, StandardTerminator, Step, Terminator, UpdateState,
                 Weighter};

mod hkweighter;
pub use hkweighter::HKWeighter;

mod master;

pub mod mcf;

// 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 std::result::Result;
use failure::Error;

pub trait MasterProblem {
    /// Unique index for a minorant.
    type MinorantIndex: Copy + Eq;

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

    /// Add or movesome variables with bounds.
    ///
    /// If an index is specified, existing variables are moved,
    /// otherwise new variables are generated.
    fn add_vars(
        &mut self,
        bounds: &[(Option<usize>, Real, Real)],
        extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector,
    ) -> Result<(), Error>;


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

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

     * d) without the influence of $\eta$.
     */
    fn update_box_multipliers(&mut self) -> bool {
        let mut updated_eta = false;
        let weight = self.master.weight();

        self.eta.resize(self.lb.len(), 0.0);
        for (&lb, &ub, x, eta) in multizip((
            self.lb.iter(),
            self.ub.iter(),
            self.primopt.iter_mut(),
            self.eta.iter_mut(),

        )) {
            let newx = if *x < lb {
                lb
            } else if *x > ub {
                ub
            } else {
                *eta = 0.0;
                continue;

    // 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.
    fn eta_cutval(&self) -> Real {
        multizip((self.lb.iter(), self.ub.iter(), self.eta.iter()))
            .map(|(&lb, &ub, &eta)| {
                if eta >= 0.0 {
                    if lb > NEG_INFINITY {
                        lb * eta
                    } else {
                        0.0
                    }
                } else if ub < INFINITY {
                    ub * eta
                } else {
                    0.0
                }




            })
            .sum()
    }

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

        self.master.weight()
    }

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

    fn add_vars(
        &mut self,
        bounds: &[(Option<usize>, Real, Real)],
        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(())
        }
    }

    #[cfg_attr(feature = "cargo-clippy", allow(cyclomatic_complexity))]
    fn solve(&mut self, center_value: Real) -> Result<(), Error> {
        debug!("Solve Master");
        debug!("  lb      ={}", self.lb);
        debug!("  ub      ={}", self.ub);

        if self.need_new_candidate {
            self.compute_candidate();
        }

        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!("  modval={}", self.master.eval_model(&self.primopt));
            debug!("  augval={}", augval);
            debug!("  cutval={}", cutval);
            debug!("  model_prec={}", model_prec);
            debug!("  old_augval={}", old_augval);
            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);
        }

        debug!("Model");

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

//! Master problem implementation using CPLEX.

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

use cplex_sys as cpx;

use std::ptr;
use std::os::raw::{c_char, c_int};
use std::f64::{self, NEG_INFINITY};
use std::result::Result;
use failure::Error;

/// A solver error.
#[derive(Debug, Fail)]
pub enum CplexMasterError {
    #[fail(display = "Solver Error: no minorants when solving the master problem")] NoMinorants,

}

pub struct CplexMaster {
    lp: *mut cpx::Lp,

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

    }

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

    fn set_num_subproblems(&mut self, n: usize) -> Result<(), Error> {
        trycpx!(cpx::setintparam(
            cpx::env(),
            cpx::Param::Qpmethod.to_c(),
            cpx::Alg::Barrier.to_c()
        ));
        trycpx!(cpx::setdblparam(
            cpx::env(),
            cpx::Param::Barepcomp.to_c(),
            1e-12
        ));

        self.min2index = vec![vec![]; n];
        self.index2min.clear();
        self.freeinds.clear();
        self.minorants = vec![vec![]; n];
        self.opt_mults = vec![dvec![]; n];

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

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

        let min_idx = self.minorants[fidx].len();
        self.minorants[fidx].push(minorant);
        self.opt_mults[fidx].push(0.0);

        self.force_update = true;

            self.min2index[fidx].push(idx);
            self.index2min.push((fidx, min_idx));
            self.updateinds.push(idx);
            Ok(idx)
        }
    }

    fn add_vars(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector,
    ) -> Result<(), Error> {

        debug_assert!(!self.minorants[0].is_empty());
        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() {
            if !mins.is_empty() {
                for (i, m) in mins.iter_mut().enumerate() {
                    let new_subg = extend_subgradient(fidx, i, &changedvars);
                    for (&j, &g) in changed.iter().zip(new_subg.iter()) {
                        m.linear[j] = g;
                    }

        for (fidx_i, mins_i) in self.minorants.iter().enumerate() {
            for (i, m_i) in mins_i.iter().enumerate() {
                let idx_i = self.min2index[fidx_i][i];
                for (fidx_j, mins_j) in self.minorants.iter().enumerate() {
                    for (j, m_j) in mins_j.iter().enumerate() {
                        let idx_j = self.min2index[fidx_j][j];
                        if idx_i <= idx_j {
                            let x: Real = (nnewvars..noldvars)
                                .map(|k| m_i.linear[k] * m_j.linear[k])
                                .sum();
                            self.qterm[idx_i][idx_j] += x;
                            self.qterm[idx_j][idx_i] = self.qterm[idx_i][idx_j];
                        }
                    }
                }
            }
        }

            let mut inds = Vec::with_capacity(nvars);
            for mins in &self.minorants {
                for m in mins {
                    inds.push(c.len() as c_int);
                    c.push(-m.constant * self.weight - m.linear.dot(eta));
                }
            }
            trycpx!(cpx::chgobj(
                cpx::env(),
                self.lp,
                nvars as c_int,
                inds.as_ptr(),
                c.as_ptr()
            ));
        }

        trycpx!(cpx::qpopt(cpx::env(), self.lp));
        let mut sol = vec![0.0; nvars];
        trycpx!(cpx::getx(
            cpx::env(),
            self.lp,
            sol.as_mut_ptr(),
            0,
            nvars as c_int - 1
        ));

        let mut idx = 0;
        let mut mults = Vec::with_capacity(nvars);
        let mut mins = Vec::with_capacity(nvars);
        for fidx in 0..self.minorants.len() {
            for i in 0..self.minorants[fidx].len() {
                self.opt_mults[fidx][i] = sol[idx];

                for _ in 0..self.minorants[i].len() {
                    rmatind.push(nvars as c_int);
                    rmatval.push(1.0);
                    nvars += 1;
                }
            }

            trycpx!(cpx::addrows(
                cpx::env(),
                self.lp,
                nvars as c_int,
                nfun as c_int,
                nvars as c_int,
                rhs.as_ptr(),
                sense.as_ptr(),
                rmatbeg.as_ptr(),
                rmatind.as_ptr(),
                rmatval.as_ptr(),
                ptr::null(),
                ptr::null()
            ));
        }

        // build quadratic term
        {
            self.qterm.resize(self.index2min.len(), dvec![]);
            for i in 0..self.qterm.len() {
                self.qterm[i].resize(self.index2min.len(), 0.0);

            }
            maxq *= 1e-8;

            // update coefficients
            for (i, &idx_i) in activeinds.iter().enumerate() {
                for (j, &idx_j) in activeinds.iter().enumerate() {
                    if i != j {
                        trycpx!(cpx::chgqpcoef(
                            cpx::env(),
                            self.lp,
                            i as c_int,
                            j as c_int,
                            self.qterm[idx_i][idx_j]
                        ));
                    } else {
                        trycpx!(cpx::chgqpcoef(
                            cpx::env(),
                            self.lp,
                            i as c_int,
                            j as c_int,
                            self.qterm[idx_i][idx_j] + maxq
                        ));
                    }
                }
            }
        }

        self.updateinds.clear();
        self.force_update = false;

        Ok(())
    }
}

// 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::UnconstrainedMasterProblem;

use std::f64::NEG_INFINITY;
use std::result::Result;

use failure::Error;


/// Minimal master problem error.
#[derive(Debug, Fail)]
pub enum MinimalMasterError {
    #[fail(display = "Solver Error: too many subproblems (got: {} must be <= 2)", nsubs)]
    NumSubproblems {
        nsubs: usize,
    },
    #[fail(display = "Solver Error: the minimal master problem allows at most two minorants")] MaxMinorants,

    #[fail(display = "Solver Error: no minorants when solving the master problem")] NoMinorants,

}

/**
 * A minimal master problem with only two minorants.
 *
 * This is the simplest possible master problem for bundle methods. It
 * has only two minorants and only one function model. The advantage

            return Err(MinimalMasterError::MaxMinorants.into());
        }
        self.minorants.push(minorant);
        self.opt_mult.push(0.0);
        Ok(self.minorants.len() - 1)
    }

    fn add_vars(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector,
    ) -> Result<(), Error> {

        if !self.minorants.is_empty() {
            let noldvars = self.minorants[0].linear.len();
            let mut changedvars = vec![];
            changedvars.extend_from_slice(changed);
            changedvars.extend(noldvars..noldvars + nnew);
            for (i, m) in self.minorants.iter_mut().enumerate() {
                let new_subg = extend_subgradient(0, i, &changedvars);
                for (&j, &g) in changed.iter().zip(new_subg.iter()) {
                    m.linear[j] = g;
                }
                m.linear.extend_from_slice(&new_subg[changed.len()..]);
            }

// 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 std::result::Result;
use failure::Error;

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

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

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

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

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<Self::MinorantIndex, Error>;

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

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

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

// 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)]
pub struct MCFFormatError {
    msg: String,
}

#[derive(Clone, Copy, Debug)]
struct ArcInfo {
    arc: usize,
    src: usize,
    snk: usize,
}

impl MMCFProblem {
    pub fn read_mnetgen(basename: &str) -> Result<MMCFProblem, Error> {
        let mut buffer = String::new();
        {
            let mut f = try!(File::open(&format!("{}.nod", basename)));
            try!(f.read_to_string(&mut buffer));
        }
        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 arc = try!(data.next().unwrap().parse::<usize>()) - 1;
            let src = try!(data.next().unwrap().parse::<usize>()) - 1;
            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,
            });
            // add arc
            try!(nets[com].add_arc(
                src,
                snk,
                cost,
                if cap < 0.0 { INFINITY } else { cap }
            ));
            // set objective
            cbase[com].push(cost); // + 1e-6 * coeff
                                   // add to mutual capacity constraint
            if mt >= 0 {
                lhsidx[mt as usize][com].push(coeff);
            }
        }

        // read rhs of coupling constraints
        {

            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 {
            primals[0]
                .iter()
                .enumerate()
                .map(|(i, x)| x * self.cbase[fidx][i])
                .sum()
        } else {
            let mut sum = 0.0;
            for (fidx, p) in primals.iter().enumerate() {
                for (i, x) in p.iter().enumerate() {
                    sum += x * self.cbase[fidx][i];
                }
            }

    }

    fn upper_bounds(&self) -> Option<Vec<Real>> {
        None
    }

    fn num_subproblems(&self) -> usize {
        if self.multimodel {
            self.nets.len()
        } else {
            1
        }
    }

    #[allow(unused_variables)]
    fn evaluate(
        &'a mut self,
        fidx: usize,
        y: &[Real],
        nullstep_bound: Real,
        relprec: Real,
    ) -> Result<Self::EvalResult, Error> {
        // compute costs
        self.rhsval = 0.0;
        for i in 0..self.c.len() {
            self.c[i].clear();
            self.c[i].extend(self.cbase[i].iter());
        }

        for (i, &y) in y.iter().enumerate().filter(|&(i, &y)| y != 0.0) {
            self.rhsval += self.rhs[i] * y;
            for (fidx, c) in self.c.iter_mut().enumerate() {
                for elem in &self.lhs[i][fidx] {
                    c[elem.ind] += y * elem.val;
                }
            }
        }

                objective = self.rhsval - try!(self.nets[fidx].objective());
            } 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<_>>())
    }
}

// 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, Real};

use cplex_sys as cpx;

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

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

use failure::{err_msg, Error};

pub struct Solver {
    net: *mut cpx::Net,
    logfile: *mut cpx::File,
}

        let mut status: c_int;
        let mut net = ptr::null_mut();
        let logfile;

        unsafe {
            #[cfg_attr(feature = "cargo-clippy", allow(never_loop))]
            loop {
                logfile = cpx::fopen(
                    const_cstr!("mcf.cpxlog").as_ptr(),
                    const_cstr!("w").as_ptr(),
                );
                if logfile.is_null() {
                    return Err(err_msg("Can't open log-file"));
                }
                status = cpx::setlogfile(cpx::env(), logfile);
                if status != 0 {
                    break;
                }

                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);
                cpx::fclose(logfile);
                return Err(
                    cpx::CplexError {
                        code: status,
                        msg: CString::from_raw(msg).to_string_lossy().into_owned(),
                    }.into(),
                );
            }
        }

        Ok(Solver {
            net: net,
            logfile: logfile,
        })
    }

    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;
        trycpx!(cpx::NETsolution(
            cpx::env(),
            self.net,
            &mut stat as *mut c_int,
            &mut objval as *mut c_double,
            sol.as_mut_ptr(),
            ptr::null_mut(),
            ptr::null_mut(),
            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()
        )))
    }
}

//
// 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 {DVector, Real};

use std::fmt;

/**
 * A linear minorant of a convex function.
 *
 * A linear minorant of a convex function $f \colon \mathbb{R}\^n \to

        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;

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

//! The main bundle method solver.

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;

    #[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 max_updates: usize,
}

impl SolverParams {
    /// Verify that all parameters are valid.
    fn check(&self) -> Result<(), SolverError> {
        if self.max_bundle_size < 2 {
            Err(SolverError::Parameter(format!(
                "max_bundle_size must be >= 2 (got: {})",
                self.max_bundle_size
            )))
        } else if self.acceptance_factor <= 0.0 || self.acceptance_factor >= 1.0 {
            Err(SolverError::Parameter(format!(
                "acceptance_factor must be in (0,1) (got: {})",
                self.acceptance_factor
            )))
        } else if self.nullstep_factor <= 0.0 || self.nullstep_factor > self.acceptance_factor {
            Err(SolverError::Parameter(format!(
                "nullstep_factor must be in (0,acceptance_factor] (got: {}, acceptance_factor:{})",
                self.nullstep_factor,
                self.acceptance_factor
            )))
        } else if self.min_weight <= 0.0 {
            Err(SolverError::Parameter(format!(
                "min_weight must be in > 0 (got: {})",
                self.min_weight
            )))
        } else if self.max_weight < self.min_weight {
            Err(SolverError::Parameter(format!(
                "max_weight must be in >= min_weight (got: {}, min_weight: {})",
                self.max_weight,
                self.min_weight
            )))
        } else if self.max_updates == 0 {
            Err(SolverError::Parameter(format!(
                "max_updates must be in > 0 (got: {})",
                self.max_updates
            )))
        } else {
            Ok(())
        }
    }
}

impl Default for SolverParams {

    }
}

/**
 * Implementation of a bundle method.
 */
pub struct Solver<P, Pr, E>
where
    P: for<'a> FirstOrderProblem<'a, Primal = Pr, EvalResult = E>,
    E: Evaluation<Pr>,
{
    /// The first order problem description.
    problem: P,

    /// The solver parameter.
    pub params: SolverParams,

    /// Accumulated information about the last iteration.
    iterinfos: Vec<IterationInfo>,
}


impl<P, Pr, E> Solver<P, Pr, E>
where
    P: for<'a> FirstOrderProblem<'a, Primal = Pr, EvalResult = E>,
    E: Evaluation<Pr>,
{
    /**
     * Create a new solver for the given problem.
     *
     * 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,
            cur_mod: 0.0,
            cur_vals: dvec![],
            cur_mods: dvec![],
            cur_valid: false,
            nxt_d: dvec![],
            nxt_y: dvec![],
            nxt_val: 0.0,
            nxt_mod: 0.0,
            nxt_vals: dvec![],
            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> {

        self.bounds.reserve(self.cur_y.len());
        for i in 0..self.cur_y.len() {
            let lb_i = lb.as_ref().map(|x| x[i]).unwrap_or(NEG_INFINITY);
            let ub_i = ub.as_ref().map(|x| x[i]).unwrap_or(INFINITY);
            if lb_i > ub_i {
                return Err(SolverError::InvalidBounds {
                    lower: lb_i,
                    upper: ub_i,
                });
            }
            if self.cur_y[i] < lb_i {
                self.cur_valid = false;
                self.cur_y[i] = lb_i;
            } else if self.cur_y[i] > ub_i {
                self.cur_valid = false;

            let changed = self.update_problem(term)?;
            // do not stop if the problem has been changed
            if changed && term == Step::Term {
                term = Step::Null
            }
            self.show_info(term);
            if term == Step::Term {
                return Ok(true);
            }
        }
        Ok(false)
    }

    /// Called to update the problem.
    ///

                        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() {
            let mut problem = &mut self.problem;
            let minorants = &self.minorants;
            self.master
                .add_vars(
                    &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::<Vec<_>>(),
                    &mut move |fidx, minidx, vars| {
                        problem
                            .extend_subgradient(minorants[fidx][minidx].primal.as_ref().unwrap(), vars)
                            .map(DVector)
                            .unwrap()
                    },
                )
                .map_err(SolverError::Master)?;
            // 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)
        }
    }

            .iter()
            .map(|m| (m.multiplier, m.primal.as_ref().unwrap()))
            .collect()
    }

    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(),
            if step == Step::Descent { "*" } else { " " },
            self.master.weight(),
            self.expected_progress,
            self.nxt_mod,
            self.nxt_val,
            self.cur_val
        );
    }

    /// Return the current center of stability.
    pub fn center(&self) -> &[Real] {
        &self.cur_y
    }

     * 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

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

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

        if self.new_cutval > self.cur_val + 1e-3 {
            warn!(
                "New minorant has higher value in center new:{} old:{}",
                self.new_cutval,
                self.cur_val
            );
            self.cur_val = self.new_cutval;
            self.iterinfos.push(IterationInfo::NewMinorantTooHigh {
                new: self.new_cutval,
                old: self.cur_val,
            });
        }

    }

    /// 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 {
        self.iter().zip(other.iter()).map(|(a, b)| a * b).sum()
        // let mut ip = 0.0;
        // for i in 0..min(self.len(), other.len()) {
        //     ip += unsafe { self.get_unchecked(i) * other.get_unchecked(i) };
        // }
        // return ip;
    }

    /// 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));
        // self.resize(x.len(), 0.0);
        // for i in 0..x.len() {
        //     unsafe { *self.get_unchecked_mut(i) = *x.get_unchecked(i) + *y.get_unchecked(i) };
        // }
    }

    /// 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());
        for (x, y) in self.iter_mut().zip(y.iter()) {
            *x += alpha * y;
        }
        // for i in 0..self.len() {
        //     unsafe { *self.get_unchecked_mut(i) += alpha * *y.get_unchecked(i) };
        // }
    }

    /// 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) {
        for (x, y) in self.iter_mut().zip(y.iter()) {
            *x += alpha * y;
        }
        let n = self.len();
        if n < y.len() {
            self.extend_from_slice(&y[n..]);
        }
        // if self.len() < y.len() {
        //     self.resize(y.len(), 0.0);
        // }
        // for i in 0..y.len() {
        //     unsafe { *self.get_unchecked_mut(i) += alpha * *y.get_unchecked(i) };
        // }
    }


    /// Combines this vector with another vector.
    pub fn combine(&self, self_factor: Real, other_factor: Real, other: &DVector) -> DVector {
        assert_eq!(self.len(), other.len());
        let mut result = vec![];
        result.extend(
            self.iter()
                .zip(other.iter())
                .map(|(a, b)| self_factor * a + other_factor * b),
        );
        DVector(result)
        // let mut result = DVector(Vec::with_capacity(self.len()));
        // for i in 0..self.len() {
        //     result.push(unsafe {
        //         self_factor * *self.get_unchecked(i) +
        //             other_factor * *other.get_unchecked(i)
        //     });
        // }
        // result
    }


    /// Return the 2-norm of this vector.
    pub fn norm2(&self) -> Real {
        self.iter().map(|x| x * x).sum::<Real>().sqrt()
        // let mut norm = 0.0;
        // for x in self.iter() {
        //     norm += x * x
        // }
        // norm.sqrt()
    }
}

     * 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)
            }
        }
    }
}