RsBundle  Check-in [b194454b53]

 */

#![allow(non_upper_case_globals)]

//! Example implementation for a capacitated facility location problem.

use better_panic;
use crossbeam::channel::unbounded as channel;
use log::{info, Level};
use rustop::opts;
use std::error::Error;
use std::fmt::Write;
use std::io::Write as _;
use std::sync::Arc;

use env_logger::{self, fmt::Color};
use ordered_float::NotNan;
use threadpool::ThreadPool;

use bundle::parallel::{
    DefaultSolver as ParallelSolver, EvalResult, FirstOrderProblem as ParallelProblem,
    NoBundleSolver as NoParallelSolver, ResultSender,
};
use bundle::{dvec, DVector, Minorant, Real};
use bundle::{DefaultSolver, FirstOrderProblem, NoBundleSolver, SimpleEvaluation};

const Nfac: usize = 3;
const Ncus: usize = 5;
const F: [Real; Nfac] = [1000.0, 1000.0, 1000.0];
const CAP: [Real; Nfac] = [500.0, 500.0, 500.0];
const C: [[Real; Ncus]; Nfac] = [
    [4.0, 5.0, 6.0, 8.0, 10.0], //
    [6.0, 4.0, 3.0, 5.0, 8.0],  //
    [9.0, 7.0, 4.0, 3.0, 4.0],  //
];
const DEMAND: [Real; 5] = [80.0, 270.0, 250.0, 160.0, 180.0];

#[derive(Debug)]
enum EvalError {
    Customer(usize),
    NoObjective(usize),
}

impl std::fmt::Display for EvalError {
    fn fmt(&self, fmt: &mut std::fmt::Formatter) -> Result<(), std::fmt::Error> {
        match self {
            EvalError::Customer(i) => writeln!(fmt, "Customer subproblem {} failed", i),
            EvalError::NoObjective(i) => writeln!(fmt, "No objective value generated for subproblem {}", i),
        }
    }
}

impl Error for EvalError {}

struct CFLProblem {

            constant: objective,
            linear: subg,
        },
        primal,
    ))
}

impl FirstOrderProblem for CFLProblem {
    type Err = EvalError;

    type Primal = DVector;

    type EvalResult = SimpleEvaluation<Self::Primal>;

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

    fn lower_bounds(&self) -> Option<Vec<Real>> {
        Some(vec![0.0; FirstOrderProblem::num_variables(self)])
    }

    fn num_subproblems(&self) -> usize {
        Nfac + Ncus
    }

    fn evaluate(
        &mut self,
        index: usize,
        lambda: &[Real],
        _nullstep_bound: Real,
        _relprec: Real,
    ) -> Result<Self::EvalResult, Self::Err> {
        let (tx, rx) = channel();
        ParallelProblem::evaluate(self, index, Arc::new(lambda.iter().cloned().collect()), index, tx)?;
        let mut objective = None;
        let mut minorants = vec![];

        for r in rx {
            match r {
                Ok(EvalResult::ObjectiveValue { value, .. }) => objective = Some(value),
                Ok(EvalResult::Minorant { minorant, primal, .. }) => minorants.push((minorant, primal)),
                _ => break,
            }
        }

        Ok(SimpleEvaluation {
            objective: objective.ok_or(EvalError::NoObjective(index))?,
            minorants,
        })
    }
}

impl ParallelProblem for CFLProblem {
    type Err = EvalError;

    type Primal = DVector;

