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

//! The main bundle method solver.

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

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

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

    }
}


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

    /**
     * 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 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 {
    fn default() -> SolverParams {
        SolverParams {
            max_bundle_size: 50,

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


/// Information about a minorant.
#[derive(Debug, Clone)]
struct MinorantInfo<Pr> {
    /// The minorant's index in the master problem
    index: usize,
    /// Current multiplier.
    multiplier: Real,
    /// Primal associated with this minorant.
    primal: Option<Pr>,
}

/// State information for the update callback.
pub struct UpdateState<'a, Pr:'a> {
    /// Current model minorants.
    minorants: &'a [Vec<MinorantInfo<Pr>>],
    /// The last step type.
    pub step: Step,
}

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

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

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

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

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

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

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

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

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

    /// 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<MinorantIndex=usize>>,

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


impl<P, Pr, E> Solver<P, Pr, E>
    where P : for <'a> FirstOrderProblem<'a, Primal=Pr,EvalResult=E>,
          E : Evaluation<Pr>
{
    /**
     * Create a new solver for the given problem.
     *
     * Note that the solver owns the problem, so you cannot use the
     * same problem description elsewhere as long as it is assigned to
     * the solver. However, it is possible to get a reference to the
     * internally stored problem using `Solver::problem()`.
     */
    pub fn new_params(problem : P, params : SolverParams) -> Result<Solver<P, Pr, E>> {
        Ok(Solver{
            problem: problem,
            params:  params,
            terminator: Box::new(StandardTerminator { termination_precision: 1e-3 }),
            weighter: Box::new(HKWeighter::new()),
            cur_y : dvec![],
            cur_val : 0.0,
            cur_mod : 0.0,
            cur_vals : dvec![],
            cur_mods : dvec![],
            cur_valid : false,
            nxt_d : dvec![],
            nxt_y : dvec![],
            nxt_val : 0.0,
            nxt_mod : 0.0,
            nxt_vals : dvec![],
            nxt_mods : dvec![],
            new_cutval: 0.0,
            sgnorm: 0.0,
            expected_progress: 0.0,
            cnt_descent: 0,
            cnt_null: 0,
            start_time: Instant::now(),
            master: match BoxedMasterProblem::<MinimalMaster>::new() {
                Ok(master) => Box::new(master),
                Err(err) => return Err(Error::Master(Box::new(err))),
            },
            minorants: vec![],
        })
    }

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

    /**
     * Set the first order problem description associated with this
     * solver.
     *
     * Note that the solver owns the problem, so you cannot use the
     * same problem description elsewhere as long as it is assigned to
     * the solver. However, it is possible to get a reference to the
     * internally stored problem using `Solver::problem()`.
     */
    pub fn set_problem(&mut self, problem : P) {
        self.problem = problem;
    }

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

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

        Ok(())
    }

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

    /// 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 state = UpdateState {minorants: &self.minorants, step: term};
        let updates = 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} => {
                    newvars.push((lower, upper));
                },
            }
        }

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

    /// Return the current aggregated primal information for a subproblem.
    ///
    /// This function returns all currently used minorants $x_i$ along
    /// with their coefficients $\alpha_i$. The aggregated primal can
    /// be computed by combining the minorants $\bar{x} =
    /// \sum_{i=1}\^m \alpha_i x_i$.
    pub fn aggregated_primals(&self, subproblem : usize) -> Vec<(Real, &Pr)> {
        self.minorants[subproblem].iter().map(|m| {
            (m.multiplier, m.primal.as_ref().unwrap())
        }).collect()
    }

    fn show_info(&self, step: Step) {
        let time = self.start_time.elapsed();
        info!("{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4} {:4} {:4}{:1}  {:9.4} {:9.4} {:12.6e}({:12.6e}) {:12.6e}",
              if step == Step::Term { "_endit" } else { "endit " },
              time.as_secs() / 3600,
              (time.as_secs() / 60) % 60,
              time.as_secs() % 60,
              time.subsec_nanos() / 10000000,
              self.cnt_descent,
              self.cnt_descent + self.cnt_null,
              self.master.cnt_updates(),
              if step == Step::Descent { "*" } else { " " },
              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<()> {
        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(|v| v.to_dense());
        let ub = self.problem.upper_bounds().map(|v| v.to_dense());

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

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

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

        self.cur_val = 0.0;
        for i in 0..m {
            let result = match self.problem.evaluate(i, &self.cur_y, INFINITY, 0.0) {
                Ok(r) => r,
                Err(err) => return Err(Error::Eval(Box::new(err))),
            };
            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.
        //
        // We could compute that subgradient directly by
        // adding up the initial minorants, but this would not include
        // the eta terms. However, this is a heuristic anyway because
        // we assume an initial weight of 1.0, which, in general, will
        // *not* be the initial weight for the first iteration.
        self.master.set_weight(1.0);
        try!(self.master.solve(self.cur_val));
        self.sgnorm = self.master.get_dualoptnorm2().sqrt();

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

        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;

        // update multiplier from master solution
        for i in 0..self.problem.num_subproblems() {
            for m in self.minorants[i].iter_mut() {
                m.multiplier = self.master.multiplier(m.index);
            }
        }

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

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

        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;

            // move center of minorant to cur_y
            nxt_minorant.move_center(-1.0, &self.nxt_d);
            self.new_cutval += nxt_minorant.constant;
            self.minorants[fidx].push(MinorantInfo{
                index: try!(self.master.add_minorant(fidx, nxt_minorant)),
                multiplier: 0.0,
                primal: Some(nxt_primal),
            });
        }

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

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