RsBundle  Check-in [4cd78f8116]

[root]
name = "bundle"
version = "0.1.0"
dependencies = [
 "env_logger 0.3.5 (registry+https://github.com/rust-lang/crates.io-index)",
 "libc 0.2.16 (registry+https://github.com/rust-lang/crates.io-index)",
 "log 0.3.6 (registry+https://github.com/rust-lang/crates.io-index)",
 "quick-error 1.1.0 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "aho-corasick"
version = "0.5.3"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "memchr 0.1.11 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "env_logger"
version = "0.3.5"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "log 0.3.6 (registry+https://github.com/rust-lang/crates.io-index)",
 "regex 0.1.77 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "kernel32-sys"
version = "0.2.2"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "winapi 0.2.8 (registry+https://github.com/rust-lang/crates.io-index)",
 "winapi-build 0.1.1 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "libc"
version = "0.2.16"
source = "registry+https://github.com/rust-lang/crates.io-index"

[[package]]
name = "log"
version = "0.3.6"
source = "registry+https://github.com/rust-lang/crates.io-index"

[[package]]
name = "memchr"
version = "0.1.11"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "libc 0.2.16 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "quick-error"
version = "1.1.0"
source = "registry+https://github.com/rust-lang/crates.io-index"

[[package]]
name = "regex"
version = "0.1.77"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "aho-corasick 0.5.3 (registry+https://github.com/rust-lang/crates.io-index)",
 "memchr 0.1.11 (registry+https://github.com/rust-lang/crates.io-index)",
 "regex-syntax 0.3.5 (registry+https://github.com/rust-lang/crates.io-index)",
 "thread_local 0.2.7 (registry+https://github.com/rust-lang/crates.io-index)",
 "utf8-ranges 0.1.3 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "regex-syntax"
version = "0.3.5"
source = "registry+https://github.com/rust-lang/crates.io-index"

[[package]]
name = "thread-id"
version = "2.0.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "kernel32-sys 0.2.2 (registry+https://github.com/rust-lang/crates.io-index)",
 "libc 0.2.16 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "thread_local"
version = "0.2.7"
source = "registry+https://github.com/rust-lang/crates.io-index"
dependencies = [
 "thread-id 2.0.0 (registry+https://github.com/rust-lang/crates.io-index)",
]

[[package]]
name = "utf8-ranges"
version = "0.1.3"
source = "registry+https://github.com/rust-lang/crates.io-index"

[[package]]
name = "winapi"
version = "0.2.8"
source = "registry+https://github.com/rust-lang/crates.io-index"

[[package]]
name = "winapi-build"
version = "0.1.1"
source = "registry+https://github.com/rust-lang/crates.io-index"

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

build = "build.rs"

[dependencies]
libc = "^0.2.6"
quick-error = "*"
log = "*"
env_logger = "*"

/*
 * 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 std::env;

fn main() {
    match env::var("GUROBI_HOME") {
        Ok(grb) => {
            println!("cargo:rustc-link-search={}/lib", grb);
            println!("cargo:rustc-link-lib=gurobi65");
        },
        Err(_) => panic!("GUROBI_HOME environment variable not set"),
    }
}

#!/bin/sh

rustdoc src/lib.rs --html-in-header src/mathjax.txt --no-defaults --passes collapse-docs --passes unindent-comments --crate-name bundle -o $PWD/target/doc -L dependency=$PWD/target/debug -L dependency=$PWD/target/debug/deps --extern libc=$(ls -1 $PWD/target/debug/deps/liblibc-*.rlib) --extern quick_error=$(ls -1 $PWD/target/debug/deps/libquick_error-*.rlib) --extern log=$(ls -1 $PWD/target/debug/deps/liblog-*.rlib)

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

#[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 EvalResult = SimpleEvaluation;

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

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

        for i in 0..2 {
            for j in 0..2 {
                g[i] += self.a[i][j] * x[j];
            }
            objective += x[i] * (g[i] + self.b[i]);
            g[i] = 2.0*g[i] + self.b[i];
        }

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

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

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

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

/*
 * 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, Vector, DVector, Minorant};

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 : IntoIterator<Item=Minorant> {
    /// 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.
 */
pub struct SimpleEvaluation {
    pub objective : Real,
    pub minorants : Vec<Minorant>,
}

impl IntoIterator for SimpleEvaluation {
    type Item = Minorant;
    type IntoIter = IntoIter<Minorant>;

    fn into_iter(self) -> Self::IntoIter {
        self.minorants.into_iter()
    }
}

impl Evaluation for SimpleEvaluation {
    fn objective(&self) -> Real {
        self.objective
    }
}


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

    /// Custom evaluation result value.
    type EvalResult : Evaluation;

    /// 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$.
     *
     * If the evaluation process reaches a lower bound on the function
     * value at $y$ and this bound is larger than $nullstep_bound$,
     * the evaluation may stop and return the lower bound and a
     * minorant. In this case the function value is guaranteed to be
     * large enough so that the new point is rejected as candidate.
     *
     * The returned objective value should be an upper bound on the
     * 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>;
}

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

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

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

#[macro_use]
extern crate log;

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

#[macro_export]
macro_rules! dvec {
    ( $ elem : expr ; $ n : expr ) => { DVector(vec![$elem; $n]) };
    ( $ ( $ x : expr ) , * ) => { DVector(vec![$($x),*]) };
    ( $ ( $ x : expr , ) * ) => { DVector(vec![$($x,)*]) };
}

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

pub mod minorant;
pub use minorant::Minorant;

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

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

mod hkweighter;
pub use hkweighter::HKWeighter;

mod master;

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

    /**
     * Update box multipliers $\eta$.
     *
     * This function solves the dual problem for fixed aggregated
     * minorant w.r.t. $\eta$. When called, the variable `self.primopt`
     * should contain the primal solution (i.e. the new candidate point
     * d) without the influence of $\eta$.
     */
    fn update_box_multipliers(&mut self) -> bool {
        let mut updated_eta = false;
        let weight = self.master.weight();

        if self.eta.len() != self.lb.len() {
            self.eta.resize(self.lb.len(), 0.0);
        }
        for i in 0..self.lb.len() {
            let mut b = self.lb[i];
            let x = self.primopt[i];
            if x >= b {
                b = self.ub[i];
                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.
    fn eta_cutval(&self) -> Real {
        let mut val = 0.0;
        for i in 0..self.lb.len() {
            if self.eta[i] >= 0.0 {
                if self.lb[i] > NEG_INFINITY {
                    val += self.lb[i] * self.eta[i];
                }
            } else {
                if self.ub[i] < INFINITY {
                    val += self.ub[i] * self.eta[i];
                }
            }
        }
        return val;
    }


    /**
     * Return $\\|G \alpha - \eta\\|_2\^2$.
     *
     * This is the norm-square of the dual optimal solution including
     * the current box-multipliers $\eta$.
     */
    fn get_norm_subg2(&self) -> Real {
        let dualopt = self.master.dualopt();
        let mut norm2 = 0.0;
        for i in 0..self.eta.len() {
            let x = dualopt[i] - self.eta[i];
            norm2 += x * x;
        }
        return norm2;
    }
}


impl<M : UnconstrainedMasterProblem> MasterProblem for BoxedMasterProblem<M> {
    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) -> usize { self.master.num_minorants() }

    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<()> {
        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 solve(&mut self, center_value: Real) -> Result<()> {
        debug!("Solve Master");
        debug!("  lb      ={}", self.lb);
        debug!("  ub      ={}", self.ub);

        // for (fidx, mins) in self.master.minorants() {
        //     for m in mins {
        //         debug!("  {}:min[{},{}] = {}",  m)
        //     }
        // }

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

            // compute value of the linearized model
            self.dualoptnorm2 = self.get_norm_subg2();
            let linval = self.master.dualopt().dot(&self.primopt) + self.master.dualopt_cutval();
            assert!((self.eta.dot(&self.primopt) - self.eta_cutval()).abs() < 1e-6); // verify complementarity
            let augval = linval + 0.5 * self.dualoptnorm2 / self.master.weight();

            let mut cutval = linval;
            if updated_eta {
                cutval = self.master.eval_model(&self.primopt)
            }

            let curval = center_value;
            let model_prec = (cutval - linval) / (curval - linval).max(1e-16);

            debug!("Eta Test");
            debug!("  dualnorm2={}", self.dualoptnorm2);
            debug!("  linval={}", linval);
            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]) {
        self.master.aggregate(fidx, mins)
    }

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

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

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

    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() {
            self.lb[i] -= alpha * d[i];
            self.ub[i] -= alpha * d[i];
        }
    }

    fn set_max_updates(&mut self, max_updates: usize) {
        BoxedMasterProblem::set_max_updates(self, max_updates);
    }

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

/*
 * 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 std::result;
use std::os::raw::c_char;
use std::ptr;
use std::ffi::{CString, CStr};

use libc::{c_int, c_double};

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

enum GRBenv {}
enum GRBmodel {}

quick_error! {
    /// A solver error.
    #[derive(Debug)]
    pub enum Error {
        Gurobi(code : c_int, msg: *const c_char) {
            description("Gurobi Error")
            display("Gurobi Error ({}): {}", code, unsafe{CStr::from_ptr(*msg).to_string_lossy().into_owned()})
        }
    }
}

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

pub struct GurobiMaster {
    env: *mut GRBenv,
    model: *mut GRBmodel,

    dualopt: DVector,
}


extern "C" {
    fn GRBloadenv(env: *mut *mut GRBenv, logfile: *const c_char) -> c_int;
    fn GRBfreeenv(env: *mut GRBenv);
    fn GRBgeterrormsg(env: *mut GRBenv) -> *const c_char;

    fn GRBnewmodel(env: *mut GRBenv,
                   model: *mut *mut GRBmodel,
                   pname: *const c_char,
                   numvars: c_int,
                   obj: *const c_double,
                   lb: *const c_double,
                   ub: *const c_double,
                   vtype: *const c_char,
                   varnames: *const *const c_char) -> c_int;
    fn GRBfreemodel(model: *mut GRBmodel) -> c_int;
}

impl Drop for GurobiMaster {
    fn drop(&mut self) {
        unsafe {
            GRBfreemodel(self.model);
            GRBfreeenv(self.env);
        }
    }
}

impl UnconstrainedMasterProblem for GurobiMaster {
    type Error = Error;

    fn new() -> Result<GurobiMaster> {
        let mut env : *mut GRBenv = ptr::null_mut();

        unsafe {
            let ret = GRBloadenv(&mut env, ptr::null());
            if ret != 0 {
                return Err(Error::Gurobi(ret, GRBgeterrormsg(env)))
            }
        }

        let mut model : *mut GRBmodel = ptr::null_mut();
        unsafe {
            let ret = GRBnewmodel(env,
                                  &mut model, CString::new("master").unwrap().as_ptr(),
                                  0,
                                  ptr::null(),
                                  ptr::null(),
                                  ptr::null(),
                                  ptr::null(),
                                  ptr::null());
            if ret != 0 {
                GRBfreeenv(env);
                return Err(Error::Gurobi(ret, GRBgeterrormsg(env)))
            }
        }

        Ok(GurobiMaster {
            env: env,
            model: model,

            dualopt: dvec![],
        })
    }

    fn add_minorant(&mut self, minorant: Minorant) -> Result<()>{
        Ok(())
    }

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

    fn set_weight(&mut self, weight: Real) {

    }

    fn solve(&mut self, fbound: Real, augbound: Real, relprec: Real) -> Result<()>{
        Ok(())
    }

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

    fn dualopt_cutval(&self) -> Real { 0.0 }

    fn eval_model(&self, y: &DVector) -> Real { 0.0 }
}

/*
 * 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! {
    /// 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 {
    /// The 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.
    fn num_minorants(&self) -> usize;

    /// Add a new minorant to the model.
    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<()>;

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

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

    /// Aggregate the given minorants according to the current solution.
    fn aggregate(&mut self, fidx: usize, mins: &[usize]);

    /// Return the primal optimal solution.
    fn get_primopt(&self) -> DVector;

    /// Return the primal optimal solution value.
    fn get_primoptval(&self) -> Real;

    /// Return $\\|d\^*\\|\^2$ of the current dual optimal solution $d\^*$.
    fn get_dualoptnorm2(&self) -> Real;

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

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

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


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(n: usize) {
            description("Too many minorants")
            display("Solver Error: maximal number of minorants too large (got: {} must be <= 2)", n)
        }

        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
 * it is only a very loose approximation of the objective function.
 *
 * Because of its properties, it can only be used if the problem to be
 * solved has a maximal number of minorants of two and only one
 * subproblem.
 */
pub struct MinimalMaster {
    weight: Real,

    /// The minorants in the model.
    minorants: Vec<Minorant>,
    /// Optimal multipliers.
    opt_mult: DVector,
    /// Optimal aggregated minorant.
    opt_minorant: Minorant,
}


impl UnconstrainedMasterProblem for MinimalMaster {
    type Error = Error;

    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 max_num_minorants(&self) -> usize { 2 }

    fn set_max_num_minorants(&mut self, n: usize) -> Result<()> {
        if n > 2 {
            Err(Error::MaxMinorants(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) -> usize {
        self.minorants.len()
    }

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

    #[allow(unused_variables)]
    fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<()> {
        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);
            let yeta = self.minorants[1].linear.dot(eta);
            let a = yy - 2.0*xy + xx;
            let b = xy - xx - yeta - xeta;

            let mut alpha2 = 0.0;
            if a > 0.0 {
                alpha2 = ((self.minorants[1].constant - self.minorants[0].constant) * self.weight - b) / a;
                alpha2 = alpha2.max(0.0).min(1.0);
            }
            self.opt_mult.resize(2, 0.0);
            self.opt_mult[0] = 1.0 - alpha2;
            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 eval_model(&self, y: &DVector) -> Real {
        let mut result = NEG_INFINITY;
        for m in &self.minorants {
            result = result.max(m.eval(y));
        }
        return result;
    }

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) {
        assert!(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);
            self.opt_mult.truncate(1);
            self.opt_mult[0] = 1.0;
            debug!("  {}", self.minorants[0]);
        }
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        for m in self.minorants.iter_mut() {
            m.move_center(alpha, d);
        }
    }
}

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

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

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

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

pub mod minimal;
pub use self::minimal::MinimalMaster;

// pub mod grb;
// pub use master::grb::GurobiMaster;

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

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

    /// 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 maximal number of minorants.
    fn max_num_minorants(&self) -> usize;

    /// Set the maximal number of minorants.
    fn set_max_num_minorants(&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.
    fn num_minorants(&self) -> usize;

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

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

    /// 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.
    fn aggregate(&mut self, fidx: usize, mins: &[usize]);

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

<script type="text/x-mathjax-config">
  MathJax.Hub.Config({
    tex2jax: {
      inlineMath: [['$','$'], ['\\(','\\)']],
      displayMath: [['$$','$$'], ['\\[','\\]']],
    },
    TeX: {
      extensions: ["AMSmath.js", "AMSsymbols.js"]
    }
  });
</script>
<script type="text/javascript" async
  src="https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS_HTML">
</script>
</head>

/*
 * 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 std::fmt;

/**
 * A linear minorant of a convex function.
 *
 * A linear minorant of a convex function $f \colon \mathbb{R}\^n \to
 * \mathbb{R}$ is a linear function of the form
 *
 *   \\[ l \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto \langle g, x
 *   \rangle + c \\]
 *
 * such that $l(x) \le f(x)$ for all $x \in \mathbb{R}\^n$.
 */
#[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(())
    }
}


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

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

    /**
     * Move the center of the minorant.
     */
    pub fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.constant += alpha * self.linear.dot(d);
    }
}