    fn num_variables(&self) -> usize {
        Nfac

    let (args, _) = opts! {
        synopsis "Solver a simple capacitated facility location problem";
        opt minimal:bool, desc:"Use the minimal master model";
    }
    .parse_or_exit();
    if !args.minimal {
        let mut slv = DefaultSolver::new(CFLProblem::new())?;
        slv.params.max_bundle_size = 5;
        slv.terminator.termination_precision = 1e-9;
        slv.solve()?;
        show_primals(|i| slv.aggregated_primals(i))?;

        let mut slv = ParallelSolver::<_>::new(CFLProblem::new());
        slv.terminator.termination_precision = 1e-9;
        slv.master.max_bundle_size = 5;
        slv.solve()?;

        show_primals(|i| slv.aggregated_primal(i).unwrap())?;
    } else {
        let mut slv = NoBundleSolver::new(CFLProblem::new())?;
        slv.params.max_bundle_size = 2;
        slv.terminator.termination_precision = 1e-5;
        slv.solve()?;
        show_primals(|i| slv.aggregated_primals(i))?;

        let mut slv = NoParallelSolver::<_>::new(CFLProblem::new());
        slv.terminator.termination_precision = 1e-5;
        slv.solve()?;

        show_primals(|i| slv.aggregated_primal(i).unwrap())?;
    }

    Ok(())
}

use rustop::opts;
use std::io::Write;

use bundle::master::{Builder as MasterBuilder, MasterProblem};
use bundle::mcf::MMCFProblem;
use bundle::parallel;
use bundle::{terminator::StandardTerminator, weighter::HKWeighter};
use bundle::{DefaultSolver, FirstOrderProblem, NoBundleSolver, Solver};
use bundle::{FullMasterBuilder, MinimalMasterBuilder};

use std::error::Error;
use std::result::Result;

fn solve_standard<M>(mut slv: Solver<MMCFProblem, StandardTerminator, HKWeighter, M>) -> Result<(), Box<dyn Error>>
where
    M: MasterBuilder + Default,
    M::MasterProblem: MasterProblem<MinorantIndex = usize>,
{
    slv.weighter.set_weight_bounds(1e-1, 100.0);
    slv.terminator.termination_precision = 1e-6;
    slv.solve()?;

    let costs: f64 = (0..slv.problem().num_subproblems())
        .map(|i| {
            let aggr_primals = slv.aggregated_primals(i);
            slv.problem().get_primal_costs(i, &aggr_primals)
        })
        .sum();
    info!("Primal costs: {}", costs);
    Ok(())
}

fn solve_parallel<M>(master: M, mmcf: MMCFProblem) -> Result<(), Box<dyn Error>>
where
    M: MasterBuilder,
    M::MasterProblem: MasterProblem<MinorantIndex = usize>,
{
    let mut slv = parallel::Solver::<_, StandardTerminator, HKWeighter, M>::with_master(mmcf, master);
    slv.weighter.set_weight_bounds(1e-1, 100.0);

    }
    .parse_or_exit();

    let filename = args.file;
    info!("Reading instance: {}", filename);

    if !args.minimal {
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            if args.aggregated {
                mmcf.multimodel = false;
                mmcf.set_separate_constraints(false);
            } else {
                mmcf.multimodel = true;
                mmcf.set_separate_constraints(args.separate);
            }

            let mut solver = DefaultSolver::new(mmcf)?;
            solver.params.max_bundle_size = if args.bundle_size <= 1 {
                if args.aggregated {
                    50
                } else {
                    5
                }
            } else {
                args.bundle_size
            };
            solve_standard(solver)?;
        }

        println!("---------------------------------");
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            mmcf.set_separate_constraints(args.separate);
            mmcf.multimodel = true;

            let mut master = FullMasterBuilder::default();
            if args.aggregated {
                master.max_bundle_size(if args.bundle_size <= 1 { 50 } else { args.bundle_size });
                master.use_full_aggregation();
            } else {
                master.max_bundle_size(if args.bundle_size <= 1 { 5 } else { args.bundle_size });
            }
            solve_parallel(master, mmcf)?;
        }
    } else {
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            mmcf.multimodel = false;

            let mut solver = NoBundleSolver::new(mmcf)?;
            solver.params.max_bundle_size = 2;
            solve_standard(solver)?;
        }

        println!("---------------------------------");
        {
            let mut mmcf = MMCFProblem::read_mnetgen(&filename)?;
            mmcf.set_separate_constraints(args.separate);
            mmcf.multimodel = true;

            let master = MinimalMasterBuilder::default();
            solve_parallel(master, mmcf)?;
        }
    }

    Ok(())
}

use std::sync::Arc;
use std::thread;

use bundle::parallel::{
    DefaultSolver as ParallelSolver, EvalResult, FirstOrderProblem as ParallelProblem,
    NoBundleSolver as NoParallelSolver, ResultSender,
};
use bundle::{DVector, DefaultSolver, FirstOrderProblem, Minorant, NoBundleSolver, Real, SimpleEvaluation};

