RsBundle  Check-in [2bc4277460]

max_width = 200

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

//! Low-level CPLEX interface.
//!
//! This module contains plain, unsafe functions to the CPLEX
//! C-interface.
//!
//! Most CPLEX functions require an environment pointer as first

use std::ptr;
use std::fmt;
use std::error;

pub use std::result::Result as cplex_std_Result;
pub use std::convert::From as cplex_std_From;

pub const CPXMSGBUFSIZE: usize = 1024;
#[allow(dead_code)]
pub const CPX_MAX: c_int = -1;
pub const CPX_MIN: c_int = 1;

pub enum CPXenv {}
pub enum CPXnet {}
pub enum CPXlp {}
pub enum CPXfile {}

pub type CPXINT = i32;


pub const CPX_PARAM_QPMETHOD: c_int = 1063;
pub const CPX_PARAM_THREADS: c_int = 1067;
pub const CPX_PARAM_BAREPCOMP: c_int = 3002;

pub const CPX_ALG_AUTOMATIC: CPXINT = 0;
pub const CPX_ALG_PRIMAL: CPXINT = 1;
pub const CPX_ALG_DUAL: CPXINT = 2;
pub const CPX_ALG_NET: CPXINT = 3;
pub const CPX_ALG_BARRIER: CPXINT = 4;

/// Globally unique environment.
static mut env_: *mut CPXenv = 0 as *mut CPXenv;

pub unsafe fn env() -> *mut CPXenv {
    if env_ == ptr::null_mut() {
        let mut status: c_int = 0;
        env_ = CPXopenCPLEX(&mut status);
        if status != 0 {
            panic!("Can't open CPLEX environment");
        }
    }
    env_
}

#[allow(dead_code)]
extern "C" {
    pub fn CPXopenCPLEX(status: *mut c_int) -> *mut CPXenv;
    pub fn CPXcloseCPLEX(env: *mut *mut CPXenv) -> c_int;
    pub fn CPXgeterrorstring(env: *const CPXenv, errcode: c_int, buffer_str: *mut c_char) -> *const c_char;
    pub fn CPXfopen(filename: *const c_char, mode: *const c_char) -> *mut CPXfile;
    pub fn CPXfclose(file: *mut CPXfile) -> c_int;
    pub fn CPXsetlogfile(env: *mut CPXenv, file: *mut CPXfile) -> c_int;
    pub fn CPXsetintparam(env: *mut CPXenv, whichparam: c_int, newvalue: CPXINT) -> c_int;
    pub fn CPXsetdblparam(env: *mut CPXenv, whichparam: c_int, newvalue: c_double) -> c_int;

    pub fn CPXcreateprob(env: *mut CPXenv, status: *mut c_int, name: *const c_char) -> *mut CPXlp;
    pub fn CPXfreeprob(env: *mut CPXenv, net: *mut *mut CPXlp) -> c_int;
    pub fn CPXgetnumrows(env: *const CPXenv, lp: *const CPXlp) -> c_int;
    pub fn CPXgetnumcols(env: *const CPXenv, lp: *const CPXlp) -> c_int;
    pub fn CPXchgobj(env: *const CPXenv, lp: *mut CPXlp, cnt: c_int, indices: *const c_int, values: *const c_double) -> c_int;
    pub fn CPXaddrows(env: *const CPXenv,
                      lp: *mut CPXlp,
                      ccnt: c_int,
                      rcnt: c_int,
                      nzcnt: c_int,
                      rhs: *const c_double,
                      sense: *const c_char,
                      rmatbeg: *const c_int,
                      rmatind: *const c_int,
                      rmatval: *const c_double,
                      colname: *const *const c_char,
                      rowname: *const *const c_char)
                      -> c_int;

    pub fn CPXchgqpcoef(env: *const CPXenv, lp: *mut CPXlp, i: c_int, j: c_int, x: c_double) -> c_int;
    pub fn CPXqpopt(env: *const CPXenv, lp: *mut CPXlp) -> c_int;
    pub fn CPXgetx(env: *const CPXenv, lp: *const CPXlp, x: *mut c_double, begin: c_int, end: c_int) -> c_int;

    pub fn CPXNETcreateprob(env: *mut CPXenv, status: *mut c_int, name: *const c_char) -> *mut CPXnet;
    pub fn CPXNETfreeprob(env: *mut CPXenv, net: *mut *mut CPXnet) -> c_int;
    pub fn CPXNETgetnumnodes(env: *const CPXenv, net: *const CPXnet) -> c_int;
    pub fn CPXNETgetnumarcs(env: *const CPXenv, net: *const CPXnet) -> c_int;
    pub fn CPXNETaddnodes(env: *const CPXenv, net: *mut CPXnet, nnodes: c_int, supply: *const c_double, names: *const *const c_char) -> c_int;
    pub fn CPXNETaddarcs(env: *const CPXenv,
                         net: *mut CPXnet,
                         narcs: c_int,
                         fromnode: *const c_int,
                         tonode: *const c_int,
                         low: *const c_double,
                         up: *const c_double,
                         obj: *const c_double,
                         names: *const *const c_char)
                         -> c_int;
    pub fn CPXNETchgobjsen(env: *const CPXenv, net: *mut CPXnet, maxormin: c_int) -> c_int;
    pub fn CPXNETchgsupply(env: *const CPXenv, net: *mut CPXnet, cnt: c_int, indices: *const c_int, supply: *const c_double) -> c_int;
    pub fn CPXNETchgobj(env: *const CPXenv, net: *mut CPXnet, cnt: c_int, indices: *const c_int, obj: *const c_double) -> c_int;

    pub fn CPXNETprimopt(env: *const CPXenv, net: *mut CPXnet) -> c_int;
    pub fn CPXNETgetobjval(env: *const CPXenv, net: *const CPXnet, objval: *mut c_double) -> c_int;
    pub fn CPXNETsolution(env: *const CPXenv, net: *const CPXnet, netstat: *mut c_int, objval: *mut c_double, x: *mut c_double, pi: *mut c_double, slack: *mut c_double, dj: *mut c_double) -> c_int;
    pub fn CPXNETwriteprob(env: *const CPXenv, net: *const CPXnet, filename: *const c_char, format: *const c_char) -> c_int;
}


