// Copyright (c) 2016-2020 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/>
//
//! A linear minorant.
use super::{Aggregatable, DVector, Real};
use num_traits::Zero;
use std::borrow::Borrow;
use crate::data::raw::BLAS1;
/// A linear minorant of a convex function.
///
/// A linear minorant of a convex function $f \colon \mathbb{R}\^n \to
/// \mathbb{R}$ is a linear function of the form
///
/// \\[ l \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto \langle g, x
/// \rangle + c \\]
///
/// such that $l(x) \le f(x)$ for all $x \in \mathbb{R}\^n$.
///
/// Each minorant may have an additional "primal" data, that is
/// linearly combined along with the minorant itself. This data can
/// be used for separation algorithms.
pub trait Minorant: Aggregatable + Clone + Default + Send + 'static {
type Primal: Aggregatable + Send;
/// Return the constant value.
fn constant(&self) -> Real;
/// Return the dimension of the linear term.
fn dim(&self) -> usize;
/// Evaluate the linear term at some point.
fn linear(&self, x: &DVector) -> Real;
/// Return the associated primal data.
fn primal(&self) -> &Self::Primal;
/// Compute the inner product with another minorant.
fn product(&self, other: &Self) -> Real;
/// Add a scaled `self.linear` to `target`.
///
/// This sets `target = target + alpha * self.linear`.
fn add_scaled_to(&self, alpha: Real, target: &mut [Real]);
/// Compute linear combination of (linear terms of) the given minorants and
/// store in a dense vector.
fn combine_to_vec<I, A>(minorants: I, target: &mut (Real, DVector))
where
I: IntoIterator<Item = (Real, A)>,
A: Borrow<Self>,
{
target.0 = Real::zero();
let mut it = minorants.into_iter();
if let Some((alpha, m)) = it.next() {
let m = m.borrow();
target.1.init0(m.dim());
m.add_scaled_to(alpha, &mut target.1);
target.0 += alpha * m.constant();
for (alpha, m) in it {
let m = m.borrow();
m.add_scaled_to(alpha, &mut target.1);
target.0 += alpha * m.constant();
}
} else {
target.1.clear();
}
}
/// Evaluate minorant at some point.
///
/// This function computes $c + \langle g, x \rangle$ for this minorant
/// \\[\ell \colon \mathbb{R}\^n \to \mathbb{R}, x \mapsto c + \langle g, x \rangle\\]
/// and the given point $x \in \mathbb{R}\^n$.
///
fn eval(&self, x: &DVector) -> Real {
self.constant() + self.linear(x)
}
/// Move the center of the minorant along direction `d`.
fn move_center(&mut self, alpha: Real, d: &DVector) {
debug_assert_eq!(d.len(), self.dim(), "Minorant and direction must have same dimension");
self.move_center_begin(alpha, d)
}
/// Move the center of the minorant.
///
/// In contrast to [`move_center`] the vector `d` might be shorter than the
/// linear term of the minorant. The missing elements are assumed to be
/// zero.
///
/// [`move_center`]: Minorant::move_center
fn move_center_begin(&mut self, alpha: Real, d: &DVector);
fn debug(&self) -> &dyn std::fmt::Debug;
}
impl<P: Aggregatable + Clone + Send + 'static> Minorant for (Real, DVector, P) {
type Primal = P;
fn dim(&self) -> usize {
self.1.len()
}
fn constant(&self) -> Real {
self.0
}
fn linear(&self, x: &DVector) -> Real {
self.1.dot(x)
}
fn primal(&self) -> &P {
&self.2
}
fn product(&self, other: &Self) -> Real {
self.1.dot(&other.1)
}
fn add_scaled_to(&self, alpha: Real, target: &mut [Real]) {
debug_assert_eq!(self.1.len(), target.len());
if !alpha.is_zero() {
BLAS1::add_scaled(target, alpha, &self.1);
}
}
fn move_center(&mut self, alpha: Real, d: &DVector) {
self.0 += alpha * self.1.dot(d);
}
fn move_center_begin(&mut self, alpha: Real, d: &DVector) {
debug_assert!(self.1.len() >= d.len());
self.0 += alpha * self.1.dot_begin(d);
}
fn debug(&self) -> &dyn std::fmt::Debug {
&self.1
}
}
impl Minorant for (Real, DVector) {
type Primal = ();
fn dim(&self) -> usize {
self.1.len()
}
fn constant(&self) -> Real {
self.0
}
fn linear(&self, x: &DVector) -> Real {
self.1.dot(x)
}
fn primal(&self) -> &Self::Primal {
&()
}
fn product(&self, other: &Self) -> Real {
self.1.dot(&other.1)
}
fn add_scaled_to(&self, alpha: Real, target: &mut [Real]) {
debug_assert_eq!(self.1.len(), target.len());
if !alpha.is_zero() {
BLAS1::add_scaled(target, alpha, &self.1);
}
}
fn move_center(&mut self, alpha: Real, d: &DVector) {
self.0 += alpha * self.1.dot(d);
}
fn move_center_begin(&mut self, alpha: Real, d: &DVector) {
debug_assert!(self.1.len() >= d.len());
self.0 += alpha * self.1.dot_begin(d);
}
fn debug(&self) -> &dyn std::fmt::Debug {
self
}
}
/// The subgradient extender is a callback used to update an existing minorant.
pub type SubgradientExtender<M, E> = dyn FnMut(usize, &mut M) -> Result<(), E> + Send;