#[derive(Clone)]
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 FirstOrderProblem for QuadraticProblem {
    type Err = Box<dyn Error + Send + Sync>;
    type Primal = ();
    type EvalResult = SimpleEvaluation<()>;

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

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

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

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

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

impl ParallelProblem for QuadraticProblem {
    type Err = Box<dyn Error + Send + Sync>;
    type Primal = ();

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

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

    fn start(&mut self) {}

    fn stop(&mut self) {}

    fn evaluate<I>(
        &mut self,
        i: usize,
        y: Arc<DVector>,
        index: I,
        tx: ResultSender<I, Self::Primal, Self::Err>,
    ) -> Result<(), Self::Err>
    where
        I: Send + Copy + 'static,
    {
        let y = y.clone();
        let mut p = self.clone();
        thread::spawn(move || match FirstOrderProblem::evaluate(&mut p, i, &y, 0.0, 0.0) {



            Ok(res) => {










                tx.send(Ok(EvalResult::ObjectiveValue {
                    index,
                    value: res.objective,
                }))
                .unwrap();
                for (minorant, primal) in res.minorants {
                    tx.send(Ok(EvalResult::Minorant {
                        index,
                        minorant,



                        primal,
                    }))
                    .unwrap();
                }
            }
            Err(err) => tx.send(Err(err)).unwrap(),
        });
        Ok(())
    }
}

fn main() -> Result<(), Box<dyn Error>> {
    better_panic::install();

    let (args, _) = opts! {
        synopsis "Solver a simple quadratic optimization problem";
        opt minimal:bool, desc:"Use the minimal master model";
    }
    .parse_or_exit();

    {
        let f = QuadraticProblem::new();
        if !args.minimal {
            let mut solver = DefaultSolver::new(f).map_err(|e| format!("{}", e))?;
            solver.weighter.set_weight_bounds(1.0, 1.0);
            solver.solve().map_err(|e| format!("{}", e))?;
        } else {
            let mut solver = NoBundleSolver::new(f).map_err(|e| format!("{}", e))?;
            solver.params.max_bundle_size = 2;
            solver.weighter.set_weight_bounds(1.0, 1.0);
            solver.solve().map_err(|e| format!("{}", e))?;
        }
    }

    println!("-------------------------");

    {
        let f = QuadraticProblem::new();
        if !args.minimal {
            let mut solver = ParallelSolver::<_>::new(f);
            solver.solve().map_err(|e| format!("{}", e))?;
        } else {
            let mut solver = NoParallelSolver::<_>::new(f);
            solver.solve().map_err(|e| format!("{}", e))?;
        }
    }

    Ok(())
}

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

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

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

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

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

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

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

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

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

/// Problem update information.
///
/// The solver calls the `update` method of the problem regularly.
/// This method can modify the problem by adding (or removing)
/// variables. The possible updates are encoded in this type.
#[derive(Debug, Clone, Copy)]
pub enum Update {
    /// Add a variable with bounds.
    ///
    /// The initial value of the variable will be the feasible value
    /// closest to 0.
    AddVariable { lower: Real, upper: Real },
    /// Add a variable with bounds and initial value.
    AddVariableValue { lower: Real, upper: Real, value: Real },
    /// Change the current value of a variable. The bounds remain
    /// unchanged.
    MoveVariable { index: usize, value: Real },
}

/**
 * Trait for implementing a first-order problem description.
 *
 */
pub trait FirstOrderProblem {
    /// Error raised by this oracle.
    type Err;

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

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

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

    /**
     * Return the lower bounds on the variables.
     *
     * If no lower bounds a specified, $-\infty$ is assumed.
     *
     * The lower bounds must be less then or equal the upper bounds.
     */
    fn lower_bounds(&self) -> Option<Vec<Real>> {
        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<Vec<Real>> {
        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(
        &mut self,
        i: usize,
        y: &[Real],
        nullstep_bound: Real,
        relprec: Real,
    ) -> Result<Self::EvalResult, Self::Err>;

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

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

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

pub mod minorant;
pub use crate::minorant::{Aggregatable, Minorant};

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

pub mod solver;
pub use crate::solver::{BundleState, FullMasterBuilder, IterationInfo, Solver, SolverParams, Step, UpdateState};

pub mod parallel;

pub mod weighter;

pub mod terminator;

pub mod master;

pub mod mcf;

/// The minimal bundle builder.
pub type MinimalMasterBuilder = master::boxed::Builder<master::minimal::Builder>;

/// The default bundle solver with general master problem.
pub type DefaultSolver<P> = Solver<P, terminator::StandardTerminator, weighter::HKWeighter, FullMasterBuilder>;

/// A bundle solver with a minimal cutting plane model.
pub type NoBundleSolver<P> = Solver<P, terminator::StandardTerminator, weighter::HKWeighter, MinimalMasterBuilder>;

//

use crate::mcf;
use crate::parallel::{
    EvalResult, FirstOrderProblem as ParallelProblem, ResultSender, Update as ParallelUpdate, UpdateSender,
    UpdateState as ParallelUpdateState,
};
use crate::{Aggregatable, DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation, Update, UpdateState};

use itertools::izip;
use log::{debug, warn};
use num_traits::Float;
use threadpool::ThreadPool;

use std::f64::INFINITY;

            .map(|elem| elem.val * primal[elem.ind])
            .sum::<Real>();
        rhs - lhs
    }
}

impl MMCFProblem {




    pub fn read_mnetgen(basename: &str) -> std::result::Result<MMCFProblem, MMCFReadError> {
        let mut buffer = String::new();
        {
            let mut f = File::open(&format!("{}.nod", basename))?;
            f.read_to_string(&mut buffer)?;
        }
        let fnod = buffer

            }
        }

        aggr
    }
}