/// Error descriping a CPLEX status code.
///
/// This is a wrapper around CPLEX status code as returned by most
/// CPLEX low-level functions. The error struct contains the status
/// code itself along with the error message string returned by
/// `CPXgeterrorstring`.
#[derive(Debug)]
pub struct CplexError {
    /// The CPLEX error status code.
    pub code: c_int,
    /// The CPLEX error message associated with the status code.
    pub msg: String,
}

impl error::Error for CplexError {
    fn description(&self) -> &str {
        "CPLEX Error"
    }
}

// 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
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.

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

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

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

use std::error;
use std::vec::IntoIter;

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


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

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

    fn into_iter(self) -> Self::IntoIter {

/// variables. The possible updates are encoded in this type.
#[derive(Debug, Clone, Copy)]
pub enum Update {
    /// Add a variable with bounds.
    ///
    /// The initial value of the variable will be the feasible value
    /// closest to 0.
    AddVariable { lower: Real, upper: Real },
    /// Add a variable with bounds and initial value.
    AddVariableValue {
        lower: Real,
        upper: Real,
        value: Real,
    },
}

/**
 * Trait for implementing a first-order problem description.
 *
 */
pub trait FirstOrderProblem<'a> {
    /// Custom error type for evaluating this oracle.
    type Error: error::Error + 'static;

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

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

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

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

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

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

    /**
     * Evaluate the i^th subproblem at the given point.
     *
     * The returned evaluation result must contain (an upper bound on)
     * the objective value at $y$ as well as at least one subgradient
     * centered at $y$.

     * 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: &DVector, nullstep_bound: Real, relprec: Real) -> Result<Self::EvalResult, Self::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.
    ///

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


//! Weight updating rule according to Helmberg and Kiwiel.
//!
//! 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 {Weighter, BundleState, SolverParams, Step};

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

const FACTOR: Real = 2.0;

/**
 * Weight updating rule according to Helmberg and Kiwiel.
 *
 * The procedure is described in
 *
 * > Helmberg, C. and Kiwiel, K.C. (2002): A spectral bundle method
 * > with bounds, Math. Programming A 93, 173--194
 */
pub struct HKWeighter {
    eps_weight: Real,
    m_r: Real,
    iter: isize,
    model_max: Real,
}

impl HKWeighter {
    /// Create a new HKWeighter with default weight $m_R = 0.5$.
    pub fn new() -> HKWeighter {
        HKWeighter::new_weight(0.5)
    }

    /// 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 Weighter for HKWeighter {
    fn weight(&mut self, state: &BundleState, params: &SolverParams) -> Real {
        assert!(params.min_weight > 0.0);
        assert!(params.max_weight >= params.min_weight);

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

            return 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.weight / 2.0
                } else if 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);
            }

            debug!("  model_max={}", self.model_max);
            debug!("  new_weight={}", new_weight);

            return new_weight;
        }
    }
}

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

//! Proximal bundle method implementation.

#[macro_use]
extern crate quick_error;
extern crate libc;

pub mod minorant;
pub use minorant::Minorant;

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

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


mod hkweighter;
pub use hkweighter::HKWeighter;

mod master;

pub mod mcf;

// 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
// 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 {Real, DVector, Minorant};
use master::master::{MasterProblem, Error, Result};
use master::UnconstrainedMasterProblem;

use std::f64::{INFINITY, NEG_INFINITY};

/**
 * 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.
 */
pub struct BoxedMasterProblem<M: UnconstrainedMasterProblem> {
    lb: DVector,
    ub: DVector,
    eta: DVector,

    /// Primal optimal solution.
    primopt: DVector,

    /// Primal optimal solution value.
    primoptval: Real,

    /// Square of norm of dual optimal solution.
    dualoptnorm2: Real,

    /// Model precision.
    model_eps: Real,

    need_new_candidate: bool,

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

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

    /// The unconstrained master problem solver.
    master: M,
}


impl<M: UnconstrainedMasterProblem> BoxedMasterProblem<M> {
    pub fn new() -> Result<BoxedMasterProblem<M>> {
        Ok(BoxedMasterProblem {
            lb: dvec![],
            ub: dvec![],
            eta: dvec![],
            primopt: dvec![],
            primoptval: 0.0,
            dualoptnorm2: 0.0,
            model_eps: 0.6,
            max_updates: 100,
            cnt_updates: 0,
            need_new_candidate: true,
            master: match M::new() {
                Ok(m) => m,
                Err(e) => return Err(Error::Solver(Box::new(e))),
            },
        })
    }

