RsBundle  Diff

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! The main bundle method solver.

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

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

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

use failure::Error;

    #[fail(display = "Parameter error: {}", _0)]
    Parameter(String),
    /// The lower bound of a variable is larger than the upper bound.
    #[fail(display = "Invalid bounds, lower:{} upper:{}", lower, upper)]
    InvalidBounds { lower: Real, upper: Real },
    /// The value of a variable is outside its bounds.
    #[fail(display = "Violated bounds, lower:{} upper:{} value:{}", lower, upper, value)]
    ViolatedBounds { lower: Real, upper: Real, value: Real },
    ViolatedBounds {
        lower: Real,
        upper: Real,
        value: Real,
    },
    /// The variable index is out of bounds.
    #[fail(display = "Variable index out of bounds, got:{} must be < {}", index, nvars)]
    InvalidVariable { index: usize, nvars: usize },
    /// Iteration limit has been reached.
    #[fail(display = "The iteration limit of {} has been reached.", limit)]
    IterationLimit { limit: usize },
}

    pub max_updates: usize,
}

impl SolverParams {
    /// Verify that all parameters are valid.
    fn check(&self) -> Result<(), SolverError> {
        if self.max_bundle_size < 2 {
            Err(SolverError::Parameter(format!("max_bundle_size must be >= 2 (got: {})",
                                               self.max_bundle_size)))
            Err(SolverError::Parameter(format!(
                "max_bundle_size must be >= 2 (got: {})",
                self.max_bundle_size
            )))
        } else if self.acceptance_factor <= 0.0 || self.acceptance_factor >= 1.0 {
            Err(SolverError::Parameter(format!("acceptance_factor must be in (0,1) (got: {})",
                                               self.acceptance_factor)))
            Err(SolverError::Parameter(format!(
                "acceptance_factor must be in (0,1) (got: {})",
                self.acceptance_factor
            )))
        } else if self.nullstep_factor <= 0.0 || self.nullstep_factor > self.acceptance_factor {
            Err(SolverError::Parameter(
                format!("nullstep_factor must be in (0,acceptance_factor] (got: {}, acceptance_factor:{})",
                        self.nullstep_factor, self.acceptance_factor)))
            Err(SolverError::Parameter(format!(
                "nullstep_factor must be in (0,acceptance_factor] (got: {}, acceptance_factor:{})",
                self.nullstep_factor,
                self.acceptance_factor
            )))
        } else if self.min_weight <= 0.0 {
            Err(SolverError::Parameter(format!("min_weight must be in > 0 (got: {})", self.min_weight)))
            Err(SolverError::Parameter(format!(
                "min_weight must be in > 0 (got: {})",
                self.min_weight
            )))
        } else if self.max_weight < self.min_weight {
            Err(SolverError::Parameter(format!(
            Err(SolverError::Parameter(format!("max_weight must be in >= min_weight (got: {}, min_weight: {})",
                                               self.max_weight, self.min_weight)))
                "max_weight must be in >= min_weight (got: {}, min_weight: {})",
                self.max_weight,
                self.min_weight
            )))
        } else if self.max_updates == 0 {
            Err(SolverError::Parameter(format!("max_updates must be in > 0 (got: {})", self.max_updates)))
            Err(SolverError::Parameter(format!(
                "max_updates must be in > 0 (got: {})",
                self.max_updates
            )))
        } else {
            Ok(())
        }
    }
}

impl Default for SolverParams {

    }
}

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

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

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