impl FirstOrderProblem for MMCFProblem {
    type Err = Error;

    type Primal = Vec<DVector>;

    type EvalResult = SimpleEvaluation<Vec<DVector>>;

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

    fn lower_bounds(&self) -> Option<Vec<Real>> {
        Some(vec![0.0; self.active_constraints.len()])
    }

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

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

    fn evaluate(&mut self, fidx: usize, y: &[Real], _nullstep_bound: Real, _relprec: Real) -> Result<Self::EvalResult> {
        let (objective, subg, sol) = if self.multimodel {
            let (objective, subg, sol) = self.subs[fidx]
                .write()
                .unwrap()
                .evaluate(y, self.active_constraints.iter().cloned())?;
            (objective, subg, vec![sol])
        } else {
            let mut objective = 0.0;
            let mut subg = dvec![0.0; y.len()];
            let mut sols = Vec::with_capacity(self.subs.len());
            for sub in &mut self.subs {
                let (obj, sg, sol) = sub
                    .write()
                    .unwrap()
                    .evaluate(y, self.active_constraints.iter().cloned())?;
                objective += obj;
                subg.add_scaled(1.0, &sg);
                sols.push(sol);
            }
            (objective, subg, sols)
        };
        Ok(SimpleEvaluation {
            objective,
            minorants: vec![(
                Minorant {
                    constant: objective,
                    linear: subg,
                },
                sol,
            )],
        })
    }

    fn update(&mut self, state: &UpdateState<Self::Primal>) -> Result<Vec<Update>> {
        if self.inactive_constraints.is_empty() {
            return Ok(vec![]);
        }

        let nold = self.active_constraints.len();
        let subs = &self.subs;

        // if state.step != Step::Descent && !self.active_constraints.is_empty() {
        //     return Ok(vec![]);
        // }

        let newconstraints = self
            .inactive_constraints
            .iter()
            .map(|&cidx| {
                subs.iter()
                    .enumerate()
                    .map(|(fidx, sub)| {
                        let primals = state.aggregated_primals(fidx);
                        let aggr = Aggregatable::combine(primals.into_iter().map(|(alpha, x)| (alpha, &x[0])));
                        sub.read().unwrap().evaluate_constraint(&aggr, cidx)
                    })
                    .sum::<Real>()
            })
            .enumerate()
            .filter_map(|(i, sg)| if sg < 1e-3 { Some(i) } else { None })
            .collect::<Vec<_>>();

        let inactive = &mut self.inactive_constraints;
        self.active_constraints
            .extend(newconstraints.into_iter().rev().map(|i| inactive.swap_remove(i)));

        // *** The following code needs `drain_filter`, which is experimental as
        // of rust 1.36 ***

        // self.active_constraints
        //     .extend(self.inactive_constraints.drain_filter(|&mut cidx| {
        //         subs.iter()
        //             .enumerate()
        //             .map(|(fidx, sub)| {
        //                 let primals = state.aggregated_primals(fidx);
        //                 let aggr = Aggregatable::combine(primals.into_iter().map(|(alpha, x)| (alpha, &x[0])));
        //                 sub.read().unwrap().evaluate_constraint(&aggr, cidx)
        //             })
        //             .sum::<Real>() < -1e-3
        //     }));

        Ok(vec![
            Update::AddVariable {
                lower: 0.0,
                upper: Real::infinity()
            };
            self.active_constraints.len() - nold
        ])
    }

    fn extend_subgradient(&mut self, fidx: usize, primal: &Self::Primal, inds: &[usize]) -> Result<Vec<Real>> {
        let sub = self.subs[fidx].read().unwrap();
        Ok(inds
            .iter()
            .map(|&i| sub.evaluate_constraint(&primal[0], self.active_constraints[i]))
            .collect())
    }
}

impl ParallelProblem for MMCFProblem {
    type Err = <Self as FirstOrderProblem>::Err;

    type Primal = DVector;

    fn num_variables(&self) -> usize {
        FirstOrderProblem::num_variables(self)
    }

    fn lower_bounds(&self) -> Option<Vec<Real>> {
        FirstOrderProblem::lower_bounds(self)
    }

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

