RsBundle  Diff

//

//! The main bundle method solver.

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

use master::{BoxedMasterProblem, Error as MasterProblemError, MasterProblem, UnconstrainedMasterProblem};
use master::{CplexMaster, MinimalMaster};

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



/// A solver error.
#[derive(Debug)]
pub enum SolverError<E> {
    /// An error occured during oracle evaluation.

    Evaluation(E),
    /// An error occured during oracle update.

    Update(E),
    /// An error has been raised by the master problem.

    Master(MasterProblemError),
    /// 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: fmt::Display> fmt::Display for SolverError<E> {
    fn fmt(&self, fmt: &mut fmt::Formatter) -> Result<(), fmt::Error> {
        use self::SolverError::*;
        match self {
            Evaluation(err) => write!(fmt, "Oracle evaluation failed: {}", err),
            Update(err) => write!(fmt, "Oracle update 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: Error> Error for SolverError<E> {
    fn cause(&self) -> Option<&Error> {
        match self {
            SolverError::Evaluation(err) => Some(err),
            SolverError::Update(err) => Some(err),
            SolverError::Master(err) => Some(err.as_ref()),
            _ => None,
        }
    }
}

impl<E> From<ParameterError> for SolverError<E> {
    fn from(err: ParameterError) -> SolverError<E> {
        SolverError::Parameter(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

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

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

     * variables.
     */
    pub max_updates: usize,
}

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

        self.minorants[fidx].last().and_then(|m| m.primal.as_ref())
    }
}

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

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

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

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

impl<P, Pr, E, Err> Solver<P, Pr, E, Err>
where
    P: for<'a> FirstOrderProblem<'a, Primal = Pr, EvalResult = E, Err = Err>,
    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, Err>, SolverError<Err>> {
        Ok(Solver {
            problem,
            params,
            terminator: Box::new(StandardTerminator {
                termination_precision: 1e-3,
            }),
            weighter: Box::new(HKWeighter::new()),

            )),
            minorants: vec![],
            iterinfos: vec![],
        })
    }

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

    /**
     * Set the first order problem description associated with this
     * solver.
     *

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

    /// Initialize the solver.
    pub fn init(&mut self) -> Result<(), SolverError<Err>> {
        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();

        self.start_time = Instant::now();

        Ok(())
    }

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

        if self.solve_iter(LIMIT)? {
            Ok(())
        } else {
            Err(SolverError::IterationLimit { limit: 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, SolverError<Err>> {
        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, SolverError<Err>> {
        let updates = {
            let state = UpdateState {
                minorants: &self.minorants,
                step: term,
                iteration_info: &self.iterinfos,
                // this is a dirty trick: when updating the center, we
                // simply swapped the `cur_*` fields with the `nxt_*`

        if !newvars.is_empty() {
            let mut problem = &mut self.problem;
            let minorants = &self.minorants;
            self.master
                .add_vars(
                    &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::<Vec<_>>(),
                    &mut move |fidx, minidx, vars| match problem

                        .extend_subgradient(minorants[fidx][minidx].primal.as_ref().unwrap(), vars)
                        .map(DVector)
                    {
                        Ok(g) => g,
                        Err(_) => unreachable!(),
                    },
                )
                .map_err(SolverError::Master)?;
            // 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;

    /**
     * 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<(), SolverError<Err>> {
        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::new(
                MinimalMaster::new().map_err(SolverError::Master)?,
            ))

        debug!("Init master completed");

        Ok(())
    }

    /// Solve the model (i.e. master problem) to compute the next candidate.
    fn solve_model(&mut self) -> Result<(), SolverError<Err>> {
        self.master.solve(self.cur_val).map_err(SolverError::Master)?;
        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!("  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<(), SolverError<Err>> {
        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();

                });
            }
        }
        Ok(())
    }

    /// Perform a descent step.
    fn descent_step(&mut self) -> Result<(), SolverError<Err>> {
        let new_weight = self.weighter.weight(&current_state!(self, Step::Descent), &self.params);
        self.master.set_weight(new_weight).map_err(SolverError::Master)?;
        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<(), SolverError<Err>> {
        let new_weight = self.weighter.weight(&current_state!(self, Step::Null), &self.params);
        self.master.set_weight(new_weight).map_err(SolverError::Master)?;
        self.cnt_null += 1;
        debug!("Null Step");
        Ok(())
    }

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

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