// Copyright (c) 2016, 2018 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
// 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 <http://www.gnu.org/licenses/>
//
//! Weight updating rule according to Helmberg and Kiwiel.
//!
//! The procedure is described in
//!
//! > Helmberg, C. and Kiwiel, K.C. (2002): A spectral bundle method
//! > with bounds, Math. Programming A 93, 173--194
//!
use Real;
use {BundleState, SolverParams, Step, Weighter};
use std::cmp::{max, min};
use std::f64::NEG_INFINITY;
const FACTOR: Real = 2.0;
/**
* Weight updating rule according to Helmberg and Kiwiel.
*
* The procedure is described in
*
* > Helmberg, C. and Kiwiel, K.C. (2002): A spectral bundle method
* > with bounds, Math. Programming A 93, 173--194
*/
pub struct HKWeighter {
eps_weight: Real,
m_r: Real,
iter: isize,
model_max: Real,
}
impl HKWeighter {
/// Create a new HKWeighter with default weight $m_R = 0.5$.
pub fn new() -> HKWeighter {
HKWeighter::new_weight(0.5)
}
/// Create new HKWeighter with weight $m_R$.
pub fn new_weight(m_r: Real) -> HKWeighter {
assert!(m_r > 0.0);
HKWeighter {
eps_weight: 1e30,
m_r,
iter: 0,
model_max: NEG_INFINITY,
}
}
}
impl Default for HKWeighter {
fn default() -> Self {
Self::new()
}
}
impl Weighter for HKWeighter {
fn weight(&mut self, state: &BundleState, params: &SolverParams) -> Real {
assert!(params.min_weight > 0.0);
assert!(params.max_weight >= params.min_weight);
debug!("HKWeighter {:?} iter:{}", state.step, self.iter);
if state.step == Step::Term {
self.eps_weight = 1e30;
self.iter = 0;
return if state.cur_y.len() == 0 || state.sgnorm < state.cur_y.len() as Real * 1e-10 {
1.0
} else {
state.sgnorm.max(1e-4)
}.max(params.min_weight)
.min(params.max_weight);
}
let cur_nxt = state.cur_val - state.nxt_val;
let cur_mod = state.cur_val - state.nxt_mod;
let w = 2.0 * state.weight * (1.0 - cur_nxt / cur_mod);
debug!(" cur_nxt={} cur_mod={} w={}", cur_nxt, cur_mod, w);
if state.step == Step::Null {
let sgnorm = state.sgnorm;
let lin_err = state.cur_val - state.new_cutval;
self.eps_weight = self.eps_weight.min(sgnorm + cur_mod - sgnorm * sgnorm / state.weight);
let new_weight = if self.iter < -3 && lin_err > self.eps_weight.max(FACTOR * cur_mod) {
w
} else {
state.weight
}.min(FACTOR * state.weight)
.min(params.max_weight);
if new_weight > state.weight {
self.iter = -1
} else {
self.iter = min(self.iter - 1, -1);
}
debug!(
" sgnorm={} cur_val={} new_cutval={} lin_err={} eps_weight={}",
sgnorm, state.cur_val, state.new_cutval, lin_err, self.eps_weight
);
debug!(" new_weight={}", new_weight);
new_weight
} else {
self.model_max = self.model_max.max(state.nxt_mod);
let new_weight = if self.iter > 0 && cur_nxt > self.m_r * cur_mod {
w
} else if self.iter > 3 || state.nxt_val < self.model_max {
state.weight / 2.0
} else {
state.weight
}.max(state.weight / FACTOR)
.max(params.min_weight);
self.eps_weight = self.eps_weight.max(2.0 * cur_mod);
if new_weight < state.weight {
self.iter = 1;
self.model_max = NEG_INFINITY;
} else {
self.iter = max(self.iter + 1, 1);
}
debug!(" model_max={}", self.model_max);
debug!(" new_weight={}", new_weight);
new_weight
}
}
}