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