    pub fn set_max_updates(&mut self, max_updates: usize) {

                if x <= b {
                    self.eta[i] = 0.0;
                    continue;
                }
            }
            self.primopt[i] = b;
            let neweta = (b - x) * weight;
            if neweta != self.eta[i] {
                updated_eta = true;
            }
            self.eta[i] = neweta;
        }

        debug!("Eta update");
        debug!("  primopt={}", self.primopt);
        debug!("  eta    ={}", self.eta);

        return updated_eta;
    }


    // Compute the new candidate point.

    //
    // This consists of two steps:

    //
    // 1. the new point is computed as $-\tfrac{1}{u}\bar{g}$, where $\bar{g}$
    //    is the aggregated minorant
    // 2. the multipliers $\eta$ are updated

    //
    // In other words, this function computes the new candidate
    // 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.

            norm2 += x * x;
        }
        return norm2;
    }
}


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

    fn set_num_subproblems(&mut self, n: usize) -> Result<()> {
        self.master.set_num_subproblems(n).map_err(|err| Error::Solver(Box::new(err)))
    }

    fn set_vars(&mut self, n: usize, lb: Option<DVector>, ub: Option<DVector>) {
        assert!(lb.as_ref().map(|x| x.len()).unwrap_or(n) == n);
        assert!(ub.as_ref().map(|x| x.len()).unwrap_or(n) == n);
        self.lb = lb.unwrap_or_else(|| dvec![NEG_INFINITY; n]);
        self.ub = ub.unwrap_or_else(|| dvec![INFINITY; n]);
    }

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

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<Self::MinorantIndex> {
        self.master.add_minorant(fidx, minorant).map_err(|err| Error::Solver(Box::new(err)))
    }

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

    fn set_weight(&mut self, weight: Real) {
        self.master.set_weight(weight);
    }

    fn add_vars(&mut self, bounds: &[(Real, Real)], extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector) {
        if !bounds.is_empty() {
            self.lb.extend(bounds.iter().map(|x| x.0));
            self.ub.extend(bounds.iter().map(|x| x.1));
            self.eta.resize(self.lb.len(), 0.0);
            self.need_new_candidate = true;
            self.master.add_vars(bounds.len(), extend_subgradient)
        }

            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");
        debug!("  cnt_update={}", cnt_updates);
        debug!("  primopt={}", self.primopt);
        debug!("  dualopt={}", self.master.dualopt());
        debug!("  etaopt={}", self.eta);
        debug!("  primoptval={}", self.primoptval);

        Ok(())
    }

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector)> {
        self.master.aggregate(fidx, mins).map_err(|err| Error::Solver(Box::new(err)))
    }

    fn get_primopt(&self) -> DVector {
        self.primopt.clone()
    }

    fn get_primoptval(&self) -> Real {
        self.primoptval
    }

    fn get_dualoptnorm2(&self) -> Real {
        self.dualoptnorm2
    }

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

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.need_new_candidate = true;
        self.master.move_center(alpha, d);
        for i in 0..self.primopt.len() {

// 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
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.

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

//! Master problem implementation using CPLEX.

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

use cplex;

        }
    }
}

pub type Result<T> = result::Result<T, Error>;

pub struct CplexMaster {
    lp: *mut CPXlp,

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

    /// List of free minorant indices.
    freeinds: Vec<usize>,

    /// List of minorant indices to be updated.
    updateinds: Vec<usize>,

    /// Mapping minorant to index.
    min2index: Vec<Vec<usize>>,

    /// Mapping index to minorant.
    index2min: Vec<(usize, usize)>,

