RsBundle  Check-in [df7544e94f]

[package]
name = "bundle"
version = "0.3.1"
authors = ["Frank Fischer <frank-fischer@shadow-soft.de>"]

[dependencies]
libc = "^0.2.6"
quick-error = "^1.1.0"

log = "^0.3.6"
const-cstr = "^0.2.1"
cplex-sys = "^0.2"

[dev-dependencies]
env_logger = "^0.4.1"

 */

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


use bundle::{Real, DVector, Minorant, SimpleEvaluation, FirstOrderProblem, Solver, SolverParams};
use std::io;


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 Error = io::Error;
    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, Self::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]);

//

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

use {Real, 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.
 *

}

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

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

        let mut primals = primals;
        primals.pop().unwrap().1
    }

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

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

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

//! Proximal bundle method implementation.

#[macro_use]
extern crate quick_error;


#[macro_use]
extern crate const_cstr;

#[macro_use]
extern crate log;

/// Type used for real numbers throughout the library.

//
// 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! {
    /// A master problem error.
    #[derive(Debug)]
    pub enum Error {
        Solver(err: Box<error::Error>) {
            cause(&**err)
            description(err.description())
            display("Master problem solver error: {}", err)
        }
    }
}


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

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

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

    /// Aggregate the given minorants according to the current
    /// solution.
    ///
    /// The (indices of the) minorants to be aggregated get invalid
    /// after this operation. The function returns the index of the
    /// aggregated minorant along with the coefficients of the convex
    /// combination. The index of the new aggregated minorant might or
    /// might not be one of indices of the original minorants.
    ///
    /// # Error The indices of the minorants `mins` must belong to
    /// subproblem `fidx`.
    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector)>;

    /// Return the (primal) optimal solution $\\|d\^*\\|$.
    fn get_primopt(&self) -> DVector;

    /// Return the value of the linear model in the optimal solution.
    ///
    /// This is the term $\langle g\^*, d\^* \rangle$ where $g^\*$ is

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


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


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

    /// 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) {
        assert!(max_updates > 0);
        self.max_updates = max_updates;
    }

    }
}


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_eq!(lb.as_ref().map(|x| x.len()).unwrap_or(n), n);
        assert_eq!(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) {

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

    #[cfg_attr(feature="cargo-clippy", allow(cyclomatic_complexity))]
    fn solve(&mut self, center_value: Real) -> Result<()> {
        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
            if let Err(err) = self.master.solve(&self.eta, center_value, old_augval, 1e-3) {
                return Err(Error::Solver(Box::new(err)));
            }

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

        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 {

//! Master problem implementation using CPLEX.

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

use cplex_sys as cpx;

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



quick_error! {
    /// A solver error.
    #[derive(Debug)]
    pub enum Error {
        Cplex(err: cpx::Error) {
            cause(err)
                description(err.description())
                display("{}", err)
                from()
        }

        NoMinorants {
            description("No minorants")
                display("Solver Error: no minorants when solving the master problem")
        }
    }
}

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

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

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

    fn drop(&mut self) {
        unsafe { cpx::freeprob(cpx::env(), &mut self.lp) };
    }
}


impl UnconstrainedMasterProblem for CplexMaster {
    type Error = Error;

    type MinorantIndex = usize;

    fn new() -> Result<CplexMaster> {
        Ok(CplexMaster {
            lp: ptr::null_mut(),
            force_update: true,
            freeinds: vec![],
            updateinds: vec![],
            min2index: vec![],
            index2min: vec![],
            qterm: vec![],
            weight: 1.0,
            minorants: vec![],
            opt_mults: vec![],
            opt_minorant: Minorant::default(),
        })
    }

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

    fn set_num_subproblems(&mut self, n: usize) -> Result<()> {
        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.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();

            }
        }

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

    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 { cpx::getnumcols(cpx::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 {
                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.is_empty(), "No minorants specified to be aggregated");

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

        // scale coefficients

            }
        }
    }
}


