RsBundle  Check-in [b0a2986e97]

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

use failure::Error;

/// A solver error.
#[derive(Debug, Fail)]
pub enum SolverError {
    /// An error occured during oracle evaluation.
    #[fail(display = "Oracle evaluation failed: {}", _0)]
    Evaluation(Error),
    /// An error occured during oracle update.
    #[fail(display = "Oracle update failed: {}", _0)]
    Update(Error),
    /// An error has been raised by the master problem.
    #[fail(display = "Master problem failed: {}", _0)]
    Master(Error),
    /// The oracle did not return a minorant.
    #[fail(display = "The oracle did not return a minorant")]
    NoMinorant,
    /// The dimension of some data is wrong.
    #[fail(display = "Dimension of lower bounds does not match number of variables")]
    Dimension,
    /// Some parameter has an invalid value.

     *
     * 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>, SolverError>
    {
        Ok(Solver {
            problem: problem,
            params: params,
            terminator: Box::new(StandardTerminator { termination_precision: 1e-3 }),
            weighter: Box::new(HKWeighter::new()),
            bounds: vec![],

            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: Box::new(BoxedMasterProblem::<MinimalMaster>::new().map_err(SolverError::Master)?),
            minorants: vec![],
            iterinfos: vec![],
        })
    }

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

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

        self.start_time = Instant::now();

        Ok(())
    }

    /// Solve the problem.
    pub fn solve(&mut self) -> Result<(), SolverError> {
        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> {
        for _ in 0..niter {
            let mut term = try!(self.step());
            let changed = try!(self.update_problem(term));
            // do not stop if the problem has been changed
            if changed && term == Step::Term {
                term = Step::Null
            }
            self.show_info(term);
            if term == Step::Term {
                return Ok(true)
            }
        }
        Ok(false)
    }

    /// Called to update the problem.
    ///
    /// Calling this function typically triggers the problem to
    /// separate new constraints depending on the current solution.
    fn update_problem(&mut self, term: Step) -> Result<bool, SolverError> {
        let updates = {
            let state = UpdateState {
                minorants: &self.minorants,
                step: term,
                iteration_info: &self.iterinfos,
                // this is a dirty trick: when updating the center, we
                // simply swapped the `cur_*` fields with the `nxt_*`
                // fields
                cur_y: if term == Step::Descent {
                    &self.nxt_y
                } else {
                    &self.cur_y
                },
                nxt_y: if term == Step::Descent {
                    &self.cur_y
                } else {
                    &self.nxt_y
                },
            };
            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 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| {
                                     problem.extend_subgradient(minorants[fidx][minidx].primal.as_ref().unwrap(), vars)
                                         .map(DVector)
                                         .unwrap()
                                 }).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;
                self.nxt_d[index] = 0.0;
            }
            // add new variables

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

        self.master = if m == 1 && self.params.max_bundle_size == 2 {
            debug!("Use minimal master problem");
            Box::new(BoxedMasterProblem::<MinimalMaster>::new().unwrap())
        } else {
            debug!("Use CPLEX master problem");
            Box::new(BoxedMasterProblem::<CplexMaster>::new().unwrap())
        };

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

        if let Some(ref x) = lb {
            if x.len() != self.problem.num_variables() {
                return Err(SolverError::Dimension);
            }
        }

        self.master.set_num_subproblems(m).map_err(SolverError::Master)?;
        self.master.set_vars(self.problem.num_variables(), lb, ub).map_err(SolverError::Master)?;
        self.master.set_max_updates(self.params.max_updates).map_err(SolverError::Master)?;

        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 = self.problem.evaluate(i, &self.cur_y, INFINITY, 0.0).map_err(SolverError::Master)?;
            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: self.master.add_minorant(i, minorant).map_err(SolverError::Master)?,
                    multiplier: 0.0,
                    primal: 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).map_err(SolverError::Master)?;
        self.master.solve(self.cur_val).map_err(SolverError::Master)?;
        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).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> {
        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;

        // update multiplier from master solution

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


    /// Reduce size of bundle.
    fn compress_bundle(&mut self) -> Result<(), SolverError> {
        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) = self.master.aggregate(i, &aggr_mins).map_err(SolverError::Master)?;
                // append aggregated minorant
                self.minorants[i].push(MinorantInfo {
                    index: aggr_min,
                    multiplier: aggr_sum,
                    primal: Some(self.problem.aggregate_primals(aggr_coeffs.into_iter()
                        .zip(aggr_primals.into_iter())
                        .collect())),
                });
            }
        }
        Ok(())
    }

    /// Perform a descent step.
    fn descent_step(&mut self) -> Result<(), SolverError> {
        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> {
        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> {
        self.iterinfos.clear();

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

        try!(self.compress_bundle());

        let mut nxt_lb = 0.0;
        let mut nxt_ub = 0.0;
        self.new_cutval = 0.0;
        for fidx in 0..self.problem.num_subproblems() {
            let result = 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;
            self.minorants[fidx].push(MinorantInfo {
                index: self.master.add_minorant(fidx, nxt_minorant).map_err(SolverError::Master)?,
                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:{}",