// Copyright (C) 2017  Davis E. King (davis@dlib.net)
// License: Boost Software License   See LICENSE.txt for the full license.
#undef DLIB_UPPER_bOUND_FUNCTION_ABSTRACT_Hh_
#ifdef DLIB_UPPER_bOUND_FUNCTION_ABSTRACT_Hh_

#include "../matrix.h"
#include <limits>

namespace dlib
{

// ----------------------------------------------------------------------------------------

    struct function_evaluation
    {
        /*!
            WHAT THIS OBJECT REPRESENTS
                This object records the output of a real valued function in response to
                some input. 

                In particular, if you have a function F(x) then the function_evaluation is
                simply a struct that records x and the scalar value F(x).
        !*/

        function_evaluation() = default;
        function_evaluation(const matrix<double,0,1>& x, double y) :x(x), y(y) {}

        matrix<double,0,1> x;
        double y = std::numeric_limits<double>::quiet_NaN();
    };

// ----------------------------------------------------------------------------------------

    class upper_bound_function
    {
        /*!
            WHAT THIS OBJECT REPRESENTS
                This object represents a non-parametric function that can be used to define
                an upper bound on some more complex and unknown function.  To describe this
                precisely, lets assume there is a function F(x) which you are capable of
                sampling from but otherwise know nothing about, and that you would like to
                find an upper bounding function U(x) such that U(x) >= F(x) for any x.  It
                would also be good if U(x)-F(x) was minimal.  I.e. we would like U(x) to be
                a tight upper bound, not something vacuous like U(x) = infinity.

                The upper_bound_function class is a tool for creating this kind of upper
                bounding function from a set of function_evaluations of F(x).  We do this
                by considering only U(x) of the form:
                    U = [](matrix<double,0,1> x) {
                       double min_ub = infinity;
                       for (size_t i = 0; i < POINTS.size(); ++i) {
                            function_evaluation p = POINTS[i]
                            double local_bound = p.y + sqrt(noise_terms[i] + trans(p.x-x)*M*(p.x-x))
                            min_ub = min(min_ub, local_bound)
                        }
                    }
                Where POINTS is an array of function_evaluation instances drawn from F(x),
                M is a diagonal matrix, and noise_terms is an array of scalars.

                To create an upper bound U(x), the upper_bound_function takes a POINTS array
                containing evaluations of F(x) as input and solves the following quadratic
                program to find the parameters of U(x):
                    
                    min_{M,noise_terms}:  sum(squared(M)) + sum(squared(noise_terms/relative_noise_magnitude))
                    s.t.   U(POINTS[i].x) >= POINTS[i].y,  for all i 
                           noise_terms[i] >= 0
                           min(M) >= 0
                           M is a diagonal matrix 
               
                Therefore, the quadratic program finds the U(x) that always upper bounds
                F(x) on the supplied POINTS, but is otherwise as small as possible.



                The inspiration for the upper_bound_function object came from the AdaLIPO
                algorithm from this excellent paper:
                    Global optimization of Lipschitz functions 
                    Malherbe, Cédric and Vayatis, Nicolas 
                    International Conference on Machine Learning - 2017
                In that paper, they propose to use a simpler U(x) where noise_terms is
                always 0 and M is a diagonal matrix where each diagonal element is the same
                value.  Therefore, there is only a single scalar parameter for U(x) in
                their formulation of the problem.  This causes difficulties if F(x) is
                stochastic or has discontinuities since, without the noise term, M will
                become really huge and the upper bound becomes vacuously large.  It is also
                problematic if the gradient of F(x) with respect to x contains elements of
                widely varying magnitude since the simpler formulation of U(x) assumes a
                uniform rate of change regardless of which dimension is varying. 
        !*/

    public:

        upper_bound_function(
        );
        /*!
            ensures
                - #num_points() == 0
                - #dimensionality() == 0
        !*/

        explicit upper_bound_function(
            const std::vector<function_evaluation>& points,
            const double relative_noise_magnitude = 0.001,
            const double solver_eps = 0.0001
        );
        /*!
            requires
                - all the x vectors in points must have the same non-zero dimensionality.
                - relative_noise_magnitude >= 0
                - solver_eps > 0
            ensures
                - Creates an upper bounding function U(x), as described above, assuming that
                  the given points are drawn from F(x).
                - Uses the provided relative_noise_magnitude when solving the QP, as
                  described above.  Note that relative_noise_magnitude can be set to 0.  If
                  you do this then all the noise terms are constrained to 0.  You should
                  only do this if you know F(x) is non-stochastic and continuous
                  everywhere.
                - When solving the QP used to find the parameters of U(x), the upper
                  bounding function, we solve the QP to solver_eps accuracy.   It's
                  possible that large enough solver_eps can lead to upper bounds that don't
                  upper bound all the supplied points.  But for reasonable epsilon values
                  this shouldn't be a problem. 
                - #num_points() == points.size()
                - #dimensionality() == points[0].x.size()
        !*/

        upper_bound_function(
            const double relative_noise_magnitude,
            const double solver_eps 
        );
        /*!
            requires
                - relative_noise_magnitude >= 0
                - solver_eps > 0
            ensures
                - #num_points() == 0
                - #dimensionality() == 0
                - This destructor is the same as calling the above constructor with points.size()==0
        !*/


        void add (
            const function_evaluation& point
        );
        /*!
            requires
                - num_points() == 0 || point.x.size() == dimensionality()
                - point.x.size() != 0
            ensures
                - Adds point to get_points().
                - Incrementally updates the upper bounding function with the given function
                  evaluation.  That is, we assume that F(point.x)==point.y and solve the QP
                  described above to find the new U(x) that upper bounds all the points
                  this object knows about (i.e. all the points in get_points() and the new point).
                - Calling add() is much faster than recreating the upper_bound_function
                  from scratch with all the points.  This is because we warm start with the
                  previous solution to the QP.  This is done by discarding any non-active
                  constraints and solving the QP again with only the previously active
                  constraints and the new constraints formed by all the pairs of the new
                  point and the old points.  This means the QP solved by add() is much
                  smaller than the QP that would be solved by a fresh call to the
                  upper_bound_function constructor.
        !*/

        const std::vector<function_evaluation>& get_points(
        ) const;
        /*!
            ensures
                - returns the points from F(x) used to define this upper bounding function.
                  These are all the function_evaluation objects given to this object via
                  its constructor and add().
        !*/

        long num_points(
        ) const;
        /*!
            ensures
                - returns the number of points used to define the upper bounding function.
                  (i.e. returns get_points().size())
        !*/

        long dimensionality(
        ) const;
        /*!
            ensures
                - returns the dimensionality of the input vectors to the upper bounding function. 
        !*/

        double operator() (
            const matrix<double,0,1>& x
        ) const;
        /*!
            requires
                - num_points() > 0
                - x.size() == dimensionality()
            ensures
                - return U(x)
                  (i.e. returns the upper bound on F(x) at x given by our upper bounding function)
        !*/

    };

// ----------------------------------------------------------------------------------------

}

#endif // DLIB_UPPER_bOUND_FUNCTION_ABSTRACT_Hh_