/*
 * 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 {FirstOrderProblem, Evaluation, HKWeighter};

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

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

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


/// 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
 * 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,
            nxt_mod : $slf.nxt_mod,
            nxt_val : $slf.nxt_val,
            new_cutval : $slf.new_cutval,
            sgnorm : $slf.sgnorm,
            weight: $slf.master.weight(),
            step: $step,
            expected_progress: $slf.expected_progress,
        }
    };
}

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

    /**
     * Maximal age of a minorant before it gets compressed.
     *
     * Minorants that have not been active for that many iterations will be compressed,
     * independent of the size of the bundle.
     */
    pub max_age: 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.
     *
     * Must be in (0,1).
     */
    pub acceptance_factor: Real,

    /**
     * Factor for doing a null step.
     *
     * Factor that guarantees a null step. This factor is used to
     * 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 Default for SolverParams {
    fn default() -> SolverParams {
        SolverParams {
            max_bundle_size: 50,
            max_age: 10,

            nullstep_factor: 0.1,
            acceptance_factor: 0.1,

            min_weight: 0.01,
            max_weight: 1000.0,

            max_updates: 50,
        }
    }
}


/// The step type that has been performed.
#[derive(Clone, Copy, PartialEq, Eq, Debug)]
pub enum Step {
    /// A null step has been performed.
    Null,
    /// A descent step has been performed.
    Descent,
    /// No step but the algorithm has been terminated.
    Term,
}

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

    /// The solver parameter.
    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,

    /// Vector of model values in current point.
    cur_mods: DVector,

    /**
     * Whether the data of the current center is valid.
     *
     * This variable is set to false of the problem data changes so
     * the function is re-evaluated at the center.
     */
    cur_valid: bool,

    /// Direction from current center to candidate.
    nxt_d: DVector,

    /// Current candidate point.
    nxt_y: DVector,

    /// (Upper bound on) function value in candidate.
    nxt_val: Real,

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

    /// DVector of subproblem function values in candidate.
    nxt_vals: DVector,

    /// Vector of model values in candidate point.
    nxt_mods: DVector,

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

    /// Norm of current aggregated subgradient.
    sgnorm: Real,

    /// Expected progress.
    expected_progress: Real,

    /// Number of descent steps.
    cnt_descent: usize,

    /// Number of null steps.
    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>,
}


