Index: src/master/boxed/unconstrained/cpx.rs ================================================================== --- src/master/boxed/unconstrained/cpx.rs +++ src/master/boxed/unconstrained/cpx.rs @@ -25,11 +25,10 @@ use cplex_sys as cpx; use cplex_sys::trycpx; use either::Either; use log::{debug, warn}; -use std; use std::collections::VecDeque; use std::f64::NEG_INFINITY; use std::iter::{once, repeat}; use std::ops::{Deref, DerefMut}; use std::os::raw::{c_char, c_int}; Index: src/mcf/solver.rs ================================================================== --- src/mcf/solver.rs +++ src/mcf/solver.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer +// Copyright (c) 2016, 2017, 2018, 2019, 2020 Frank Fischer // // 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. @@ -20,11 +20,10 @@ use c_str_macro::c_str; use cplex_sys as cpx; use cplex_sys::trycpx; -use std; use std::ffi::CString; use std::ptr; use std::result; use std::os::raw::{c_char, c_double, c_int}; Index: src/solver/asyn.rs ================================================================== --- src/solver/asyn.rs +++ src/solver/asyn.rs @@ -22,12 +22,11 @@ #[cfg(not(feature = "crossbeam"))] use std::sync::mpsc::{channel, RecvError}; use float_pretty_print::PrettyPrintFloat; use log::{debug, info, warn}; -use num_cpus; -use num_traits::{Float, ToPrimitive}; +use num_traits::{Float, ToPrimitive, Zero}; use std::iter::repeat; use std::sync::Arc; use std::time::Instant; use threadpool::ThreadPool; @@ -41,12 +40,12 @@ use crate::master::{Builder as MasterBuilder, MasterProblem}; use crate::problem::{FirstOrderProblem, UpdateState}; use crate::terminator::{StandardTerminatable, StandardTerminator, Terminator}; use crate::weighter::{HKWeightable, HKWeighter, Weighter}; -mod subzero; -use subzero::SubZero; +pub mod guessmodels; +use guessmodels::{Guess, GuessModel, NearestValue}; /// The default iteration limit. pub const DEFAULT_ITERATION_LIMIT: usize = 10_000; /// The default solver. @@ -187,10 +186,41 @@ /// The index of the subproblem. subproblem: usize, /// The index of the candidate at which the subproblem is evaluated. candidate_index: usize, } + +/// An evaluation point. +#[derive(Clone)] +pub struct Point { + /// The globally unique index of the evaluation point. + index: usize, + + /// The evaluation point itself. + point: Arc, +} + +impl Point { + fn distance(&self, p: &Point) -> Real { + if self.index != p.index { + let mut d = self.point.as_ref().clone(); + d.add_scaled(-1.0, &p.point); + d.norm2() + } else { + Real::zero() + } + } +} + +impl Default for Point { + fn default() -> Point { + Point { + index: 0, + point: Arc::new(dvec![]), + } + } +} /// Parameters for tuning the solver. #[derive(Debug, Clone)] pub struct Parameters { /// The descent step acceptance factors, must be in (0,1). @@ -250,14 +280,20 @@ Term, } pub struct SolverData { /// Current center of stability. - cur_y: Arc, + cur_y: Point, /// Function value in the current point. cur_val: Real, + + /// Step direction (i.e. nxt_y - cur_y). + nxt_d: Arc, + + /// Current candidate. + nxt_y: Point, /// Function value at the current candidate. nxt_val: Real, /// Model value at the current candidate. @@ -290,30 +326,24 @@ cnt_updates: usize, /// Number of descent steps. cnt_descent: usize, - /// Index of the current center. - center_index: usize, - - /// Index of the current candidate. - candidate_index: usize, - /// The number of subproblems with insufficient evaluation data. num_insufficient_candidates: usize, + + /// The number of subproblems that have not been evaluated exactly in the center. + num_inexact_center: usize, + + /// Whether the next step should be a forced descent step. + force_descent: bool, /// Subproblem data. subs: Vec, - /// Step direction (i.e. nxt_y - cur_y). - nxt_d: Arc, - - /// Current candidate. - nxt_y: Arc, - /// The list of all evaluation points. - candidates: Vec, + candidates: Vec, /// Whether we need a new update need_update: bool, /// Whether a problem update is currently in progress. @@ -323,51 +353,53 @@ impl SolverData { /// Reset solver data to initial values. /// /// This means that almost everything is set to +infinity so that /// a null-step is forced after the first evaluation. - fn init(&mut self, y: DVector) { + fn init(&mut self, y: Point) { self.cnt_descent = 0; - self.cur_y = Arc::new(y); + self.cur_y = y.clone(); self.cur_val = Real::infinity(); + self.nxt_y = y; self.nxt_val = Real::infinity(); self.nxt_mod = -Real::infinity(); self.nxt_submods = vec![-Real::infinity(); self.nxt_submods.len()]; self.expected_progress = Real::infinity(); self.error_bound = Real::infinity(); self.candidates.clear(); - self.candidate_index = 0; self.sgnorm = Real::infinity(); self.cur_weight = 1.0; self.num_insufficient_candidates = 0; + self.num_inexact_center = self.nxt_submods.len(); + self.force_descent = false; } } impl Default for SolverData { fn default() -> SolverData { SolverData { - cur_y: Arc::new(dvec![]), + cur_y: Point::default(), cur_val: 0.0, nxt_val: 0.0, nxt_mod: 0.0, nxt_submods: vec![], expected_progress: 0.0, sgnorm: 0.0, error_bound: Real::infinity(), cur_weight: 1.0, - center_index: 0, - candidate_index: 0, num_insufficient_candidates: 0, + num_inexact_center: 0, + force_descent: false, candidates: vec![], max_iter: 0, cnt_descent: 0, cnt_updates: 0, subs: vec![], nxt_d: Arc::new(dvec![]), - nxt_y: Arc::new(dvec![]), + nxt_y: Point::default(), need_update: true, update_in_progress: false, } } } @@ -386,11 +418,11 @@ fn current_weight(&self) -> Real { self.cur_weight } fn center(&self) -> &DVector { - &self.cur_y + &self.cur_y.point } fn center_value(&self) -> Real { self.cur_val } @@ -446,163 +478,10 @@ fn aggregated_primal(&self, i: usize) -> &Pr { &self.primal_aggrs[i] } } -/// Information about an evaluation point. -#[derive(Clone)] -struct EvalPoint { - /// The index of the evaluation point. - index: usize, - /// The evaluation point itself. - point: Arc, - - /// The current center point. - center: Arc, - - /// The direction from the current center to the current candidate. - nxt_d: Arc, - - /// The direction from the current center to this evaluation point. - center_d: DVector, - - /// The l2-distance to the current candidate. - candidate_dist: Real, - - /// The index of the current center. - center_index: usize, - - /// The index of the current candidate. - candidate_index: usize, -} - -impl EvalPoint { - fn new( - index: usize, - y: Arc, - center_index: usize, - center: &Arc, - candidate_index: usize, - candidate: &DVector, - nxt_d: &Arc, - ) -> EvalPoint { - // Initialize evaluation point with some dummy data. - let mut p = EvalPoint { - index, - point: y, - center: center.clone(), - nxt_d: nxt_d.clone(), - center_d: DVector::default(), - candidate_dist: 0.0, - center_index: center_index + 1, // +1 ensures that the call to `update` will take effect - candidate_index: candidate_index + 1, - }; - // Initialize the distances/directions. - p.update(center_index, center, candidate_index, candidate, nxt_d); - p - } - - /// Possibly data if center or candidate changes. - fn update( - &mut self, - center_index: usize, - center: &Arc, - candidate_index: usize, - candidate: &DVector, - nxt_d: &Arc, - ) { - // candidate changed -> update candidate_dist = |y - y_cand| - if self.candidate_index != candidate_index { - let mut cand_d = self.point.as_ref().clone(); - cand_d.add_scaled(-1.0, candidate); - self.candidate_dist = cand_d.norm2(); - } - - // center changed -> update direction center_d = y - y_center - if self.center_index != center_index { - self.center_d = self.point.as_ref().clone(); - self.center_d.add_scaled(-1.0, center); - self.center = center.clone(); - } - - // candidate or center changed -> update direction nxt_d = y_cand - y_center - if center_index != self.center_index || candidate_index != self.candidate_index { - self.nxt_d = nxt_d.clone(); - } - - self.center_index = center_index; - self.candidate_index = candidate_index; - } -} - -/// Update of the guess value in the current candidate. -enum SubCandidateUpdate { - Unchanged, - New { dist: Real, value: Real }, - Diff { dist: Real, diff: Real }, -} - -/// Update of the center value in the current candidate. -enum SubCenterUpdate { - Unchanged, - New { value: Real }, - Diff { diff: Real }, -} - -/// A subproblem model for guessing candidate and center values. -trait SubProblem { - /// Add a function value at the given evaluation point. - /// - /// ## Parameters - /// - `y`: the evaluation point - /// - `value`: the function value - /// - /// The function returns an information update about the model's (estimated) - /// value in the candidate. - fn new_function_value(&mut self, y: &EvalPoint, value: Real) -> SubCandidateUpdate; - - /// Add a new minorant at the given evaluation point. - /// - /// The minorant is always centered at the current center. - /// - /// ## Parameters - /// - `y`: the evaluation point - /// - `minorant`: the new minorant - /// - /// The function returns the difference in both, the guess of the - /// current candidate AND the current center lower bound. - fn new_minorant(&mut self, y: &EvalPoint, minorant: &Minorant) -> (SubCenterUpdate, SubCandidateUpdate); - - /// Set the new candidate. - /// - /// The function gets the index, the candidate point and the model value. - /// - /// The function returns an initial guess for the new candidate or `None` if - /// there is no initial guess. - fn set_candidate( - &mut self, - index: usize, - y: &Arc, - nxt_d: &Arc, - value: Real, - ) -> SubCandidateUpdate; - - /// Move the new center to the current candidate. - /// - /// The function returns an (optional) initial lower bound for the center value. - fn move_center(&mut self, d: &Arc) -> Option; - - /// Return the current center guess value of this subproblem. - fn cur_guess_value(&self) -> Real; - - /// Return the current center cut value (lower bound) of this subproblem. - fn cur_cut_value(&self) -> Real; - - /// Return the distance measure to the current candidate. - fn eval_distance(&self) -> Real; -} - /// Model data of a single subproblem. /// /// This struct does not handle the subproblem model itself. However, it handles /// the asynchronous precision data, i.e. the guessed Lipschitz-constant and the /// distance of the evaluation points to the candidate. @@ -611,126 +490,134 @@ /// and the center must be provided by an implementation of `SubProblem`. struct SubData { /// The index associated with this subproblem. fidx: usize, /// The subproblem. - sub: Box, - /// The current candidate index. - candidate_index: usize, + sub: Box, + /// The current center. + center: Point, + /// The current candidate. + candidate: Point, /// The last index at which the evaluation has been started. last_eval_index: usize, /// Whether a subproblem evaluation is currently running. is_running: bool, /// Whether the last evaluation has been sufficiently close. is_close_enough: bool, + /// The original guess value and its evaluation distance in the current center. + center_guess: Guess, /// The current guess of the Lipschitz constant. l_guess: Real, } impl SubData { - fn new(fidx: usize, sub: Box) -> SubData { + fn new(fidx: usize, sub: Box, y: &Point) -> SubData { SubData { fidx, sub, last_eval_index: 0, - candidate_index: 0, + center: y.clone(), + candidate: y.clone(), is_running: false, is_close_enough: false, + center_guess: Guess::default(), l_guess: 0.0, } } - fn new_function_value(&mut self, y: &EvalPoint, value: Real, accept_factor: Real) -> SubCandidateUpdate { - let update = self.sub.new_function_value(y, value); - self.update_close_enough(&update, y.index, accept_factor); - - match update { - SubCandidateUpdate::Diff { dist, .. } | SubCandidateUpdate::New { dist, .. } => debug!( - "Improved candidate fidx:{} eval:{} new-dist:{} l:{}", - self.fidx, y.index, dist, self.l_guess - ), - _ => (), - }; - - update - } - - fn new_minorant(&mut self, y: &EvalPoint, minorant: &Minorant) -> (SubCenterUpdate, SubCandidateUpdate) { - self.sub.new_minorant(y, minorant) - } - - fn set_candidate( - &mut self, - index: usize, - y: &Arc, - nxt_d: &Arc, - value: Real, - accept_factor: Real, - ) -> SubCandidateUpdate { - let update = self.sub.set_candidate(index, y, nxt_d, value); - let dist = match update { - SubCandidateUpdate::New { dist, .. } | SubCandidateUpdate::Diff { dist, .. } => dist, - _ => 0.0, - }; - - self.update_close_enough(&update, index, accept_factor); - self.candidate_index = index; - - debug!( - "Old evaluation {}sufficient fidx:{} at:{} dist:{}", - if self.is_close_enough { "" } else { "in " }, - self.fidx, - index, - dist, - ); - - update - } - - /// Move the center to the current candidate. + /// Set the center of this model. /// - /// If `update_l_guess` is `true`, also update the guess of the Lipschitz - /// constant. - fn move_center(&mut self, d: &Arc, update_l_guess: bool) -> Option { - let eval_dist = self.sub.eval_distance(); - let cur_guess_value = self.sub.cur_guess_value(); - let cur_cutvalue = self.sub.move_center(d); + /// If `update_l_guess` is true also update the guess of the Lipschitz constant. + fn move_center(&mut self, y: &Point, update_l_guess: bool) { + assert_eq!(y.index, self.candidate.index, "Must move to current candidate"); + + // The guess value used in the current (i.e. old) center + let old_guess = self.center_guess; + // The cut value now known for the center. + let old_cutvalue = self.sub.get_lower_bound(&self.center); // There has been a previous evaluation, so first update the Lipschitz guess ... - if let Some(cur_cutvalue) = cur_cutvalue { - if update_l_guess && eval_dist > 0.0 { - let new_l_guess = (cur_cutvalue - cur_guess_value) / eval_dist; - if new_l_guess > self.l_guess { - debug!( - "New l_guess fidx:{} old-L:{} L:{}", - self.fidx, self.l_guess, new_l_guess - ); - self.l_guess = new_l_guess; - } - } - } - - cur_cutvalue - } - - fn update_close_enough(&mut self, update: &SubCandidateUpdate, eval_index: usize, accept_factor: Real) { - match update { - SubCandidateUpdate::New { dist, .. } | SubCandidateUpdate::Diff { dist, .. } => { - self.is_close_enough = eval_index == self.candidate_index || dist * self.l_guess <= accept_factor; - } - _ => (), - } - } - - /// Return the current guessed value in the center. - fn cur_guess_value(&self) -> Real { - self.sub.cur_guess_value() + if update_l_guess && old_guess.dist > Real::zero() { + debug!( + "L-guess fidx:{} guess:{} cut:{}", + self.fidx, old_guess.value, old_cutvalue + ); + let new_l_guess = if old_guess.dist.is_finite() { + (old_cutvalue - old_guess.value) / old_guess.dist + } else { + Real::zero() + }; + + if new_l_guess > self.l_guess { + debug!( + "New l_guess fidx:{} old-L:{} L:{}", + self.fidx, self.l_guess, new_l_guess + ); + self.l_guess = new_l_guess; + } + } + + // Save the new center + self.center = y.clone(); + + // Save guess value of the candidate/new center + self.center_guess = self.sub.get_guess_value(&self.center); + } + + /// Set the candidate of this model. + fn update_candidate(&mut self, y: &Point, accept_factor: Real) { + self.candidate = y.clone(); + self.update_close_enough(accept_factor); + } + + /// Add a function value to this model. + /// + /// The `accept_factor` is a parameter for possibly accepting + /// the candidate guess value as "good enough" (if it has been + /// changed by the new minorant). + fn add_function_value(&mut self, y: &Point, value: Real, accept_factor: Real) { + self.sub.add_function_value(y, value); + self.update_close_enough(accept_factor) + } + + /// Add a minorant to this model. + /// + /// The `accept_factor` is a parameter for possibly accepting + /// the candidate guess value as "good enough" (if it has been + /// changed by the new minorant). + /// + /// The minorant must be centered at the global 0. + fn add_minorant(&mut self, y: &Point, m: &Arc, accept_factor: Real) { + self.sub.add_minorant(y, m); + self.update_close_enough(accept_factor) + } + + /// Return the current guess value at the given point. + fn get_guess_value(&mut self, y: &Point) -> Guess { + self.sub.get_guess_value(y) + } + + /// Return the lower bound at the given point. + fn get_lower_bound(&mut self, y: &Point) -> Real { + self.sub.get_lower_bound(y) + } + + fn update_close_enough(&mut self, accept_factor: Real) { + let g = self.sub.get_guess_value(&self.candidate); + self.is_close_enough = g.is_exact() || g.dist * self.l_guess <= accept_factor } /// Return the error estimation in the center. - fn error_estimate(&self) -> Real { - self.sub.cur_cut_value() - self.sub.cur_guess_value() + /// + /// This is the difference between the (current) lower bound and the used + /// guess value. + fn error_estimate(&mut self) -> Real { + self.sub.get_lower_bound(&self.center) - self.center_guess.value + } + + fn center_guess_value(&self) -> Real { + self.center_guess.value } } /// Implementation of a parallel bundle method. pub struct Solver @@ -765,11 +652,11 @@ /// Whether there is currently a master computation running. master_running: bool, /// Whether the master problem has been changed. - master_need_resolve: bool, + master_has_changed: bool, /// The channel to receive the evaluation results from subproblems. client_tx: Option>, /// The channel to receive the evaluation results from subproblems. @@ -810,11 +697,11 @@ data: SolverData::default(), threadpool: ThreadPool::with_name("Parallel bundle solver".to_string(), ncpus), master, master_proc: None, - master_need_resolve: false, + master_has_changed: false, master_running: false, client_tx: None, client_rx: None, @@ -858,11 +745,14 @@ debug!("Initialize solver"); let n = self.problem.num_variables(); let m = self.problem.num_subproblems(); - self.data.init(dvec![0.0; n]); + self.data.init(Point { + index: 1, + point: Arc::new(dvec![Real::zero(); n]), + }); let (tx, rx) = channel(); self.client_tx = Some(tx.clone()); self.client_rx = Some(rx); @@ -897,36 +787,22 @@ tx, &mut self.threadpool, )); debug!("Initial problem evaluation"); + // The initial evaluation point. + self.data.nxt_y = self.data.cur_y.clone(); // We need an initial evaluation of all oracles for the first center. self.data.subs = (0..m) - .map(|fidx| SubData::new(fidx, Box::new(SubZero::new()))) + .map(|fidx| SubData::new(fidx, Box::new(NearestValue::new()), &self.data.cur_y)) .collect(); - self.data.nxt_y = self.data.cur_y.clone(); - - // The initial evaluation point. - let evalpoint = EvalPoint::new( - self.data.candidate_index, - self.data.nxt_y.clone(), - self.data.candidate_index, - &self.data.cur_y, - self.data.candidate_index, - &self.data.nxt_y, - &self.data.nxt_d, - ); // This could be done better: the initial evaluation has index 1! - self.data.candidates.push(evalpoint.clone()); - self.data.candidates.push(evalpoint.clone()); - self.data.candidate_index = 1; + self.data.candidates.push(self.data.nxt_y.clone()); + self.data.candidates.push(self.data.nxt_y.clone()); + for i in 0..m { - match self.data.subs[i].set_candidate(0, &self.data.nxt_y, &self.data.nxt_d, -Real::infinity(), -1.0) { - SubCandidateUpdate::Unchanged => (), - _ => panic!("first candidate must have no guess"), - } self.evaluate_subproblem(i)?; } self.start_time = Instant::now(); @@ -934,43 +810,38 @@ let mut cnt_remaining = self.problem.num_subproblems(); let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; let client_rx = self.client_rx.as_ref().ok_or(Error::NotInitialized)?; - self.data.cur_val = 0.0; - self.data.nxt_val = 0.0; - while cnt_remaining > 0 { let msg = client_rx.recv(); match msg? { Message::Eval(m) => match m { EvalResult::ObjectiveValue { index, value } => { assert_eq!( - index.candidate_index, self.data.candidate_index, + index.candidate_index, self.data.nxt_y.index, "Receive objective value for unexpected candidate" ); - self.data.nxt_val += - match self.data.subs[index.subproblem].new_function_value(&evalpoint, value, -1.0) { - SubCandidateUpdate::Unchanged => 0.0, - SubCandidateUpdate::Diff { diff, .. } => diff, - SubCandidateUpdate::New { value, .. } => value, - }; + self.data.subs[index.subproblem].add_function_value(&self.data.nxt_y, value, 0.0); } EvalResult::Minorant { index, - minorant, + mut minorant, primal, } => { assert_eq!( - index.candidate_index, self.data.candidate_index, + index.candidate_index, self.data.nxt_y.index, "Receive objective value for unexpected candidate" ); - match self.data.subs[index.subproblem].new_minorant(&evalpoint, &minorant) { - (SubCenterUpdate::New { value }, _) => self.data.cur_val += value, - _ => panic!("unexpected minorant update"), - } - master.add_minorant(index.subproblem, minorant, primal)?; + // Add the minorant to the master problem. + // The minorant is centered at the candidate == center, so it does + // not need to be moved. + master.add_minorant(index.subproblem, minorant.clone(), primal)?; + self.master_has_changed = true; + // Center the minorant at 0. + minorant.move_center(-1.0, &self.data.nxt_y.point); + self.data.subs[index.subproblem].add_minorant(&self.data.nxt_y, &Arc::new(minorant), 0.0); } EvalResult::Done { index } => { self.data.subs[index.subproblem].is_running = false; cnt_remaining -= 1; } @@ -986,10 +857,23 @@ unreachable!("Receive master response during initialization"); } } } + // Set the initial values. + // For the center this is the current lower bound (cut value), + // for the candidate it is the model candidate value. + // + // Note that both are the same ... + self.data.cur_val = 0.0; + self.data.nxt_val = 0.0; + for s in &mut self.data.subs { + self.data.cur_val += s.get_lower_bound(&self.data.nxt_y); + self.data.nxt_val += s.get_guess_value(&self.data.nxt_y).value; + } + assert!((self.data.cur_val - self.data.nxt_val).abs() < 1e-6); + // This is the first evaluation. We effectively get // the function value at the current center but // we do not have a model estimate yet. Hence, we do not know // a good guess for the weight. self.data.cur_weight = Real::infinity(); @@ -1000,11 +884,11 @@ debug!("First Step"); debug!(" cur_val={}", self.data.cur_val); self.data.cnt_descent += 1; // compute the initial candidate - self.update_candidate(true)?; + self.update_candidate()?; self.show_info(Step::Descent); debug!("Initialization complete"); @@ -1016,12 +900,18 @@ let num_variables = self.problem.num_variables(); self.data.max_iter = max_iter; self.data.cnt_updates = 0; self.data.nxt_d = Arc::new(dvec![0.0; num_variables]); - self.data.nxt_y = Arc::new(dvec![]); - self.data.nxt_val = 0.0; + self.data.nxt_y = self.data.cur_y.clone(); + let nxt_y = &self.data.nxt_y; + self.data.nxt_val = self + .data + .subs + .iter_mut() + .map(|s| s.get_guess_value(&nxt_y).value) + .sum::(); self.data.need_update = true; self.data.update_in_progress = false; } /// Solve the problem with the default maximal iteration limit [`DEFAULT_ITERATION_LIMIT`]. @@ -1059,19 +949,19 @@ let msg = self.client_rx.as_ref().ok_or(Error::NotInitialized)?.recv()?; match msg { Message::Eval(m) => { // Receive a evaluation result if self.handle_client_response(m)? { - return Ok(false); + return Ok(true); } } Message::Update(msg) => { debug!("Receive update response"); if self.handle_update_response(msg)? { // The master problem has been changed so we need a new // candidate as well. - self.update_candidate(true)?; + self.update_candidate()?; } } Message::Master(mresponse) => { debug!("Receive master response"); // Receive result (new candidate) from the master @@ -1127,47 +1017,39 @@ /// A new objective value has been computed. fn handle_new_objective(&mut self, id: EvalId, value: Real) -> Result { debug!( "Receive objective from subproblem:{} candidate:{} current:{} obj:{}", - id.subproblem, id.candidate_index, self.data.candidate_index, value + id.subproblem, id.candidate_index, self.data.nxt_y.index, value ); // check if new evaluation is closer to the current candidate let sub = &mut self.data.subs[id.subproblem]; // Whether the subproblem evaluation had been good enough. let was_close_enough = sub.is_close_enough; - if let Some(evalpoint) = self.data.candidates.get_mut(id.candidate_index) { + if let Some(nxt) = self.data.candidates.get_mut(id.candidate_index) { let accept_factor = self.params.imprecision_factor * self.params.acceptance_factor * self.data.expected_progress / self.problem.num_subproblems().to_f64().unwrap(); - // possibly update internal data of the EvalPoint - evalpoint.update( - self.data.center_index, - &self.data.cur_y, - self.data.candidate_index, - &self.data.nxt_y, - &self.data.nxt_d, - ); - - match sub.new_function_value(evalpoint, value, accept_factor) { - // candidate is not closer -> ignore - SubCandidateUpdate::Unchanged => return Ok(false), - SubCandidateUpdate::New { value, .. } => self.data.nxt_val += value, - SubCandidateUpdate::Diff { diff, .. } => self.data.nxt_val += diff, - }; + let old_can_val = sub.get_guess_value(&self.data.nxt_y).value; + let old_cen_val = sub.get_lower_bound(&self.data.cur_y); + sub.add_function_value(nxt, value, accept_factor); + let new_can_val = sub.get_guess_value(&self.data.nxt_y).value; + let new_cen_val = sub.get_lower_bound(&self.data.cur_y); + + self.data.nxt_val += new_can_val - old_can_val; + self.data.cur_val += new_cen_val - old_cen_val; // This is just a safe-guard: if the function has been evaluated at // the current candidate, the evaluation *must* be good enough. assert!( - sub.is_close_enough || sub.last_eval_index < self.data.candidate_index, - "Unexpected insufficiency fidx:{} l:{} dist:{}", + sub.is_close_enough || sub.last_eval_index < self.data.nxt_y.index, + "Unexpected insufficiency fidx:{} l:{}", id.subproblem, sub.l_guess, - evalpoint.candidate_dist ); } else { // unknown candidate -> ignore objective value warn!("Ignore unknown candidate index:{}", id.candidate_index); return Ok(false); @@ -1174,140 +1056,67 @@ } // Test if the new candidate is close enough for the asynchronous // precision test. if !was_close_enough && sub.is_close_enough { - self.data.num_insufficient_candidates -= 1; - debug!( - "Accept result fidx:{} index:{} candidate:{} (remaining insufficient: {})", - id.subproblem, id.candidate_index, self.data.candidate_index, self.data.num_insufficient_candidates - ); - } - - // If not all subproblems have reached sufficient precision, stop - // (eventually all subproblems will be evaluated at the center). - if self.data.num_insufficient_candidates > 0 { - return Ok(false); - } - - self.do_step() - } - - fn do_step(&mut self) -> Result { - let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; - let descent_bnd = Self::get_descent_bound(self.params.acceptance_factor, &self.data); - - // Test whether we do a descent step - if self.data.nxt_val <= descent_bnd { - debug!("Descent Step"); - debug!(" cur_val ={}", self.data.cur_val); - debug!(" nxt_mod ={}", self.data.nxt_mod); - debug!(" nxt_ub ={}", self.data.nxt_val); - debug!(" descent_bnd={}", descent_bnd); - - self.data.cnt_descent += 1; - self.data.center_index = self.data.candidate_index; - - // Note that we must update the weight *before* we - // change the internal data, so the old information - // that caused the descent step is still available. - self.data.cur_weight = self.weighter.descent_weight(&self.data); - self.data.cur_y = self.data.nxt_y.clone(); - // The new value in the center is the model value in the candidate. - // In particular, it is a lower bound on the real function value. - // - // Note that we do not use the model value `nxt_mod`, but the - // sum of the single model values, because the latter might be higher - // in case of an aggregated model. - self.data.cur_val = self.data.nxt_submods.iter().sum(); - self.data.nxt_val = Real::infinity(); - - // Check if the progress of the last decent step was large enough - // when using the lower bound in the center instead of the former - // guess value. - let error = self.data.subs.iter().map(SubData::error_estimate).sum::(); - let update_l_guess = error > self.data.error_bound; - - // save new error bound - self.data.error_bound = self.data.expected_progress * self.params.acceptance_factor; - - // Move all subproblems. - for sub in &mut self.data.subs { - sub.move_center(&self.data.nxt_d, update_l_guess); - } - self.data.need_update = true; - - master.move_center(1.0, self.data.nxt_d.clone(), self.data.center_index)?; - master.set_weight(self.data.cur_weight)?; - - self.show_info(Step::Descent); - self.update_problem(Step::Descent)?; - - // We need a new candidate. - self.update_candidate(true)?; - Ok(self.data.cnt_descent >= self.data.max_iter) - } else { - // No descent-step, so this is declared a null step - self.data.cur_weight = self.weighter.null_weight(&self.data); - self.show_info(Step::Null); - self.update_problem(Step::Null)?; - - // After a null step we need a new candidate, too. However, in this - // case any new candidate will do, so we only start a new master - // problem evaluation if there is no running computation. - // - // TODO: does this make sense? - self.update_candidate(false)?; - Ok(false) - } - } - - /// Add a new minorant. - fn handle_new_minorant(&mut self, id: EvalId, minorant: Minorant, primal: P::Primal) -> Result { - debug!( - "Receive minorant subproblem:{} candidate:{} current:{} center:{}", - id.subproblem, id.candidate_index, self.data.candidate_index, self.data.center_index, - ); - - let sub = &mut self.data.subs[id.subproblem]; - let mut minorant = minorant; - if let Some(evalpoint) = self.data.candidates.get_mut(id.candidate_index) { - // possibly update internal data of the EvalPoint - evalpoint.update( - self.data.center_index, - &self.data.cur_y, - self.data.candidate_index, - &self.data.nxt_y, - &self.data.nxt_d, - ); - - // move center of minorant to cur_y - minorant.move_center(-1.0, &evalpoint.center_d); - // add minorant to submodel - let (cur_cutvalue_diff, nxt_guess_diff) = sub.new_minorant(&evalpoint, &minorant); - - self.data.nxt_val += match nxt_guess_diff { - SubCandidateUpdate::Diff { diff, .. } => diff, - SubCandidateUpdate::New { value, .. } => value, - SubCandidateUpdate::Unchanged => 0.0, - }; - - self.data.cur_val += match cur_cutvalue_diff { - SubCenterUpdate::Diff { diff } => diff, - SubCenterUpdate::New { value } => value, - SubCenterUpdate::Unchanged => 0.0, - }; + debug!( + "Accept result fidx:{} index:{} candidate:{} (remaining insufficient: {})", + id.subproblem, id.candidate_index, self.data.nxt_y.index, self.data.num_insufficient_candidates + ); + } + + self.maybe_do_step(false) + } + + /// Add a new minorant. + /// + /// The minorant is added to the submodel as well as the master problem. The modified submodel + /// may then lead to a potential null/descent-step. + /// + /// Return values + /// - `Ok(true)` if the termination criterion has been satisfied, + /// - `Ok(false)` if the termination criterion has not been satisfied, + /// - `Err(_)` on error. + fn handle_new_minorant(&mut self, id: EvalId, minorant: Minorant, primal: P::Primal) -> Result { + debug!( + "Receive minorant subproblem:{} candidate:{} current:{} center:{}", + id.subproblem, id.candidate_index, self.data.nxt_y.index, self.data.cur_y.index, + ); + + let accept_factor = + self.params.imprecision_factor * self.params.acceptance_factor * self.data.expected_progress + / self.problem.num_subproblems().to_f64().unwrap(); + + let sub = &mut self.data.subs[id.subproblem]; + let mut minorant = minorant; + if let Some(point) = self.data.candidates.get_mut(id.candidate_index) { + // center the minorant at 0 + minorant.move_center(-1.0, &point.point); + + // add minorant to submodel + let old_can_val = sub.get_guess_value(&self.data.nxt_y).value; + let old_cen_val = sub.get_lower_bound(&self.data.cur_y); + sub.add_minorant(&point, &Arc::new(minorant.clone()), accept_factor); + let new_can_val = sub.get_guess_value(&self.data.nxt_y).value; + let new_cen_val = sub.get_lower_bound(&self.data.cur_y); + + self.data.nxt_val += new_can_val - old_can_val; + self.data.cur_val += new_cen_val - old_cen_val; + + // center the minorant at the current center + minorant.move_center(1.0, &self.data.cur_y.point); } else { warn!("Ignore unknown candidate index:{}", id.candidate_index); return Ok(false); } // add minorant to master problem let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; master.add_minorant(id.subproblem, minorant, primal)?; + self.master_has_changed = true; - Ok(false) + self.maybe_do_step(false) } /// Handle a response `master_res` from the master problem process. /// /// The master response is the new candidate point. The method updates the @@ -1316,10 +1125,12 @@ /// a termination criterion is satisfied. This is only the case if there is /// no pending problem update. /// /// Finally the master problem starts the evaluation of all subproblems at /// the new candidate. + /// + /// The new candidate is immediately checked for a potential new test. /// /// Return values /// - `Ok(true)` if the termination criterion has been satisfied, /// - `Ok(false)` if the termination criterion has not been satisfied, /// - `Err(_)` on error. @@ -1341,11 +1152,11 @@ self.data.nxt_submods.clear(); self.data.nxt_submods.extend(nxt_submods); debug!( "Master Response current_center:{} current_candidate:{} res_center:{} nxt_mod:{}", - self.data.center_index, self.data.candidate_index, center_index, self.data.nxt_mod + self.data.cur_y.index, self.data.nxt_y.index, center_index, self.data.nxt_mod ); } }; self.data.expected_progress = self.data.cur_val - self.data.nxt_mod; @@ -1357,145 +1168,268 @@ // we use its result as to make a good guess for the initial weight // of the proximal term and resolve. if self.data.cur_weight.is_infinite() { self.data.cur_weight = self.weighter.initial_weight(&self.data); master.set_weight(self.data.cur_weight)?; - self.update_candidate(true)?; + self.master_has_changed = true; + self.update_candidate()?; return Ok(false); } - if self.terminator.terminate(&self.data) - && !self.data.update_in_progress - && self.data.cnt_descent > 2 - && !self.data.need_update - { - self.show_info(Step::Term); - info!("Termination criterion satisfied"); - return Ok(true); - } - // Compress bundle master.compress()?; + + // Check if new variables had been added. In this case, resize cur_y. + if self.data.nxt_d.len() > self.data.cur_y.point.len() { + let nnew = self.data.nxt_d.len() - self.data.cur_y.point.len(); + if nnew != self.data.cur_y.point.len() { + let mut cur_y = self.data.cur_y.point.as_ref().clone(); + cur_y.extend(repeat(0.0).take(nnew)); + self.data.cur_y.point = Arc::new(cur_y); + } + } // Compute new candidate. let mut next_y = dvec![]; - self.data.candidate_index += 1; - - // Check if new variables had been added. In this case, resize cur_y. - if self.data.nxt_d.len() > self.data.cur_y.len() { - let nnew = self.data.nxt_d.len() - self.data.cur_y.len(); - if nnew != self.data.cur_y.len() { - let mut cur_y = self.data.cur_y.as_ref().clone(); - cur_y.extend(repeat(0.0).take(nnew)); - self.data.cur_y = Arc::new(cur_y); - } - } - - next_y.add(&self.data.cur_y, &self.data.nxt_d); + next_y.add(&self.data.cur_y.point, &self.data.nxt_d); #[cfg(debug_assertions)] { - if self.data.nxt_y.len() == next_y.len() { - let mut diff = self.data.nxt_y.as_ref().clone(); + if self.data.nxt_y.point.len() == next_y.len() { + let mut diff = self.data.nxt_y.point.as_ref().clone(); diff.add_scaled(-1.0, &next_y); debug!("Candidates move distance:{}", diff.norm2()); } } - self.data.nxt_y = Arc::new(next_y); + self.data.nxt_y.point = Arc::new(next_y); + self.data.nxt_y.index += 1; - // Reset evaluation data. - self.data.nxt_val = 0.0; + // Add the new candidate to the list of candidates. + debug_assert_eq!(self.data.nxt_y.index, self.data.candidates.len()); + self.data.candidates.push(self.data.nxt_y.clone()); - // Compute a new guess for the function value at the new candidate. + // Update the candidate in all submodels and + // compute first guess value for new candidate. let accept_factor = self.params.imprecision_factor * self.params.acceptance_factor * self.data.expected_progress / self.problem.num_subproblems().to_f64().unwrap(); - - self.data.num_insufficient_candidates = 0; - - // Create a new evaluation point with the current center and (new) candidate. - let candidate_index = self.data.candidates.len(); - self.data.candidates.push(EvalPoint::new( - candidate_index, - self.data.nxt_y.clone(), - self.data.center_index, - &self.data.cur_y, - self.data.candidate_index, - &self.data.nxt_y, - &self.data.nxt_d, - )); - - for (fidx, sub) in self.data.subs.iter_mut().enumerate() { - self.data.nxt_val += match sub.set_candidate( - self.data.candidate_index, - &self.data.nxt_y, - &self.data.nxt_d, - self.data.nxt_submods[fidx], - accept_factor, - ) { - SubCandidateUpdate::Unchanged => 0.0, - SubCandidateUpdate::New { value, .. } => value, - SubCandidateUpdate::Diff { .. } => todo!("Only `New` is supported currently"), - }; - - if !sub.is_close_enough { - self.data.num_insufficient_candidates += 1; - } - } - debug!( - "Number of insufficient subproblems: {}", - self.data.num_insufficient_candidates - ); + self.data.nxt_val = Real::zero(); + for sub in self.data.subs.iter_mut() { + sub.update_candidate(&self.data.nxt_y, accept_factor); + self.data.nxt_val += sub.get_guess_value(&self.data.nxt_y).value; + } // Start evaluation of all subproblems at the new candidate. for i in 0..self.data.subs.len() { self.evaluate_subproblem(i)?; } + + self.maybe_do_step(true) + } + + /// Do a descent or null step if precision is sufficient. + /// + /// Also checks if the termination criterion is satisfied. + /// + /// Return values + /// - `Ok(true)` if the termination criterion has been satisfied, + /// - `Ok(false)` if the termination criterion has not been satisfied, + /// - `Err(_)` on error. + fn maybe_do_step(&mut self, check_termination: bool) -> Result { + // No step if there is no real new candidate + if self.data.nxt_y.index == self.data.cur_y.index { + return Ok(false); + } + + self.data.num_insufficient_candidates = self.data.subs.iter().filter(|s| !s.is_close_enough).count(); + + let nxt_y = &self.data.nxt_y; + let num_exact = self + .data + .subs + .iter_mut() + .map(|s| s.get_guess_value(nxt_y).is_exact()) + .filter(|&is_exact| is_exact) + .count(); + + let sum_dist = self + .data + .subs + .iter_mut() + .map(|s| s.get_guess_value(nxt_y).dist) + .sum::(); + + debug!( + "Number of insufficient subproblems: {} num exact: {} sum dist: {}", + self.data.num_insufficient_candidates, num_exact, sum_dist, + ); + + // test if we should terminate + if check_termination + && self.terminator.terminate(&self.data) + && !self.data.update_in_progress + && self.data.cnt_descent > 2 + && !self.data.need_update + { + if self.data.num_inexact_center.is_zero() { + self.show_info(Step::Term); + info!("Termination criterion satisfied"); + return Ok(true); + } + + // The termination criterion has been satisfied, but the current + // center evaluations are not exact. We force the current + // *candidate* evaluations to be exact and do a forced descent step. + // This causes the next center to be exact and hopefully satisfying + // the termination criterion, too. + let num_inexact = self.data.subs.len() - num_exact; + + debug!( + "Termination criterion satisfied with {} inexact evaluations", + num_inexact + ); + + // Current candidate is exact, force a descent step. + self.data.force_descent = true; + if num_inexact.is_zero() { + return self.do_step(); + } + + // Otherwise just wait for the exact evaluations. + return Ok(false); + } + + if check_termination { + self.data.force_descent = false; + } if self.data.num_insufficient_candidates == 0 { self.do_step() } else { Ok(false) } } + + /// Do a null-step or descent-step based on the current candidate data. + fn do_step(&mut self) -> Result { + let master = self.master_proc.as_mut().ok_or(Error::NotInitialized)?; + let descent_bnd = Self::get_descent_bound(self.params.acceptance_factor, &self.data); + + debug!( + "Try step from center:{} to candidate:{}", + self.data.cur_y.index, self.data.nxt_y.index + ); + + // Test whether we do a descent step + if self.data.nxt_val <= descent_bnd || self.data.force_descent { + debug!("Descent Step{}", if self.data.force_descent { " (forced)" } else { "" }); + debug!(" cur_val ={}", self.data.cur_val); + debug!(" nxt_mod ={}", self.data.nxt_mod); + debug!(" nxt_ub ={}", self.data.nxt_val); + debug!(" descent_bnd={}", descent_bnd); + + self.data.force_descent = false; + self.data.cnt_descent += 1; + self.data.cur_y = self.data.nxt_y.clone(); + let cur_y = &self.data.cur_y; + self.data.num_inexact_center = self + .data + .subs + .iter_mut() + .map(|s| !s.get_guess_value(&cur_y).is_exact()) + .filter(|&is_exact| is_exact) + .count(); + + // Note that we must update the weight *before* we + // change the internal data, so the old information + // that caused the descent step is still available. + self.data.cur_weight = self.weighter.descent_weight(&self.data); + // The new value in the center is the model value in the candidate. + // In particular, it is a lower bound on the real function value. + // + // Note that we do not use the model value `nxt_mod`, but the + // sum of the single model values, because the latter might be higher + // in case of an aggregated model. + self.data.cur_val = self.data.nxt_submods.iter().sum(); + + // Check if the progress of the last decent step was large enough + // when using the lower bound in the center instead of the former + // guess value. + let error = self.data.subs.iter_mut().map(SubData::error_estimate).sum::(); + let update_l_guess = error > self.data.error_bound; + + // save new error bound + self.data.error_bound = self.data.expected_progress * self.params.acceptance_factor; + + // Move all subproblems. + self.data.nxt_val = Real::zero(); + for sub in &mut self.data.subs { + sub.move_center(&self.data.cur_y, update_l_guess); + self.data.nxt_val += sub.get_guess_value(&self.data.nxt_y).value; + } + self.data.need_update = true; + + master.move_center(1.0, self.data.nxt_d.clone(), self.data.cur_y.index)?; + master.set_weight(self.data.cur_weight)?; + self.master_has_changed = true; + + self.show_info(Step::Descent); + self.update_problem(Step::Descent)?; + + // We need a new candidate. + self.update_candidate()?; + Ok(self.data.cnt_descent >= self.data.max_iter) + } else { + debug!("Null Step nxt_val:{} descent_bnd:{}", self.data.nxt_val, descent_bnd); + // No descent-step, so this is declared a null step + self.data.cur_weight = self.weighter.null_weight(&self.data); + self.show_info(Step::Null); + self.update_problem(Step::Null)?; + + // After a null step we need a new candidate, too. + self.update_candidate()?; + Ok(false) + } + } /// Start evaluation of a subproblem if it is not running. /// /// The evaluation is started at the current candidate. The candidate /// is added to the subproblem's candidate list. /// /// Returns `true` iff a new evaluations has been started. fn evaluate_subproblem(&mut self, subproblem: usize) -> Result { let sub = &mut self.data.subs[subproblem]; - if !sub.is_running && sub.last_eval_index < self.data.candidate_index { + if !sub.is_running && sub.last_eval_index < self.data.nxt_y.index { sub.is_running = true; - sub.last_eval_index = sub.last_eval_index.max(self.data.candidate_index); + sub.last_eval_index = sub.last_eval_index.max(self.data.nxt_y.index); self.problem .evaluate( subproblem, - self.data.nxt_y.clone(), + self.data.nxt_y.point.clone(), ChannelResultSender::new( EvalId { subproblem, - candidate_index: self.data.candidate_index, + candidate_index: self.data.nxt_y.index, }, self.client_tx.as_ref().ok_or(Error::NotInitialized)?.clone(), ), ) .map_err(Error::Evaluation)?; - debug!("Evaluate fidx:{} candidate:{}", subproblem, self.data.candidate_index); + debug!("Evaluate fidx:{} candidate:{}", subproblem, self.data.nxt_y.index); Ok(true) } else { assert!(sub.is_running || sub.is_close_enough); Ok(false) } } /// Possibly start a new master process computation. - fn update_candidate(&mut self, need_resolve: bool) -> Result<(), P, M> { - self.master_need_resolve = self.master_need_resolve || need_resolve; - if !self.master_running && self.master_need_resolve { + fn update_candidate(&mut self) -> Result<(), P, M> { + if !self.master_running && self.master_has_changed { + debug!("Start master problem"); self.master_running = true; + self.master_has_changed = false; self.master_proc .as_mut() .ok_or(Error::NotInitialized)? .solve(self.data.cur_val)?; } @@ -1519,17 +1453,20 @@ return Ok(false); } // Ok, we are doing a new update now ... self.data.need_update = false; + + // TODO: fix this + return Ok(false); let master_proc = self.master_proc.as_mut().unwrap(); self.problem .update( UpdateData { - cur_y: self.data.cur_y.clone(), - nxt_y: self.data.nxt_y.clone(), + cur_y: self.data.cur_y.point.clone(), + nxt_y: self.data.nxt_y.point.clone(), step, primal_aggrs: (0..self.problem.num_subproblems()) .map(|i| { master_proc .get_aggregated_primal(i) @@ -1539,11 +1476,11 @@ .collect(), }, ChannelUpdateSender::new( EvalId { subproblem: 0, - candidate_index: self.data.center_index, + candidate_index: self.data.cur_y.index, }, self.client_tx.clone().ok_or(Error::NotInitialized)?, ), ) .map_err(Error::Update)?; @@ -1611,11 +1548,11 @@ PrettyPrintFloat(self.data.expected_progress()), PrettyPrintFloat(self.data.cur_val - self.data.nxt_val), PrettyPrintFloat(self.data.nxt_mod), PrettyPrintFloat(self.data.nxt_val), PrettyPrintFloat(self.data.cur_val), - PrettyPrintFloat(self.data.subs.iter().map(|s| s.cur_guess_value()).sum::()), + PrettyPrintFloat(self.data.subs.iter().map(|s| s.center_guess_value()).sum::()), ); } /// Return the aggregated primal of the given subproblem. pub fn aggregated_primal(&self, subproblem: usize) -> Result { ADDED src/solver/asyn/guessmodels.rs Index: src/solver/asyn/guessmodels.rs ================================================================== --- /dev/null +++ src/solver/asyn/guessmodels.rs @@ -0,0 +1,93 @@ +/* + * Copyright (c) 2020 Frank Fischer + * + * 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 + */ + +use super::Point; +use crate::{data::Minorant, Real}; + +use num_traits::{Float, Zero}; +use std::sync::Arc; + +mod nearestvalue; +pub use nearestvalue::NearestValue; + +/// A guessed function value. +#[derive(Clone, Copy)] +pub struct Guess { + /// The guessed value. + pub value: Real, + /// The accuracy distance. + pub dist: Real, +} + +impl Guess { + /// Create a new guess value. + /// + /// If `dist` is zero the value is assumed to be exact. + pub fn new(value: Real, dist: Real) -> Guess { + Guess { value, dist } + } + + /// Create an approximate guess value. + pub fn new_approx(value: Real, dist: Real) -> Guess { + Guess { value, dist } + } + + /// Create an exact guess value. + pub fn new_exact(value: Real) -> Guess { + Guess { + value, + dist: Real::zero(), + } + } + + /// Return `true` if this is an exact guess value. + /// + /// In other words, the value is not `guessed` anymore. + pub fn is_exact(&self) -> bool { + self.dist.is_zero() + } +} + +impl Default for Guess { + fn default() -> Guess { + Guess { + value: Real::zero(), + dist: Real::infinity(), + } + } +} + +/// A subproblem model for guessing candidate and center values. +pub trait GuessModel { + /// Add a function value to this model. + fn add_function_value(&mut self, y: &Point, value: Real); + + /// Add a minorant to this model. + /// + /// The minorant must be centered at the global 0. + fn add_minorant(&mut self, y: &Point, m: &Arc); + + /// Return a guess value at the given point. + /// + /// A guess value is an approximation of the function value in + /// the given point. If the returned guess value has distance + /// zero, it must be exact. + fn get_guess_value(&mut self, y: &Point) -> Guess; + + /// Return a lower bound at the given point. + fn get_lower_bound(&mut self, y: &Point) -> Real; +} ADDED src/solver/asyn/guessmodels/nearestvalue.rs Index: src/solver/asyn/guessmodels/nearestvalue.rs ================================================================== --- /dev/null +++ src/solver/asyn/guessmodels/nearestvalue.rs @@ -0,0 +1,135 @@ +/* + * Copyright (c) 2020 Frank Fischer + * + * 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 + */ + +//! Asynchronous subproblem with zero-order approximation. + +use num_traits::{Float, Zero}; +use std::collections::VecDeque; +use std::sync::Arc; + +use crate::data::Minorant; +use crate::Real; + +use super::{Guess, GuessModel, Point}; + +/// The maximal number of last evaluation points to be kept in the model. +#[allow(non_upper_case_globals)] +const MaxEvalPoints: usize = 5; +/// The maximal number of last minorants to be kept in the model. +#[allow(non_upper_case_globals)] +const MaxMinorants: usize = 5; + +/// Information associated with one subproblem. +pub struct NearestValue { + /// The last evaluation points (y, function_value). + eval_points: VecDeque<(Point, Real)>, + + /// The last minorants. + minorants: VecDeque>, + + /// The last computed guess value (point, guess). + last_guess: Option<(Point, Guess)>, + + /// The last computed lower bound (point, value). + last_lower_bound: Option<(Point, Real)>, +} + +impl NearestValue { + pub fn new() -> NearestValue { + NearestValue { + eval_points: VecDeque::new(), + minorants: VecDeque::new(), + last_guess: None, + last_lower_bound: None, + } + } +} + +impl Default for NearestValue { + fn default() -> NearestValue { + NearestValue::new() + } +} + +impl GuessModel for NearestValue { + fn add_function_value(&mut self, y: &Point, value: Real) { + // update value at last guess point + if let Some(ref mut g) = self.last_guess { + let dist = y.distance(&g.0); + if dist.is_zero() || dist < g.1.dist { + g.1 = Guess::new(value, dist); + } + } + + // Add evaluation point to list. + self.eval_points.push_back((y.clone(), value)); + while self.eval_points.len() > MaxEvalPoints { + self.eval_points.pop_front(); + } + } + + fn add_minorant(&mut self, _y: &Point, m: &Arc) { + // update value at last lower-bound point + if let Some(ref mut lb) = self.last_lower_bound { + lb.1 = lb.1.max(m.eval(&lb.0.point)); + } + + // Add minorant to list. + self.minorants.push_back(m.clone()); + while self.minorants.len() > MaxMinorants { + self.minorants.pop_front(); + } + } + + fn get_guess_value(&mut self, y: &Point) -> Guess { + if let Some(ref g) = self.last_guess { + if g.0.index == y.index { + return g.1; + } + } + + let g = self + .eval_points + .iter() + .map(|(x, value)| Guess::new(*value, x.distance(y))) + .min_by(|a, b| a.dist.partial_cmp(&b.dist).unwrap()) + .unwrap_or_else(Guess::default); + self.last_guess = Some((y.clone(), g)); + g + } + + fn get_lower_bound(&mut self, y: &Point) -> Real { + // check if y equals the last evaluation point + if let Some(ref lb) = self.last_lower_bound { + if lb.0.index == y.index { + // ... yes, so just return that value + return lb.1; + } + } + + // full computation of lower bound at y + let lb = self + .minorants + .iter() + .map(|m| m.eval(&y.point)) + .max_by(|a, b| a.partial_cmp(b).unwrap()) + .unwrap_or_else(|| -Real::infinity()); + // save last evaluation + self.last_lower_bound = Some((y.clone(), lb)); + lb + } +} DELETED src/solver/asyn/subzero.rs Index: src/solver/asyn/subzero.rs ================================================================== --- src/solver/asyn/subzero.rs +++ /dev/null @@ -1,151 +0,0 @@ -/* - * Copyright (c) 2020 Frank Fischer - * - * 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 - */ - -//! Asynchronous subproblem with zero-order approximation. - -use num_traits::Float; -use std::sync::Arc; - -use crate::data::Minorant; -use crate::{DVector, Real}; - -use super::{EvalPoint, SubCandidateUpdate, SubCenterUpdate, SubProblem}; - -/// Information associated with one subproblem. -pub struct SubZero { - /// The current candidate point (index, point, model-value). - candidate: Option<(usize, Arc, Real)>, - - /// Closest evaluation point (index, point, value, distance) - closest_eval: Option<(usize, Arc, Real, Real)>, - - /// The current best possible cut value (lower bound) in the center. - cur_cutvalue: Option, - - /// The guess value that had been used in the current center. - cur_guess_value: Real, -} - -impl SubZero { - pub fn new() -> SubZero { - SubZero { - candidate: None, - closest_eval: None, - cur_cutvalue: None, - cur_guess_value: Real::infinity(), - } - } -} - -impl SubProblem for SubZero { - fn new_function_value(&mut self, y: &EvalPoint, value: Real) -> SubCandidateUpdate { - if self.candidate.is_some() { - if let Some(ref mut eval) = self.closest_eval { - if y.candidate_dist >= eval.3 { - // New evaluation is not closer - return SubCandidateUpdate::Unchanged; - } - let old_value = eval.2; - self.closest_eval = Some((y.index, y.point.clone(), value, y.candidate_dist)); - return SubCandidateUpdate::Diff { - dist: y.candidate_dist, - diff: value - old_value, - }; - } else { - self.closest_eval = Some((y.index, y.point.clone(), value, y.candidate_dist)); - return SubCandidateUpdate::New { - dist: y.candidate_dist, - value, - }; - } - } - SubCandidateUpdate::Unchanged - } - - fn new_minorant(&mut self, _y: &EvalPoint, minorant: &Minorant) -> (SubCenterUpdate, SubCandidateUpdate) { - if let Some(cur_cutvalue) = self.cur_cutvalue { - if cur_cutvalue < minorant.constant { - self.cur_cutvalue = Some(minorant.constant); - ( - SubCenterUpdate::Diff { - diff: minorant.constant - cur_cutvalue, - }, - SubCandidateUpdate::Unchanged, - ) - } else { - (SubCenterUpdate::Unchanged, SubCandidateUpdate::Unchanged) - } - } else { - self.cur_cutvalue = Some(minorant.constant); - ( - SubCenterUpdate::New { - value: minorant.constant, - }, - SubCandidateUpdate::Unchanged, - ) - } - } - - fn set_candidate( - &mut self, - index: usize, - y: &Arc, - _nxt_d: &Arc, - value: Real, - ) -> SubCandidateUpdate { - if let Some(ref mut eval) = self.closest_eval { - // we simply reuse the latest evaluation point -> update distance - let mut d = eval.1.as_ref().clone(); - d.add_scaled(-1.0, &y); - let dist = d.norm2(); - self.candidate = Some((index, y.clone(), value)); - eval.3 = dist; - SubCandidateUpdate::New { dist, value: eval.2 } - } else { - assert_eq!(index, 0); - self.candidate = Some((index, y.clone(), value)); - // there is no previous evaluation point, so there is no guess - SubCandidateUpdate::Unchanged - } - } - - fn move_center(&mut self, _d: &Arc) -> Option { - let cand = self.candidate.as_ref().expect("No candidate available"); - if let Some(ref mut eval) = self.closest_eval { - // ... and store the new guess value we used ... - self.cur_guess_value = eval.2; - // ... and use the model value as first lower bound in the center - self.cur_cutvalue = Some(cand.2); - self.cur_cutvalue - } else { - // This should never happen. - unreachable!() - } - } - - fn cur_guess_value(&self) -> Real { - self.cur_guess_value - } - - fn cur_cut_value(&self) -> Real { - self.cur_cutvalue.unwrap_or(-Real::infinity()) - } - - fn eval_distance(&self) -> Real { - self.closest_eval.as_ref().map(|e| e.3).unwrap_or(Real::infinity()) - } -} Index: src/solver/sync.rs ================================================================== --- src/solver/sync.rs +++ src/solver/sync.rs @@ -1,7 +1,7 @@ /* - * Copyright (c) 2019 Frank Fischer + * Copyright (c) 2019, 2020 Frank Fischer * * 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. @@ -21,11 +21,10 @@ use rs_crossbeam::channel::{unbounded as channel, RecvError}; #[cfg(not(feature = "crossbeam"))] use std::sync::mpsc::{channel, RecvError}; use log::{debug, info, warn}; -use num_cpus; use num_traits::Float; use std::sync::Arc; use std::time::Instant; use threadpool::ThreadPool;