impl CplexMaster {
    fn init_qp(&mut self) -> Result<()> {
        if !self.lp.is_null() {
            trycpx!(cpx::freeprob(cpx::env(), &mut self.lp));
        }
        trycpx!({
            let mut status = 0;
            self.lp = cpx::createprob(cpx::env(), &mut status, const_cstr!("mcf").as_ptr());
            status

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




quick_error! {
    /// A solver error.
    #[derive(Debug)]
    pub enum Error {
        NumSubproblems(n: usize) {
            description("Too many subproblems")
            display("Solver Error: too many subproblems (got: {} must be <= 2)", n)
        }

        MaxMinorants {
            description("Too many minorants")
            display("Solver Error: the minimal master problem allows at most two minorants")
        }

        NoMinorants {
            description("No minorants")
            display("Solver Error: no minorants when solving the master problem")

        }
    }
}

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

/**
 * 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
 * is that this model can be solved explicitely and very quickly, but

    opt_mult: DVector,
    /// Optimal aggregated minorant.
    opt_minorant: Minorant,
}


impl UnconstrainedMasterProblem for MinimalMaster {
    type Error = Error;

    type MinorantIndex = usize;

    fn new() -> Result<MinimalMaster> {
        Ok(MinimalMaster {
            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 weight(&self) -> Real {
        self.weight
    }

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

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

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

                }
                m.linear.extend_from_slice(&new_subg[changed.len()..]);
            }
        }
    }

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

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

        let mut result = NEG_INFINITY;
        for m in &self.minorants {
            result = result.max(m.eval(y));
        }
        result
    }

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector)> {
        assert_eq!(fidx, 0);
        if mins.len() == 2 {
            debug!("Aggregate");
            debug!("  {} * {}", self.opt_mult[0], self.minorants[0]);
            debug!("  {} * {}", self.opt_mult[1], self.minorants[1]);
            self.minorants[0] = self.minorants[0].combine(self.opt_mult[0], self.opt_mult[1], &self.minorants[1]);
            self.minorants.truncate(1);

//! * changing the weight parameter $w$,
//! * modifying $\hat{f}$ by adding or removing linear functions $\ell_i$,
//! * moving the center of the linear functions $\ell_i$ (and the
//!   bounds), i.e. replacing $\hat{f}$ by $d \mapsto \hat{f}(d -
//!   \hat{d})$ for some given $\hat{d} \in \mathbb{R}\^n$.

mod base;
pub use self::base::{MasterProblem, Error, Result};

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

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

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

/**
 * Trait for master problems without box constraints.
 *
 * Implementors of this trait are supposed to solve quadratic
 * optimization problems of the form
 *

 * \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 {
    /// 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 weight.
    fn weight(&self) -> Real;

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

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

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

    /// Solve the master problem.
    fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> 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;

    /// after this operation. The function returns the index of the
    /// aggregated minorant along with the coefficients of the convex
    /// combination. The index of the new aggregated minorant might or
    /// might not be one of indices of the original minorants.
    ///
    /// # Error
    /// The indices of the minorants `mins` must belong to subproblem `fidx`.
    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector), Self::Error>;

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

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

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

use std::fs::File;
use std::io::{self, Read};

use std::result;
use std::num::{ParseIntError, ParseFloatError};
use std::f64::INFINITY;

quick_error! {
    /// A solver error.
    #[derive(Debug)]
    pub enum Error {
        Io(err: io::Error) {
            cause(err)
            description(err.description())
            display("Io Error: {}", err)
            from()
        }

        Solver(err: mcf::solver::Error) {
            cause(err)
            description(err.description())
            display("Solver Error: {}", err)
            from()
        }

        Format(msg: String) {
            description("Format error")
            display("Format error: {}", msg)
            from(err: ParseIntError) -> (format!("{}", err))
            from(err: ParseFloatError) -> (format!("{}", err))
        }
    }
}

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

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

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

        }

        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 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::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());
        }

use {Real, DVector};

use cplex_sys as cpx;

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

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

quick_error! {
    #[derive(Debug)]
    pub enum Error {
        Cplex(err: cpx::Error) {
            cause(err)
                description(err.description())
                display("{}", err)
                from()
        }
    }
}

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

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


impl Drop for Solver {
    fn drop(&mut self) {
        unsafe {
            cpx::NETfreeprob(cpx::env(), &mut self.net);
            cpx::fclose(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 {
            #[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(Error::Cplex(cpx::Error {
                        code: 0,
                        msg: "Can't open log-file".to_string(),
                    }));
                }
                status = cpx::setlogfile(cpx::env(), logfile);
                if status != 0 {
                    break;
                }

                net = cpx::NETcreateprob(cpx::env(), &mut status, const_cstr!("mcf").as_ptr());

            }

            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(Error::Cplex(cpx::Error {
                    code: status,
                    msg: CString::from_raw(msg).to_string_lossy().into_owned(),
                }));
            }
        }

        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<()> {
        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<()> {
        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<()> {
        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<()> {
        Ok(trycpx!(cpx::NETprimopt(cpx::env(), self.net)))
    }

    pub fn objective(&self) -> Result<Real> {
        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> {
        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<()> {
        let fname = CString::new(filename).unwrap();
        Ok(trycpx!(cpx::NETwriteprob(cpx::env(),
                                     self.net,
                                     fname.as_ptr(),
                                     ptr::null_mut())))
    }
}

//

//! The main bundle method solver.

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

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

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



quick_error! {
    /// A solver error.
    #[derive(Debug)]
    pub enum Error {
        /// An error occurred during evaluation of the oracle.
        Eval(err: Box<error::Error>) {
            cause(&**err)
            description(err.description())
            display("Evaluation error: {}", err)
        }

        /// An error occurred during update of the oracle.
        Update(err: Box<error::Error>) {
            cause(&**err)
            description(err.description())
            display("Update error: {}", err)
        }

        /// Error solving the master problem.
        Master(err: Box<error::Error>) {
            cause(&**err)
            description(err.description())
            display("Master problem error: {}", err)
            from(err: master::Error) -> (Box::new(err))
        }

        /// The oracle did not return a minorant.
        NoMinorant {
            description("No minorant")
            display("The oracle did not return a minorant")
        }

        /// The dimension of some data is wrong.

        Dimension(msg: &'static str) {
            description("Dimension error")

            display("Dimension error: {}", msg)
        }














        /// Some parameter has an invalid value.
        Parameter(msg: String) {
            description("Parameter error")
            display("Parameter error: {}", msg)
        }

        /// The lower bound of a variable is larger than the upper bound.
        InvalidBounds(lower: Real, upper: Real) {
            description("invalid bounds")
            display("Invalid bounds, lower:{} upper:{}", lower, upper)
        }

        /// The value of a variable is outside its bounds.
        ViolatedBounds(lower: Real, upper: Real, value: Real) {
            description("violated bounds")
            display("Violated bounds, lower:{} upper:{} value:{}", lower, upper, value)
        }

        /// The variable index is out of bounds.
        InvalidVariable(index: usize, nvars: usize) {
            description("invalid variable")
            display("Variable index out of bounds, got:{} must be < {}", index, nvars)
        }

        /// Iteration limit has been reached.
        IterationLimit(limit: usize) {
            description("iteration limit reached")
            display("The iteration limit of {} has been reached.", limit)
        }
    }
}


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

/**
 * The current state of the bundle method.
 *
 * 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

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

impl Default for SolverParams {

     * 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()),
            bounds: vec![],
            cur_y: 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.
     *

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

    /// Initialize the solver.
    pub fn init(&mut self) -> Result<()> {
        try!(self.params.check());
        if self.cur_y.len() != self.problem.num_variables() {
            self.cur_valid = false;
            self.cur_y.init0(self.problem.num_variables());
        }

        let lb = self.problem.lower_bounds();
        let ub = self.problem.upper_bounds();
        self.bounds.clear();
        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(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;

        self.start_time = Instant::now();

        Ok(())
    }

    /// Solve the problem.
    pub fn solve(&mut self) -> Result<()> {
        const LIMIT: usize = 10_000;

        if self.solve_iter(LIMIT)? {
            Ok(())
        } else {
            Err(Error::IterationLimit(LIMIT))
        }
    }

    /// Solve the problem but stop after `niter` iterations.
    ///
    /// The function returns `Ok(true)` if the termination criterion
    /// has been satisfied. Otherwise it returns `Ok(false)` or an
    /// error code.
    ///
    /// If this function is called again, the solution process is
    /// continued from the previous point. Because of this one must
    /// call `init()` before the first call to this function.
    pub fn solve_iter(&mut self, niter: usize) -> Result<bool> {
        for _ in 0..niter {
            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 {
                return Ok(true)
            }
        }
        Ok(false)
    }

    /// Called to update the problem.
    ///
    /// Calling this function typically triggers the problem to
    /// separate new constraints depending on the current solution.
    fn update_problem(&mut self, term: Step) -> Result<bool> {
        let updates = {
            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
                    };
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, 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));
                    }
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, value));
                }
                Update::MoveVariable { index, value } => {
                    if index >= self.bounds.len() {
                        return Err(Error::InvalidVariable(index, self.bounds.len()));


                    }
                    let (lower, upper) = self.bounds[index];
                    if value < lower || value > upper {
                        return Err(Error::ViolatedBounds(lower, upper, value));
                    }
                    newvars.push((Some(index), lower - value, upper - value, value));
                }
            }
        }

        if !newvars.is_empty() {

    /**
     * Initializes the master problem.
     *
     * The oracle is evaluated once at the initial center and the
     * master problem is initialized with the returned subgradient
     * information.
     */
    fn init_master(&mut self) -> Result<()> {
        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::<MinimalMaster>::new().unwrap())
        } else {
            debug!("Use CPLEX master problem");
            Box::new(BoxedMasterProblem::<CplexMaster>::new().unwrap())
        };

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

        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))),
            };
            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: try!(self.master.add_minorant(i, minorant)),
                    multiplier: 0.0,
                    primal: Some(primal),
                });
            } else {
                return Err(Error::NoMinorant);
            }
        }

        self.cur_valid = true;

        // Solve the master problem once to compute the initial
        // subgradient.

        debug!("Init master completed");

        Ok(())
    }


    /// Solve the model (i.e. master problem) to compute the next candidate.
    fn solve_model(&mut self) -> Result<()> {
        try!(self.master.solve(self.cur_val));
        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;

        debug!("  nxt_mod ={}", self.nxt_mod);
        debug!("  expected={}", self.expected_progress);
        Ok(())
    }


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

        self.master.set_weight(new_weight);
        self.cnt_null += 1;
        debug!("Null Step");
    }

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

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

        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 result = match self.problem.evaluate(fidx, &self.nxt_y, nullstep_bnd, relprec) {
                Ok(r) => r,
                Err(err) => return Err(Error::Eval(Box::new(err))),
            };

            let fun_ub = result.objective();

            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;