    /// The quadratic term.
    qterm: Vec<DVector>,

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

    }

    fn set_weight(&mut self, weight: Real) {
        assert!(weight > 0.0);
        self.weight = weight;
    }

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

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<usize> {
        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;

        }
    }

    fn add_vars(&mut self, nvars: usize, extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector) {
        debug_assert!(!self.minorants[0].is_empty());
        let noldvars = self.minorants[0][0].linear.len();
        let nnewvars = noldvars + nvars;
        let newvars = (noldvars..nnewvars).collect::<Vec<_>>();
        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, &newvars);
                    m.linear.extend_from_slice(&new_subg);
                }
            }
        }
        // update qterm
        if self.force_update {
            return;
        }
        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 = (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];
                        }
                    }
                }
            }
        }

        // WORST CASE: DO THIS
        // self.force_update = true;
    }

    #[allow(unused_variables)]
    fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<()> {
        if self.force_update || !self.updateinds.is_empty() {
            try!(self.init_qp());
        }

        let nvars = unsafe { CPXgetnumcols(env(), self.lp) as usize };
        if nvars == 0 {
            return Err(Error::NoMinorants);
        }
        // update linear costs
        {
            let mut c = Vec::with_capacity(nvars);
            let mut inds = Vec::with_capacity(nvars);
            for mins in self.minorants.iter() {
                for m in mins {
                    inds.push(c.len() as c_int);
                    c.push(-m.constant * self.weight - m.linear.dot(eta));
                }
            }
            trycpx!(CPXchgobj(env(), self.lp, nvars as c_int, inds.as_ptr(), c.as_ptr()));
        }

        trycpx!(CPXqpopt(env(), self.lp));
        let mut sol = vec![0.0; nvars];
        trycpx!(CPXgetx(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];

        &self.opt_minorant.linear
    }

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

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

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

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector)> {
        assert!(mins.len() > 0, "No minorants specified to be aggregated");

        if mins.len() == 1 {
            return Ok((mins[0], dvec![1.0]));
        }

        // scale coefficients
        let mut sum_coeffs = 0.0;
        for &i in mins {
            sum_coeffs += self.opt_mults[fidx][self.index2min[i].1];
        }
        let aggr_coeffs = if sum_coeffs != 0.0 {
            mins.iter()
                .map(|&i| self.opt_mults[fidx][self.index2min[i].1] / sum_coeffs)
                .collect::<DVector>()
        } else {
            dvec![0.0; mins.len()]
        };

        // compute aggregated diagonal term
        let mut aggr_diag = 0.0;
        for (idx_i, &i) in mins.iter().enumerate() {

        {
            let mut aggr_mins = Vec::with_capacity(mins.len());

            for &i in mins {
                let (min_fidx, min_idx) = self.index2min[i];
                debug_assert!(min_fidx == fidx);

                let m = self.minorants[fidx].swap_remove(min_idx);
                let idx = self.min2index[fidx].swap_remove(min_idx);
                self.opt_mults[fidx].swap_remove(min_idx);
                self.freeinds.push(idx);
                debug_assert!(idx == i);

                aggr_mins.push(m);

                // update index2min table for moved minorant
                if min_idx < self.minorants[fidx].len() {
                    self.index2min[self.min2index[fidx][min_idx]].1 = min_idx;
                }
            }
            aggr.combine_all(&aggr_coeffs, &aggr_mins);
        }

        // save aggregated minorant
        let aggr_idx = self.freeinds.pop().unwrap();
        self.minorants[fidx].push(aggr);
        self.opt_mults[fidx].push(sum_coeffs);
        self.min2index[fidx].push(aggr_idx);
        self.index2min[aggr_idx] = (fidx, self.minorants[fidx].len() - 1);

        // update qterm
        for fidx_i in 0..self.minorants.len() {
            for idx_i in 0..self.minorants[fidx_i].len() {
                let i = self.min2index[fidx_i][idx_i];
                self.qterm[i][aggr_idx] = aggr_qterm[i];
                self.qterm[aggr_idx][i] = aggr_qterm[i];

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

            trycpx!(CPXaddrows(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);
            }

            // the global indices for each minorant in order
            let mut activeinds = vec![];
            for (fidx, mins_i) in self.minorants.iter().enumerate() {
                for (i, m_i) in mins_i.iter().enumerate() {
                    let idx_i = self.min2index[fidx][i];
                    activeinds.push(idx_i);

                        debug_assert!((self.qterm[idx_i][idx_j] - m_i.linear.dot(&m_j.linear)).abs() < 1e-6);
                    }
                }
            }

            // main diagonal plus small identity to ensure Q being semi-definite
            let mut maxq = 0.0;
            for &i in &activeinds {
                maxq = f64::max(maxq, self.qterm[i][i])
            }
            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!(CPXchgqpcoef(env(),
                                             self.lp,
                                             i as c_int,
                                             j as c_int,
                                             self.qterm[idx_i][idx_j]));
                    } else {
                        trycpx!(CPXchgqpcoef(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(())
    }
}

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

use std::error;
use std::result;

quick_error! {

/// Result type for master problems.
pub type Result<T> = result::Result<T, Error>;

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

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

    /// Set the lower and upper bounds of the variables.
    fn set_vars(&mut self, nvars: usize, lb: Option<DVector>, ub: Option<DVector>);

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

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

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

    /// Set the maximal number of inner iterations.
    fn set_max_updates(&mut self, max_updates: usize);

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

    /// Add some variables with bounds.
    fn add_vars(&mut self, bounds: &[(Real, Real)], extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector);

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

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

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

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

// 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
// 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 {Real, DVector, Minorant};
use master::UnconstrainedMasterProblem;

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

            weight: 1.0,
            minorants: vec![],
            opt_mult: dvec![],
            opt_minorant: Minorant::default(),
        })
    }

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

    fn set_num_subproblems(&mut self, n: usize) -> Result<()> {
        if n != 1 {
            Err(Error::NumSubproblems(n))
        } else {
            Ok(())
        }

    }

    fn num_minorants(&self, fidx: usize) -> usize {
        assert!(fidx == 0);
        self.minorants.len()
    }

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<usize> {
        assert!(fidx == 0);
        if self.minorants.len() >= 2 {
            return Err(Error::MaxMinorants);
        }
        self.minorants.push(minorant);
        self.opt_mult.push(0.0);
        Ok(self.minorants.len() - 1)
    }

    fn add_vars(&mut self, nvars: usize, extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> DVector) {
        if !self.minorants.is_empty() {
            let noldvars = self.minorants[0].linear.len();
            let newvars = (noldvars..noldvars + nvars).collect::<Vec<_>>();
            for (i, m) in self.minorants.iter_mut().enumerate() {
                let new_subg = extend_subgradient(0, i, &newvars);
                m.linear.extend_from_slice(&new_subg);
            }
        }
    }

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

        if self.minorants.len() == 2 {
            let xx = self.minorants[0].linear.dot(&self.minorants[0].linear);
            let yy = self.minorants[1].linear.dot(&self.minorants[1].linear);
            let xy = self.minorants[0].linear.dot(&self.minorants[1].linear);
            let xeta = self.minorants[0].linear.dot(eta);

            self.opt_mult[1] = alpha2;
            self.opt_minorant = self.minorants[0].combine(1.0 - alpha2, alpha2, &self.minorants[1]);
        } else if self.minorants.len() == 1 {
            self.opt_minorant = self.minorants[0].clone();
            self.opt_mult.resize(1, 1.0);
            self.opt_mult[0] = 1.0;
        } else {
            return Err(Error::NoMinorants);
        }

        debug!("Unrestricted");
        debug!("  opt_minorant={}", self.opt_minorant);
        if self.opt_mult.len() == 2 {
            debug!("  opt_mult={}", self.opt_mult);
        }

        Ok(())
    }

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

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

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

    fn eval_model(&self, y: &DVector) -> Real {
        let mut result = NEG_INFINITY;
        for m in &self.minorants {
            result = result.max(m.eval(y));

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

//! Bundle master problem solver.
//!
//! This module contains solvers for the bundle master problem, i.e.
//! for solving convex optimization problems of the form
//!
//! \\[ \min \left\\{ \hat{f}(d) + \frac{w}{2} \\|d\\|\^2 \colon d \in [l,u] \right\\}, \\]

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

use std::error;
use std::result;

/**

 * 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 {
    /// Error type.
    type Error: error::Error + 'static;

    /// Unique index for a minorant.
    type MinorantIndex: Copy + Eq;

    /// Return a new instance of the unconstrained master problem.
    fn new() -> result::Result<Self, 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::Result<(), Self::Error>;

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

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

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

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

    /// Aggregate the given minorants according to the current solution.
    ///
    /// The (indices of the) minorants to be aggregated get invalid

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

//! Multi-commodity min-cost-flow subproblems.

mod solver;
pub use mcf::solver::Solver;

mod problem;

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

#[allow(dead_code)]

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

use std::fs::File;

        }
    }
}

pub type Result<T> = result::Result<T, Error>;

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

#[derive(Clone, Copy, Debug)]
struct Elem {
    ind: usize,
    val: Real,
}

pub struct MMCFProblem {
    pub multimodel: bool,

    nets: Vec<mcf::Solver>,
    lhs: Vec<Vec<Vec<Elem>>>,
    rhs: DVector,
    rhsval: Real,
    cbase: Vec<DVector>,
    c: Vec<DVector>,
}

impl MMCFProblem {
    pub fn read_mnetgen(basename: &str) -> Result<MMCFProblem> {
        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(Error::Format(format!("Expected 4 numbers in {}.nod, but got {}",
                                             basename,
                                             fnod.len())));
        }

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

        // read nodes
        let mut nets = Vec::with_capacity(ncom);
        for _ in 0..ncom {
            nets.push(try!(mcf::Solver::new(nnodes)))
        }
        {
            let mut f = try!(File::open(&format!("{}.sup", basename)));
            buffer.clear();
            try!(f.read_to_string(&mut buffer));
        }
        for line in buffer.lines() {
            let mut data = line.split_whitespace();
            let node = try!(data.next().unwrap().parse::<usize>());
            let com = try!(data.next().unwrap().parse::<usize>());
            let supply = try!(data.next().unwrap().parse::<Real>());
            try!(nets[com - 1].set_balance(node - 1, supply));
        }

        // read arcs
        let mut arcmap = vec![vec![]; ncom];
        let mut cbase = vec![dvec![]; ncom];

        // lhs nonzeros
        let mut lhsidx = vec![vec![vec![]; ncom]; ncaps];
        {
            let mut f = try!(File::open(&format!("{}.arc", basename)));
            buffer.clear();
            try!(f.read_to_string(&mut buffer));
        }
        for line in buffer.lines() {
            let mut data = line.split_whitespace();
            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
        {
            let mut f = try!(File::open(&format!("{}.mut", basename)));
            buffer.clear();
            try!(f.read_to_string(&mut buffer));
        }
        let mut rhs = dvec![0.0; ncaps];
        for line in buffer.lines() {
            let mut data = line.split_whitespace();
            let mt = try!(data.next().unwrap().parse::<usize>()) - 1;
            let cap = try!(data.next().unwrap().parse::<Real>());
            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: &Vec<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];
                }
            }
            sum
        }
    }

    /// Aggregate primal vectors.
    pub fn aggregate_primals_ref(&self, primals: &[(Real, &Vec<DVector>)]) -> Vec<DVector> {
        let mut aggr = primals[0]
            .1
            .iter()
            .map(|x| {
                let mut r = dvec![];
                r.scal(primals[0].0, x);
                r
            })
            .collect::<Vec<_>>();

        for &(alpha, primal) in &primals[1..] {
            for (j, x) in primal.iter().enumerate() {
                aggr[j].add_scaled(alpha, x);
            }
        }

        aggr
    }
}

impl<'a> FirstOrderProblem<'a> for MMCFProblem {
    type Error = Error;

    type Primal = Vec<DVector>;

    type EvalResult = SimpleEvaluation<Vec<DVector>>;

    fn num_variables(&self) -> usize {
        self.lhs.len()
    }

    fn lower_bounds(&self) -> Option<Vector> {
        Some(Vector::new_sparse(self.lhs.len(), &[], &[]))
    }

    fn upper_bounds(&self) -> Option<Vector> {
        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: &DVector, nullstep_bound: Real, relprec: Real) -> result::Result<Self::EvalResult, Self::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 in 0..self.lhs.len() {

        debug!("y={}", y);
        for i in 0..self.nets.len() {
            debug!("c[{}]={}", i, self.c[i]);
            try!(self.nets[i].set_objective(&self.c[i]));
        }

        // solve subproblems
        for (i, net) in self.nets.iter_mut().enumerate() {
            try!(net.solve());
            debug!("c[{}]={}", i, try!(net.objective()));
        }

        // compute minorant
        if self.multimodel {
            let objective;

                for elem in &self.lhs[i][fidx] {
                    subg[i] -= elem.val * sol[elem.ind];
                }
            }

            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 in 0..self.lhs.len() {
                for (fidx, flhs) in self.lhs[i].iter().enumerate() {
                    for elem in flhs {
                        subg[i] -= elem.val * sols[fidx][elem.ind];
                    }
                }
            }

            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<_>>())
    }
}

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

use cplex;
use cplex::*;

use std::ptr;

        }
    }
}

pub type Result<T> = result::Result<T, Error>;

pub struct Solver {
    net: *mut CPXnet,
    logfile: *mut CPXfile,
}


impl Drop for Solver {
    fn drop(&mut self) {
        unsafe {
            CPXNETfreeprob(cplex::env(), &mut self.net);
            CPXfclose(self.logfile);
        }
    }
}


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

        unsafe {
            loop {
                logfile = CPXfopen(CString::new("mcf.cpxlog").unwrap().as_ptr(),
                                   CString::new("w").unwrap().as_ptr());
                if logfile == ptr::null_mut() {
                    return Err(Error::Cplex(CplexError {
                        code: 0,
                        msg: "Can't open log-file".to_string(),
                    }));
                }
                status = CPXsetlogfile(env(), logfile);
                if status != 0 {
                    break;
                }

                net = CPXNETcreateprob(env(), &mut status, CString::new("mcf").unwrap().as_ptr());
                if status != 0 {
                    break;
                }
                status = CPXNETaddnodes(env(), net, nnodes as c_int, ptr::null(), ptr::null());
                if status != 0 {
                    break;
                }
                status = CPXNETchgobjsen(env(), net, CPX_MIN);
                if status != 0 {
                    break;
                }
                break;
            }

            if status != 0 {
                let msg = CString::new(vec![0; CPXMSGBUFSIZE]).unwrap().into_raw();
                CPXgeterrorstring(env(), status, msg);
                CPXNETfreeprob(env(), &mut net);
                CPXfclose(logfile);
                return Err(Error::Cplex(CplexError {
                    code: status,
                    msg: CString::from_raw(msg).to_string_lossy().into_owned(),
                }));
            }
        }

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

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

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

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

    pub fn set_objective(&mut self, obj: &DVector) -> Result<()> {
        let inds = (0..obj.len() as c_int).collect::<Vec<_>>();
        Ok(trycpx!(CPXNETchgobj(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<()> {
        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!(CPXNETaddarcs(env(),
                                 self.net,
                                 1,
                                 &f,
                                 &t,
                                 ptr::null(),
                                 &u,
                                 &c,
                                 &cname as *const *const c_char)))
    }

    pub fn solve(&mut self) -> Result<()> {
        Ok(trycpx!(CPXNETprimopt(env(), self.net)))
    }

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

    pub fn get_solution(&self) -> Result<DVector> {
        let mut sol = dvec![0.0; self.num_arcs()];
        let mut stat: c_int = 0;
        let mut objval: c_double = 0.0;
        trycpx!(CPXNETsolution(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()));
        return Ok(sol);

// 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
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.

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

//! A linear minorant.

use {Real, DVector};

use std::fmt;

 *   \rangle + c \\]
 *
 * such that $l(x) \le f(x)$ for all $x \in \mathbb{R}\^n$.
 */
#[derive(Clone, Debug, Default)]
pub struct Minorant {
    /// The constant term.
    pub constant: Real,

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


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

     */
    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]) {

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

//! The main bundle method solver.

use {Real, DVector};
use {FirstOrderProblem, Update, Evaluation, HKWeighter};

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

 * Captures the current state of the bundle method during the run of
 * the algorithm. This state is passed to certain callbacks like
 * Terminator or Weighter so that they can compute their result
 * depending on the state.
 */
pub struct BundleState<'a> {
    /// Current center of stability.
    pub cur_y: &'a DVector,

    /// Function value in current center.
    pub cur_val: Real,

    /// Current candidate, point of last evaluation.
    pub nxt_y: &'a DVector,

    /// Function value in candidate.
    pub nxt_val: Real,

    /// Model value in candidate.
    pub nxt_mod: Real,

    /// Cut value of new subgradient in current center.
    pub new_cutval: Real,

    /// The current aggregated subgradient norm.
    pub sgnorm: Real,

    /// The expected progress of the current model.
    pub expected_progress: Real,

    /// Currently used weight of quadratic term.
    pub weight: Real,

    /**
     * The type of the current step.
     *
     * If the current step is Step::Term, the weighter should be reset.
     */
    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,

 * Termination predicate.
 *
 * Given the current state of the bundle method, this function returns
 * whether the solution process should be stopped.
 */
pub trait Terminator {
    /// Return true if the method should stop.
    fn terminate(&mut self, state: &BundleState, params: &SolverParams) -> bool;
}


/**
 * Terminates if expected progress is small enough.
 */
pub struct StandardTerminator {
    pub termination_precision: Real,
}

impl Terminator for StandardTerminator {
    #[allow(unused_variables)]
    fn terminate(&mut self, state: &BundleState, params: &SolverParams) -> bool {
        assert!(self.termination_precision >= 0.0);
        state.expected_progress <= self.termination_precision * (state.cur_val.abs() + 1.0)
    }
}

/**
 * Bundle weight controller.
 *
 * Given the current state of the bundle method, this function determines the
 * weight factor of the quadratic term for the next iteration.
 */
pub trait Weighter {
    /// Return the new weight of the quadratic term.
    fn weight(&mut self, state: &BundleState, params: &SolverParams) -> Real;
}

/// Parameters for tuning the solver.
#[derive(Clone, Debug)]
pub struct SolverParams {
    /// Maximal individual bundle size.
    pub max_bundle_size: usize,

    /**
     * Factor for doing a descent step.
     *
     * If the proportion of actual decrease to predicted decrease is
     * at least that high, a descent step will be done.
     *

     * compute a bound for the function oracle, that guarantees a null
     * step. If the function is evaluated by some iterative method that ensures
     * an objective value that is at least as large as this bound, the
     * oracle can stop returning an appropriate $\varepsilon$-subgradient.
     *
     * Must be in (0, acceptance_factor).
     */
    pub nullstep_factor: Real,

    /// Minimal allowed bundle weight. Must be > 0 and < max_weight.
    pub min_weight: Real,

    /// Maximal allowed bundle weight. Must be > min_weight,
    pub max_weight: Real,

    /**
     * Maximal number of updates of box multipliers.
     *
     * This is the maximal number of iterations for updating the box
     * multipliers when solving the master problem with box
     * constraints. This is a technical parameter that should probably
     * never be changed. If you experience an unexpectedly high number
     * of inner iterations, consider removing/fixing the corresponding
     * variables.
     */
    pub max_updates: usize,
}

impl SolverParams {
    /// Verify that all parameters are valid.
    fn check(&self) -> Result<()> {
        if self.max_bundle_size < 2 {
            Err(Error::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(Error::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(Error::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(Error::Parameter(format!("min_weight must be in > 0 (got: {})", self.min_weight)))
        } else if self.max_weight < self.min_weight {
            Err(Error::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(Error::Parameter(format!("max_updates must be in > 0 (got: {})", self.max_updates)))
        } else {
            Ok(())
        }
    }
}

    /// Primal associated with this minorant.
    primal: Option<Pr>,
}

/// Information about the last iteration.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum IterationInfo {
    NewMinorantTooHigh { new: Real, old: Real },
    UpperBoundNullStep,
    ShallowCut,
}


/// State information for the update callback.
pub struct UpdateState<'a, Pr: 'a> {
    /// Current model minorants.
    minorants: &'a [Vec<MinorantInfo<Pr>>],
    /// The last step type.
    pub step: Step,
    /// Iteration information.
    pub iteration_info: &'a [IterationInfo],
    /// The current candidate. If the step was a descent step, this is
    /// the new center.
    pub nxt_y: &'a DVector,
    /// The center. IF the step was a descent step, this is the old
    /// center.
    pub cur_y: &'a DVector,
}

impl<'a, Pr: 'a> UpdateState<'a, Pr> {
    pub fn aggregated_primals(&self, subproblem: usize) -> Vec<(Real, &Pr)> {
        self.minorants[subproblem]
            .iter()
            .map(|m| (m.multiplier, m.primal.as_ref().unwrap()))
            .collect()
    }

    /// Return the last primal for a given subproblem.
    ///
    /// This is the last primal generated by the oracle.
    pub fn last_primal(&self, fidx: usize) -> Option<&Pr> {
        self.minorants[fidx].last().and_then(|m| m.primal.as_ref())
    }
}

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

    /// Termination predicate.
    pub terminator: Box<Terminator>,

    /// Weighter heuristic.
    pub weighter: Box<Weighter>,

    /// Current center of stability.
    cur_y: DVector,

    /// Function value in current point.
    cur_val: Real,

    /// Model value in current point.
    cur_mod: Real,

    /// Vector of subproblem function values in current point.
    cur_vals: DVector,

    cnt_null: usize,

    /**
     * Time when the solution process started.
     *
     * This is actually the time of the last call to `Solver::init`.
     */
    start_time: Instant,

    /// The master problem.
    master: Box<MasterProblem<MinorantIndex = usize>>,

    /// The active minorant indices for each subproblem.
    minorants: Vec<Vec<MinorantInfo<Pr>>>,

    /// 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>> {
        Ok(Solver {
            problem: problem,
            params: params,
            terminator: Box::new(StandardTerminator { termination_precision: 1e-3 }),
            weighter: Box::new(HKWeighter::new()),
            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: match BoxedMasterProblem::<MinimalMaster>::new() {
                Ok(master) => Box::new(master),
                Err(err) => return Err(Error::Master(Box::new(err))),
            },
            minorants: vec![],
            iterinfos: vec![],
        })
    }

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

    /**
     * Set the first order problem description associated with this
     * solver.
     *
     * 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 set_problem(&mut self, problem: P) {
        self.problem = problem;
    }

    /// Returns a reference to the solver's current problem.
    pub fn problem(&self) -> &P {
        &self.problem
    }

        }

        let lb = self.problem.lower_bounds().map(|x| x.to_dense());
        let ub = self.problem.upper_bounds().map(|x| x.to_dense());
        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(Error::InvalidBounds(lb_i, 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;
                self.cur_y[i] = ub_i;
            }

    /// Solve the problem.
    pub fn solve(&mut self) -> Result<()> {
        try!(self.init());
        for _ in 0..100000 {
            let mut term = try!(self.step());
            let changed = try!(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 {
                break;
            }
        }
        Ok(())
    }

            let state = UpdateState {
                minorants: &self.minorants,
                step: term,
                iteration_info: &self.iterinfos,
                // this is a dirty trick: when updating the center, we
                // simply swapped the `cur_*` fields with the `nxt_*`
                // fields
                cur_y: if term == Step::Descent {
                    &self.nxt_y
                } else {
                    &self.cur_y
                },
                nxt_y: if term == Step::Descent {
                    &self.cur_y
                } else {
                    &self.nxt_y
                },
            };
            match self.problem.update(&state) {
                Ok(updates) => updates,
                Err(err) => return Err(Error::Update(Box::new(err))),
            }
        };

        let mut newvars = Vec::with_capacity(updates.len());
        for u in updates {
            match u {
                Update::AddVariable { lower, upper } => {
                    if lower > upper {
                        return Err(Error::InvalidBounds(lower, upper));
                    }
                    let value = if lower > 0.0 {
                        lower
                    } else if upper < 0.0 {
                        upper
                    } else {
                        0.0
                    };
                    newvars.push((lower - value, upper - value));

                }
                Update::AddVariableValue { lower, upper, value } => {
                    if lower > upper {
                        return Err(Error::InvalidBounds(lower, upper));
                    }
                    if value < lower || value > upper {
                        return Err(Error::ViolatedBounds(lower, upper, value));
                    }
                    newvars.push((lower - value, upper - value));

                }
            }
        }

        if !newvars.is_empty() {
            let mut problem = &mut self.problem;
            let minorants = &self.minorants;
            self.master.add_vars(&newvars,
                                 &mut move |fidx, minidx, vars| {
                                     problem.extend_subgradient(minorants[fidx][minidx].primal.as_ref().unwrap(), vars)
                                         .unwrap()
                                 });
            let newn = self.cur_y.len() + newvars.len();
            self.cur_y.resize(newn, 0.0);
            self.nxt_d.resize(newn, 0.0);
            self.nxt_y.resize(newn, 0.0);
            Ok(true)
        } else {
            Ok(false)
        }
    }

    /// Return the current aggregated primal information for a subproblem.
    ///
    /// This function returns all currently used minorants $x_i$ along
    /// with their coefficients $\alpha_i$. The aggregated primal can
    /// be computed by combining the minorants $\bar{x} =
    /// \sum_{i=1}\^m \alpha_i x_i$.
    pub fn aggregated_primals(&self, subproblem: usize) -> Vec<(Real, &Pr)> {
        self.minorants[subproblem]
            .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() / 10000000,
              self.cnt_descent,
              self.cnt_descent + self.cnt_null,
              self.master.cnt_updates(),

        };

        let lb = self.problem.lower_bounds().map(|v| v.to_dense());
        let ub = self.problem.upper_bounds().map(|v| v.to_dense());

        if let Some(ref x) = lb {
            if x.len() != self.problem.num_variables() {
                return Err(Error::Dimension("Dimension of lower bounds does not match number of \
                                             variables"));
            }
        }

        try!(self.master.set_num_subproblems(m));
        self.master.set_vars(self.problem.num_variables(), lb, ub);
        self.master.set_max_updates(self.params.max_updates);

        self.minorants = Vec::with_capacity(m);
        for _ in 0..m {
            self.minorants.push(vec![]);
        }

        self.cur_val = 0.0;
        for i in 0..m {
            let result = match self.problem.evaluate(i, &self.cur_y, INFINITY, 0.0) {
                Ok(r) => r,
                Err(err) => return Err(Error::Eval(Box::new(err))),
            };

    /// Reduce size of bundle.
    fn compress_bundle(&mut self) -> Result<()> {
        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) = try!(self.master.aggregate(i, &aggr_mins));
                // 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.

        try!(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
        };

        try!(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 mut minorants = result.into_iter();
            let mut nxt_minorant;
            let nxt_primal;
            match minorants.next() {
                Some((m, p)) => {
                    nxt_minorant = m;
                    nxt_primal = p;

                }
                None => return Err(Error::NoMinorant),
            }
            let fun_lb = nxt_minorant.constant;

            nxt_lb += fun_lb;
            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: try!(self.master.add_minorant(fidx, nxt_minorant)),
                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,
            });
        }

        self.nxt_val = nxt_ub;

        // check for potential problems with relative precision of all kinds
        if nxt_lb <= descent_bnd {
            // lower bound gives descent step

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

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

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

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

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

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

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

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

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

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


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

            self[i] = 0.0;
        }
    }

    /// Set self = factor * y.
    pub fn scal(&mut self, factor: Real, y: &DVector) {
        self.resize(y.len(), 0.0);
        for (i, x) in y.iter().enumerate() {
            self[i] = factor * x;
        }
    }

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

        self.resize(x.len(), 0.0);
        for i in 0..x.len() {
            self[i] = x[i] + y[i];
        }
    }

    /// Add two vectors and store result in this vector.
    pub fn add_scaled(&mut self, alpha: Real, y: &DVector) {
        assert!(self.len() == y.len());
        for i in 0..self.len() {
            self[i] += alpha * y[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) {
        if self.len() < y.len() {
            self.resize(y.len(), 0.0);
        }
        for i in 0..y.len() {
            self[i] += alpha * y[i];
        }
    }


    /// Combines this vector with another vector.
    pub fn combine(&self, self_factor: Real, other_factor: Real, other: &DVector) -> DVector {
        assert!(self.len() == other.len());
        let mut result = DVector(Vec::with_capacity(self.len()));
        for i in 0..self.len() {
            result.push(self_factor * self[i] + other_factor * other[i]);
        }
        result
    }


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


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

        if indices.len() == 0 {
            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 values[i] != 0.0 {
                    if indices[i] != last_idx {
                        elems.push((indices[i], values[i]));
                        last_idx = indices[i];
                    } else {
                        elems.last_mut().unwrap().1 += values[i];
                        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 {
                    v[i] = x;
                }
                DVector(v)
            }
        }
    }
}