RsBundle  Diff

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

}

impl<E> From<ParameterError> for SolverError<E> {
    fn from(err: ParameterError) -> SolverError<E> {
        SolverError::Parameter(err)
    }
}

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

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

/**
 * Implementation of a bundle method.
 */
pub struct Solver<P: FirstOrderProblem> {

     */
    start_time: Instant,

    /// The master problem.
    master: Box<dyn MasterProblem<MinorantIndex = usize>>,

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

impl<P: FirstOrderProblem> Solver<P>
where

            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::new(MinimalMaster::new()?)),


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

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

        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<(), SolverError<P::Err>> {
        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<(), SolverError<P::Err>> {
        // 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

    ///
    /// 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<P::Err>> {
        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

                    if subproblem >= self.minorants.len() {
                        return Err(SolverError::InvalidSubproblem {
                            subproblem: subproblem,
                            nsubs: self.minorants.len(),
                        });
                    }
                    for m in &mut self.minorants[subproblem] {
                        if let Some(ref mut p) = self.primals[m.index] {
                            if let Err(err) = modify(p) {
                                return Err(SolverError::Update(err));
                            }
                        }
                    }
                }
            }
        }

        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

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

     * information.
     */
    fn init_master(&mut self) -> Result<(), SolverError<P::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()?))


        } else {
            debug!("Use CPLEX master problem");
            Box::new(BoxedMasterProblem::new(CplexMaster::new()?))


        };

        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.master.set_max_updates(self.params.max_updates)?;


        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.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<(), SolverError<P::Err>> {
        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() {
            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(
                    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<P::Err>> {
        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);
        Ok(())
    }

    /// Perform a null step.
    fn null_step(&mut self) -> Result<(), SolverError<P::Err>> {
        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");
        Ok(())
    }

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

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

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