impl<P> Solver<P>
    where P : for <'a> FirstOrderProblem<'a>
{
    /**
     * 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>> {
        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))),
            }
        })
    }

    /// A new solver with default parameter.
    pub fn new(problem : P) -> Result<Solver<P>> {
        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
    }

    /// Initialize the solver.
    fn init(&mut self) {
        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().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 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;
            }
        }

        let m = self.problem.num_subproblems();
        self.cur_vals.init0(m);
        self.cur_mods.init0(m);
        self.nxt_vals.init0(m);
        self.nxt_mods.init0(m);

        self.start_time = Instant::now();
    }

    /// Solve the problem.
    pub fn solve(&mut self) -> Result<()> {
        self.init();
        for _ in 0..100000 {
            let term = try!(self.step());
            self.show_info(term);
            if term == Step::Term {
                break;
            }
        }
        Ok(())
    }

    fn show_info(&self, step: Step) {
        let time = self.start_time.elapsed();
        info!("{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4}{:1} {:4} {:4}  {:9.4} {:9.4} {:12.6}({:12.6}) {:12.6}",
              if step == Step::Term { "_endit" } else { "endit" },
              time.as_secs() / 3600,
              (time.as_secs() / 60) % 60,
              time.as_secs() % 60,
              time.subsec_nanos() / 1000000,
              self.cnt_descent,
              if step == Step::Descent { "*" } else { " " },
              self.cnt_descent + self.cnt_null,
              self.master.cnt_updates(),
              self.master.weight(),
              self.expected_progress,
              self.nxt_mod,
              self.nxt_val,
              self.cur_val);
    }

    /**
     * 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<()> {
        self.master = Box::new(BoxedMasterProblem::<MinimalMaster>::new().unwrap());

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

        self.master.set_vars(self.problem.num_variables(), lb, ub);
        self.master.set_max_updates(self.params.max_updates);

        let m = self.problem.num_subproblems();

        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) = minorants.next() {
                self.cur_mods[i] = minorant.constant;
                self.cur_mod += self.cur_mods[i];
                try!(self.master.add_minorant(i, minorant));
            } else {
                return Err(Error::NoMinorant);
            }
        }

        self.cur_valid = true;
        self.master.set_weight(1.0);

        let state = current_state!(self, Step::Term);
        let new_weight = self.weighter.weight(&state, &self.params);
        self.master.set_weight(new_weight);

        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!("Model result");
        debug!("  cur_val ={}", self.cur_val);
        debug!("  nxt_mod ={}", self.nxt_mod);
        debug!("  expected={}", self.expected_progress);
        Ok(())
    }


    /// Reduce size of bundle.
    fn compress_bundle(&mut self) {
        if self.master.num_minorants() > 1 {
            let n = self.master.num_minorants();
            self.master.aggregate(0, &(0..n).collect::<Vec<_>>());
        }
    }

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

    /// Perform a null step.
    fn null_step(&mut self) {
        let new_weight = self.weighter.weight(&current_state!(self, Step::Null), &self.params);
        self.master.set_weight(new_weight);
        self.cnt_null += 1;
        debug!("Null Step");
    }

    /// Perform one bundle iteration.
    pub fn step(&mut self) -> Result<Step> {
        if !self.cur_valid {
            // current point needs new evaluation
            try!(self.init_master());
        }

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

        self.compress_bundle();

        let mut nxt_lb = 0.0;
        let mut nxt_ub = 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 = match minorants.next() {
                Some(m) => m,
                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;
            try!(self.master.add_minorant(fidx, nxt_minorant));
        }

        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.nxt_val = nxt_ub;

        // check for potential problems with relative precision of all kinds
        if nxt_lb <= descent_bnd {
            // lower bound gives descent step
            if nxt_ub > descent_bnd {
                // upper bound will produce null-step
                if self.cur_val - nxt_lb > (self.cur_val - self.nxt_mod) * self.params.nullstep_factor.max(0.5) {
                    warn!("Relative precision of returned objective interval enforces null-step.")
                }
            }
        } else {
            // lower bound gives already a null step
            if self.cur_val - nxt_lb > 0.8 * (self.cur_val - self.nxt_mod) {
                // subgradient won't yield much improvement
                warn!("Shallow cut (subgradient won't yield much improvement)");
            }
        }

        debug!("Step");
        debug!("  cur_val    ={}", self.cur_val);
        debug!("  nxt_mod    ={}", self.nxt_mod);
        debug!("  nxt_ub     ={}", self.nxt_val);
        debug!("  descent_bnd={}", descent_bnd);

        // do a descent step or null step
        if nxt_ub <= descent_bnd {
            self.descent_step();
            return Ok(Step::Descent);
        } else {
            self.null_step();
            return Ok(Step::Null);
        }
    }

    /**
     * Return the bound on the function value that enforces a
     * nullstep.
     *
     * If the oracle guarantees that $f(\bar{y}) \ge$ this bound, the
     * bundle method will perform a nullstep.
     *
     * This value is $f(\hat{y}) + \varrho' \cdot \Delta$ where
     * $\Delta = f(\hat{y}) - \hat{f}(\bar{y})$ is the expected
     * progress and $\varrho'$ is the `nullstep_factor`.
     */
    fn get_nullstep_bound(&self) -> Real {
        self.cur_val - self.params.nullstep_factor * (self.cur_val - self.nxt_mod)
    }


    /**
     * Return the bound the function value must be below of to enforce a descent step.
     *
     * If the oracle guarantees that $f(\bar{y}) \le$ this bound, the
     * bundle method will perform a descent step.
     *
     * This value is $f(\hat{y}) + \varrho \cdot \Delta$ where
     * $\Delta = f(\hat{y}) - \hat{f}(\bar{y})$ is the expected
     * progress and $\varrho$ is the `acceptance_factor`.
     */
    fn get_descent_bound(&self) -> Real {
        self.cur_val - self.params.acceptance_factor * (self.cur_val - self.nxt_mod)
    }

    /**
     * Return the required relative precision for the computation.
     */
    fn get_relative_precision(&self) -> Real {
        (0.1 * (self.cur_val - self.get_nullstep_bound()) / (self.cur_val.abs() + 1.0)).min(1e-3)
    }
}

/*
 * 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;
use std::fmt;
use std::ops::{Deref, DerefMut};
use std::cmp::min;

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

/// 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, ")")
            }
        }
    }
}


impl DVector {
    /// Set all elements to 0.
    pub fn init0(&mut self, size: usize) {
        let n = self.len();
        self.resize(size, 0.0);
        for i in 0..min(n, size) {
            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 the inner product with another vector.
    pub fn dot(&self, other: &DVector) -> Real {
        let mut ip = 0.0;
        for i in 0..self.len() {
            ip += self[i] * other[i];
        }
        return ip;
    }

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


impl Vector {
    /**
     * 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)
            }
        }
    }
}