    fn start(&mut self) {

use threadpool::ThreadPool;

use crate::{DVector, Real};

use super::masterprocess::{MasterConfig, MasterProcess, MasterResponse};
use super::problem::{EvalResult, FirstOrderProblem, Update, UpdateState};
use crate::master::{self, MasterProblem};
use crate::solver::{SolverParams, Step};
use crate::terminator::{StandardTerminatable, StandardTerminator, Terminator};
use crate::weighter::{HKWeightable, HKWeighter, Weighter};

/// The default iteration limit.
pub const DEFAULT_ITERATION_LIMIT: usize = 10_000;

/// Error raised by the parallel bundle [`Solver`].

type ClientSender<P> =
    Sender<std::result::Result<EvalResult<usize, <P as FirstOrderProblem>::Primal>, <P as FirstOrderProblem>::Err>>;

type ClientReceiver<P> =
    Receiver<std::result::Result<EvalResult<usize, <P as FirstOrderProblem>::Primal>, <P as FirstOrderProblem>::Err>>;

/// Parameters for tuning the solver.














pub type Parameters = SolverParams;
























pub struct SolverData {
    /// Current center of stability.
    cur_y: DVector,

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

// Copyright (c) 2016, 2017, 2018, 2019 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 crate::{Aggregatable, DVector, Real};
use crate::{Evaluation, FirstOrderProblem, Update};

use crate::master::{self, MasterProblem};
use crate::terminator::{StandardTerminatable, StandardTerminator, Terminator};
use crate::weighter::{HKWeightable, HKWeighter, Weighter};

use log::{debug, info, warn};

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

/// A solver error.
#[derive(Debug)]
pub enum SolverError<E, MErr> {
    /// An error occurred during oracle evaluation.
    Evaluation(E),
    /// An error occurred during oracle update.
    Update(E),
    /// An error has been raised by the master problem.
    BuildMaster(Box<dyn Error>),
    /// An error has been raised by the master problem.
    Master(MErr),
    /// The oracle did not return a minorant.
    NoMinorant,
    /// The dimension of some data is wrong.
    Dimension,
    /// Some parameter has an invalid value.
    Parameter(ParameterError),
    /// The lower bound of a variable is larger than the upper bound.
    InvalidBounds { lower: Real, upper: Real },
    /// The value of a variable is outside its bounds.
    ViolatedBounds { lower: Real, upper: Real, value: Real },
    /// The variable index is out of bounds.
    InvalidVariable { index: usize, nvars: usize },
    /// Iteration limit has been reached.
    IterationLimit { limit: usize },
}

impl<E, MErr> fmt::Display for SolverError<E, MErr>
where
    E: fmt::Display,
    MErr: fmt::Display,
{
    fn fmt(&self, fmt: &mut fmt::Formatter) -> std::result::Result<(), fmt::Error> {
        use self::SolverError::*;
        match self {
            Evaluation(err) => write!(fmt, "Oracle evaluation failed: {}", err),
            Update(err) => write!(fmt, "Oracle update failed: {}", err),
            BuildMaster(err) => write!(fmt, "Creation of master problem failed: {}", err),
            Master(err) => write!(fmt, "Master problem failed: {}", err),
            NoMinorant => write!(fmt, "The oracle did not return a minorant"),
            Dimension => write!(fmt, "Dimension of lower bounds does not match number of variables"),
            Parameter(msg) => write!(fmt, "Parameter error: {}", msg),
            InvalidBounds { lower, upper } => write!(fmt, "Invalid bounds, lower:{}, upper:{}", lower, upper),
            ViolatedBounds { lower, upper, value } => write!(
                fmt,
                "Violated bounds, lower:{}, upper:{}, value:{}",
                lower, upper, value
            ),
            InvalidVariable { index, nvars } => {
                write!(fmt, "Variable index out of bounds, got:{} must be < {}", index, nvars)
            }
            IterationLimit { limit } => write!(fmt, "The iteration limit of {} has been reached.", limit),
        }
    }
}

impl<E, MErr> Error for SolverError<E, MErr>
where
    E: Error + 'static,
    MErr: Error + 'static,
{
    fn source(&self) -> Option<&(dyn Error + 'static)> {
        match self {
            SolverError::Evaluation(err) => Some(err),
            SolverError::Update(err) => Some(err),
            SolverError::Master(err) => Some(err),
            SolverError::BuildMaster(err) => Some(err.as_ref()),
            _ => None,
        }
    }
}

impl<E, MErr> From<MErr> for SolverError<E, MErr> {
    fn from(err: MErr) -> SolverError<E, MErr> {
        SolverError::Master(err)
    }
}

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

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

impl<'a> HKWeightable for BundleState<'a> {
    fn current_weight(&self) -> Real {
        self.weight
    }

    fn center(&self) -> &DVector {
        self.cur_y
    }

    fn center_value(&self) -> Real {
        self.cur_val
    }

    fn candidate_value(&self) -> Real {
        self.nxt_val
    }

    fn candidate_model(&self) -> Real {
        self.nxt_mod
    }

    fn new_cutvalue(&self) -> Real {
        self.new_cutval
    }

    fn sgnorm(&self) -> Real {
        self.sgnorm
    }
}

impl<'a> StandardTerminatable for BundleState<'a> {
    fn expected_progress(&self) -> Real {
        self.expected_progress
    }

    fn center_value(&self) -> Real {
        self.cur_val
    }
}

/// An invalid value for some parameter has been passes.
#[derive(Debug)]
pub struct ParameterError(String);

impl fmt::Display for ParameterError {
    fn fmt(&self, fmt: &mut fmt::Formatter) -> std::result::Result<(), fmt::Error> {
        write!(fmt, "{}", self.0)
    }
}

impl Error for ParameterError {}

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

impl SolverParams {
    /// Verify that all parameters are valid.
    fn check(&self) -> std::result::Result<(), ParameterError> {
        if self.max_bundle_size < 2 {
            Err(ParameterError(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(ParameterError(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(ParameterError(format!(
                "nullstep_factor must be in (0,acceptance_factor] (got: {}, acceptance_factor:{})",
                self.nullstep_factor, self.acceptance_factor
            )))
        } else {
            Ok(())
        }
    }
}

impl Default for SolverParams {
    fn default() -> SolverParams {
        SolverParams {
            max_bundle_size: 50,

            nullstep_factor: 0.1,
            acceptance_factor: 0.1,
        }
    }
}

/// 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 {
    /// The minorant's index in the master problem
    index: usize,
    /// Current multiplier.
    multiplier: Real,
}

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

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

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

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

/// The default builder.
pub type FullMasterBuilder = master::boxed::Builder<master::cpx::Builder>;

/**
 * Implementation of a bundle method.
 */
pub struct Solver<P, T = StandardTerminator, W = HKWeighter, M = FullMasterBuilder>
where
    P: FirstOrderProblem,
    M: master::Builder,
{
    /// The first order problem description.
    problem: P,

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

    /// Termination predicate.
    pub terminator: T,

    /// Weighter heuristic.
    pub weighter: W,

    /// Lower and upper bounds of all variables.
    bounds: Vec<(Real, Real)>,

    /// 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: M::MasterProblem,

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

    /// The primals associated with each global minorant index.
    primals: Vec<Option<P::Primal>>,

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

pub type Result<T, P, M> = std::result::Result<
    T,
    SolverError<<P as FirstOrderProblem>::Err, <<M as master::Builder>::MasterProblem as MasterProblem>::Err>,
>;

impl<P, T, W, M> Solver<P, T, W, M>
where
    P: FirstOrderProblem,
    P::Err: Into<Box<dyn std::error::Error + Sync + Send>>,
    T: for<'a> Terminator<BundleState<'a>> + Default,
    W: for<'a> Weighter<BundleState<'a>> + Default,
    M: master::Builder + Default,
    M::MasterProblem: MasterProblem<MinorantIndex = usize>,
{
    /**
     * 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()`.
     */
    #[allow(clippy::type_complexity)]
    pub fn new_params(problem: P, params: SolverParams) -> Result<Solver<P, T, W, M>, P, M> {
        Ok(Solver {
            problem,
            params,
            terminator: T::default(),
            weighter: W::default(),
            bounds: vec![],
            cur_y: dvec![],
            cur_val: 0.0,
            cur_mod: 0.0,
            cur_vals: dvec![],
            cur_mods: dvec![],
            cur_valid: false,
            nxt_d: dvec![],
            nxt_y: dvec![],
            nxt_val: 0.0,
            nxt_mod: 0.0,
            nxt_vals: dvec![],
            nxt_mods: dvec![],
            new_cutval: 0.0,
            sgnorm: 0.0,
            expected_progress: 0.0,
            cnt_descent: 0,
            cnt_null: 0,
            start_time: Instant::now(),
            master: M::default().build().map_err(|e| SolverError::BuildMaster(e.into()))?,
            minorants: vec![],
            primals: vec![],
            iterinfos: vec![],
        })
    }

    /// A new solver with default parameter.
    #[allow(clippy::type_complexity)]
    pub fn new(problem: P) -> Result<Solver<P, T, W, M>, P, M> {
        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<(), P, M> {
        self.params.check().map_err(SolverError::Parameter)?;
        if self.cur_y.len() != self.problem.num_variables() {
            self.cur_valid = false;
            self.cur_y.init0(self.problem.num_variables());
        }

        let lb = self.problem.lower_bounds();
        let ub = self.problem.upper_bounds();
        self.bounds.clear();
        self.bounds.reserve(self.cur_y.len());
        for i in 0..self.cur_y.len() {
            let lb_i = lb.as_ref().map(|x| x[i]).unwrap_or(NEG_INFINITY);
            let ub_i = ub.as_ref().map(|x| x[i]).unwrap_or(INFINITY);
            if lb_i > ub_i {
                return Err(SolverError::InvalidBounds {
                    lower: lb_i,
                    upper: ub_i,
                });
            }
            if self.cur_y[i] < lb_i {
                self.cur_valid = false;
                self.cur_y[i] = lb_i;
            } else if self.cur_y[i] > ub_i {
                self.cur_valid = false;
                self.cur_y[i] = ub_i;
            }
            self.bounds.push((lb_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 with at most 10_000 iterations.
    ///
    /// Use `solve_with_limit` for an explicit iteration limit.
    pub fn solve(&mut self) -> Result<(), P, M> {
        const LIMIT: usize = 10_000;
        self.solve_with_limit(LIMIT)
    }

    /// Solve the problem with explicit iteration limit.
    pub fn solve_with_limit(&mut self, iter_limit: usize) -> Result<(), P, M> {
        // First initialize the internal data structures.
        self.init()?;

        if self.solve_iter(iter_limit)? {
            Ok(())
        } else {
            Err(SolverError::IterationLimit { limit: iter_limit })
        }
    }

    /// Solve the problem but stop after `niter` iterations.
    ///
    /// The function returns `Ok(true)` if the termination criterion
    /// has been satisfied. Otherwise it returns `Ok(false)` or an
    /// error code.
    ///
    /// If this function is called again, the solution process is
    /// continued from the previous point. Because of this one must
    /// call `init()` before the first call to this function.
    pub fn solve_iter(&mut self, niter: usize) -> Result<bool, P, M> {
        for _ in 0..niter {
            let mut term = self.step()?;
            let changed = self.update_problem(term)?;
            // do not stop if the problem has been changed
            if changed && term == Step::Term {
                term = Step::Null
            }
            self.show_info(term);
            if term == Step::Term {
                return Ok(true);
            }
        }
        Ok(false)
    }

    /// Called to update the problem.
    ///
    /// 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, P, M> {
        let updates = {
            let state = UpdateState {
                minorants: &self.minorants,
                primals: &self.primals,
                step: term,
                iteration_info: &self.iterinfos,
                // this is a dirty trick: when updating the center, we
                // simply swapped the `cur_*` fields with the `nxt_*`
                // fields
                cur_y: if term == Step::Descent {
                    &self.nxt_y
                } else {
                    &self.cur_y
                },
                nxt_y: if term == Step::Descent {
                    &self.cur_y
                } else {
                    &self.nxt_y
                },
            };
            self.problem.update(&state).map_err(SolverError::Update)?
        };

        let mut newvars = Vec::with_capacity(updates.len());
        for u in updates {
            match u {
                Update::AddVariable { lower, upper } => {
                    if lower > upper {
                        return Err(SolverError::InvalidBounds { lower, upper });
                    }
                    let value = if lower > 0.0 {
                        lower
                    } else if upper < 0.0 {
                        upper
                    } else {
                        0.0
                    };
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, value));
                }
                Update::AddVariableValue { lower, upper, value } => {
                    if lower > upper {
                        return Err(SolverError::InvalidBounds { lower, upper });
                    }
                    if value < lower || value > upper {
                        return Err(SolverError::ViolatedBounds { lower, upper, value });
                    }
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, value));
                }
                Update::MoveVariable { index, value } => {
                    if index >= self.bounds.len() {
                        return Err(SolverError::InvalidVariable {
                            index,
                            nvars: self.bounds.len(),
                        });
                    }
                    let (lower, upper) = self.bounds[index];
                    if value < lower || value > upper {
                        return Err(SolverError::ViolatedBounds { lower, upper, value });
                    }
                    newvars.push((Some(index), lower - value, upper - value, value));
                }
            }
        }

        if !newvars.is_empty() {
            let problem = &mut self.problem;
            let primals = &self.primals;
            self.master.add_vars(
                &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::<Vec<_>>(),
                &mut |fidx, minidx, vars| {
                    problem
                        .extend_subgradient(fidx, primals[minidx].as_ref().unwrap(), vars)
                        .map(DVector)
                        .map_err(|e| e.into())
                },
            )?;
            // modify moved variables
            for (index, val) in newvars.iter().filter_map(|v| v.0.map(|i| (i, v.3))) {
                self.cur_y[index] = val;
                self.nxt_y[index] = val;
                self.nxt_d[index] = 0.0;
            }
            // add new variables
            self.cur_y.extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3));
            self.nxt_y.extend(newvars.iter().filter(|v| v.0.is_none()).map(|v| v.3));
            self.nxt_d.resize(self.nxt_y.len(), 0.0);
            Ok(true)
        } else {
            Ok(false)
        }
    }

    /// 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) -> P::Primal {
        Aggregatable::combine(
            self.minorants[subproblem]
                .iter()
                .map(|m| (m.multiplier, self.primals[m.index].as_ref().unwrap())),
        )
    }

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

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

    /// Return the last candidate point.
    pub fn candidate(&self) -> &[Real] {
        &self.nxt_y
    }

    /**
     * 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<(), P, M> {
        let m = self.problem.num_subproblems();

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

        if lb
            .as_ref()
            .map(|lb| lb.len() != self.problem.num_variables())
            .unwrap_or(false)
        {
            return Err(SolverError::Dimension);
        }
        if ub
            .as_ref()
            .map(|ub| ub.len() != self.problem.num_variables())
            .unwrap_or(false)
        {
            return Err(SolverError::Dimension);
        }

        self.master.set_num_subproblems(m)?;
        self.master.set_vars(self.problem.num_variables(), lb, ub)?;

        self.minorants = (0..m).map(|_| vec![]).collect();

        self.cur_val = 0.0;
        for i in 0..m {
            let result = self
                .problem
                .evaluate(i, &self.cur_y, INFINITY, 0.0)
                .map_err(SolverError::Evaluation)?;
            self.cur_vals[i] = result.objective();
            self.cur_val += self.cur_vals[i];

            let mut minorants = result.into_iter();
            if let Some((minorant, primal)) = minorants.next() {
                self.cur_mods[i] = minorant.constant;
                self.cur_mod += self.cur_mods[i];
                let minidx = self.master.add_minorant(i, minorant)?;
                self.minorants[i].push(MinorantInfo {
                    index: minidx,
                    multiplier: 0.0,
                });
                if minidx >= self.primals.len() {
                    self.primals.resize_with(minidx + 1, || None);
                }
                self.primals[minidx] = Some(primal);
            } else {
                return Err(SolverError::NoMinorant);
            }
        }

        self.cur_valid = true;

        // Solve the master problem once to compute the initial
        // subgradient.
        //
        // We could compute that subgradient directly by
        // adding up the initial minorants, but this would not include
        // the eta terms. However, this is a heuristic anyway because
        // we assume an initial weight of 1.0, which, in general, will
        // *not* be the initial weight for the first iteration.
        self.master.set_weight(1.0)?;
        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.initial_weight(&state);
        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<(), P, M> {
        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 &mut self.minorants[i] {
                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<(), P, M> {
        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, self.primals[m.index].take().unwrap()))
                    .unzip();
                let (aggr_min, aggr_coeffs) = self.master.aggregate(i, &aggr_mins)?;
                // append aggregated minorant
                self.minorants[i].push(MinorantInfo {
                    index: aggr_min,
                    multiplier: aggr_sum,
                });
                self.primals[aggr_min] = Some(Aggregatable::combine(aggr_coeffs.into_iter().zip(&aggr_primals)));
            }
        }
        Ok(())
    }

    /// Perform a descent step.
    fn descent_step(&mut self) -> Result<(), P, M> {
        let new_weight = self.weighter.descent_weight(&current_state!(self, Step::Descent));
        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);
        Ok(())
    }

    /// Perform a null step.
    fn null_step(&mut self) -> Result<(), P, M> {
        let new_weight = self.weighter.null_weight(&current_state!(self, Step::Null));
        self.master.set_weight(new_weight)?;
        self.cnt_null += 1;
        debug!("Null Step");
        Ok(())
    }

    /// Perform one bundle iteration.
    #[allow(clippy::collapsible_if)]
    pub fn step(&mut self) -> Result<Step, P, M> {
        self.iterinfos.clear();

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

        self.solve_model()?;
        if self.terminator.terminate(&current_state!(self, Step::Term)) {
            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;
        self.new_cutval = 0.0;
        for fidx in 0..self.problem.num_subproblems() {
            let result = self
                .problem
                .evaluate(fidx, &self.nxt_y, nullstep_bnd, relprec)
                .map_err(SolverError::Evaluation)?;
            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(SolverError::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;
            let minidx = self.master.add_minorant(fidx, nxt_minorant)?;
            self.minorants[fidx].push(MinorantInfo {
                index: minidx,
                multiplier: 0.0,
            });
            if minidx >= self.primals.len() {
                self.primals.resize_with(minidx + 1, || None);
            }
            self.primals[minidx] = Some(nxt_primal);
        }

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

        self.nxt_val = nxt_ub;

        // check for potential problems with relative precision of all kinds
        if nxt_lb <= descent_bnd {
            // lower bound gives descent step
            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.");
                    self.iterinfos.push(IterationInfo::UpperBoundNullStep);
                }
            }
        } else if self.cur_val - nxt_lb > 0.8 * (self.cur_val - self.nxt_mod) {
            // TODO: double check with ConicBundle if this test makes sense.
            // lower bound gives already a null step
            // subgradient won't yield much improvement
            warn!("Shallow cut (subgradient won't yield much improvement)");
            self.iterinfos.push(IterationInfo::ShallowCut);
        }

        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()?;
            Ok(Step::Descent)
        } else {
            self.null_step()?;
            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)
    }
}