impl<P, Pr, E> Solver<P, Pr, E>
where
    where P: for<'a> FirstOrderProblem<'a, Primal = Pr, EvalResult = E>,
          E: Evaluation<Pr>
    P: for<'a> FirstOrderProblem<'a, Primal = Pr, EvalResult = E>,
    E: Evaluation<Pr>,
{
    /**
     * Create a new solver for the given problem.
     *
     * Note that the solver owns the problem, so you cannot use the
     * same problem description elsewhere as long as it is assigned to
     * the solver. However, it is possible to get a reference to the
     * internally stored problem using `Solver::problem()`.
     */
    pub fn new_params(problem: P, params: SolverParams)
    pub fn new_params(problem: P, params: SolverParams) -> Result<Solver<P, Pr, E>, SolverError> {
                      -> Result<Solver<P, Pr, E>, SolverError>
    {
        Ok(Solver {
            problem: problem,
            params: params,
            terminator: Box::new(StandardTerminator { termination_precision: 1e-3 }),
            terminator: Box::new(StandardTerminator {
                termination_precision: 1e-3,
            }),
            weighter: Box::new(HKWeighter::new()),
            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: Box::new(BoxedMasterProblem::new(MinimalMaster::new().map_err(SolverError::Master)?)),
            master: Box::new(BoxedMasterProblem::new(
                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> {

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

            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)
                return Ok(true);
            }
        }
        Ok(false)
    }

    /// Called to update the problem.
    ///

                        upper
                    } else {
                        0.0
                    };
                    self.bounds.push((lower, upper));
                    newvars.push((None, lower - value, upper - value, value));
                }
                Update::AddVariableValue { lower, upper, 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 });
                        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,
                            index, nvars: self.bounds.len()
                            nvars: self.bounds.len(),
                        });
                    }
                    let (lower, upper) = self.bounds[index];
                    if value < lower || value > upper {
                        return Err(SolverError::ViolatedBounds { lower, upper, value });
                        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(
            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)?;
                    &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
            self.cur_y
            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));
                .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)
        }
    }

            .iter()
            .map(|m| (m.multiplier, m.primal.as_ref().unwrap()))
            .collect()
    }

    fn show_info(&self, step: Step) {
        let time = self.start_time.elapsed();
        info!(
        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);
            "{} {: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
    }

     * 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::new(MinimalMaster::new().map_err(SolverError::Master)?))
            Box::new(BoxedMasterProblem::new(
                MinimalMaster::new().map_err(SolverError::Master)?,
            ))
        } else {
            debug!("Use CPLEX master problem");
            Box::new(BoxedMasterProblem::new(CplexMaster::new().map_err(SolverError::Master)?))
            Box::new(BoxedMasterProblem::new(
                CplexMaster::new().map_err(SolverError::Master)?,
            ))
        };

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

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

        self.master
            .set_num_subproblems(m)
        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)?;
            .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 = (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];
                self.minorants[i].push(MinorantInfo {
                    index: self.master.add_minorant(i, minorant).map_err(SolverError::Master)?,
                    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)
        self.master.solve(self.cur_val).map_err(SolverError::Master)?;
            .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)
        self.master.set_weight(new_weight).map_err(SolverError::Master)?;
            .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)
        self.master.solve(self.cur_val).map_err(SolverError::Master)?;
            .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

                // 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)?;
                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(
                    primal: Some(self.problem.aggregate_primals(aggr_coeffs.into_iter()
                        .zip(aggr_primals.into_iter())
                        .collect())),
                        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
        let new_weight = self.weighter.weight(&current_state!(self, Step::Descent), &self.params);
        self.master.set_weight(new_weight).map_err(SolverError::Master)?;
            .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
        let new_weight = self.weighter.weight(&current_state!(self, Step::Null), &self.params);
        self.master.set_weight(new_weight).map_err(SolverError::Master)?;
            .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
            self.init_master()?;
        }

        self.solve_model()?;
        if self.terminator
        if self.terminator.terminate(&current_state!(self, Step::Term), &self.params) {
            .terminate(&current_state!(self, Step::Term), &self.params)
        {
            return Ok(Step::Term);
        }

        let m = self.problem.num_subproblems();
        let descent_bnd = self.get_descent_bound();
        let nullstep_bnd = if m == 1 {
            self.get_nullstep_bound()

            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
                index: self.master.add_minorant(fidx, nxt_minorant).map_err(SolverError::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!(
            warn!("New minorant has higher value in center new:{} old:{}",
                  self.new_cutval,
                  self.cur_val);
                "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,
            });
        }