//===========================================================================
/*!
 *
 *
 * \brief       GridSearch.h
 *
 *
 *
 * \author      O. Krause
 * \date        2010
 *
 *
 * \par Copyright 1995-2015 Shark Development Team
 *
 * <BR><HR>
 * This file is part of Shark.
 * <http://image.diku.dk/shark/>
 *
 * Shark is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * Shark 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 Lesser General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with Shark.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
//===========================================================================

#ifndef SHARK_ALGORITHMS_GRIDSEARCH_H
#define SHARK_ALGORITHMS_GRIDSEARCH_H


#include <Rcpp.h>
#include <shark/Algorithms/AbstractSingleObjectiveOptimizer.h>
#include <shark/Rng/GlobalRng.h>



namespace shark {

//!
//! \brief Optimize by trying out a grid of configurations
//!
//! \par
//! The GridSearch class allows for the definition of a grid in
//! parameter space. It does a simple one-step optimization over
//! the grid by trying out every possible parameter combination.
//! Please note that the computation effort grows exponentially
//! with the number of parameters.
//!
//! \par
//! If you only want to try a subset of the grid, consider using
//! the PointSearch class instead.
//! A more sophisticated (less exhaustive) grid search variant is
//! available with the NestedGridSearch class.
//!
class GridSearch : public AbstractSingleObjectiveOptimizer<RealVector >
{
public:
        GridSearch() {
                m_configured=false;
        }

        /// \brief From INameable: return the class name.
        std::string name() const
        {
                return "GridSearch";
        }

        //! uniform initialization for all parameters
        //! \param  params  number of model parameters to optimize
        //! \param  min     smallest parameter value
        //! \param  max     largest parameter value
        //! \param  numSections   total number of values in the interval
        void configure(size_t params, double min, double max, size_t numSections)
        {
                SIZE_CHECK(params>=1);
                RANGE_CHECK(min<=max);
                SIZE_CHECK(numSections>=1);
                m_nodeValues.resize(params);
                for (size_t i = 0; i < numSections; i++)
                {
                        double section = min + i * (max - min) / (numSections - 1.0);
                        for(auto& node: m_nodeValues)
                                node.push_back(section);
                }
                m_configured=true;
        }

        //! individual definition for every parameter
        //! \param  min     smallest value for every parameter
        //! \param  max     largest value for every parameter
        //! \param  sections   total number of values for every parameter
        void configure(const std::vector<double>& min, const std::vector<double>& max, const std::vector<size_t>& sections)
        {
                size_t params = min.size();
                SIZE_CHECK(sections.size() == params);
                SIZE_CHECK(max.size() == params);
                RANGE_CHECK(min <= max);

                m_nodeValues.resize(params);
                for (size_t dimension = 0; dimension < params; dimension++)
                {
                        size_t numSections = sections[dimension];
                        std::vector<double>& node = m_nodeValues[dimension];
                        node.resize(numSections);

                        if ( numSections == 1 )
                        {
                                node[0] = (( min[dimension] + max[dimension] ) / 2.0);
                        }
                        else for (size_t section = 0; section < numSections; section++)
                                {
                                        node[section] = (min[dimension] + section * (max[dimension] - min[dimension]) / (numSections - 1.0));
                                }
                }
                m_configured=true;
        }


        //! special case for 2D grid, individual definition for every parameter
        //! \param  min1     smallest value for first parameter
        //! \param  max1     largest value for first parameter
        //! \param  sections1   total number of values for first parameter
        //! \param  min2     smallest value for second parameter
        //! \param  max2     largest value for second parameter
        //! \param  sections2   total number of values for second parameter
        void configure(double min1, double max1, size_t sections1, double min2, double max2, size_t sections2)
        {
                RANGE_CHECK(min1<=max1);
                RANGE_CHECK(min2<=max2);
                RANGE_CHECK(sections1 > 0);
                RANGE_CHECK(sections2 > 0);

                m_nodeValues.resize(2u);

                if ( sections1 == 1 ) {
                        m_nodeValues[0].push_back(( min1 + max1 ) / 2.0);
                } else {
                        for (size_t section = 0; section < sections1; section++)
                                m_nodeValues[0].push_back(min1 + section * (max1 - min1) / (sections1 - 1.0));
                }

                if ( sections2 == 1 ) {
                        m_nodeValues[1].push_back(( min2 + max2 ) / 2.0);
                } else {
                        for (size_t section = 0; section < sections2; section++)
                                m_nodeValues[1].push_back(min2 + section * (max2 - min2) / (sections2 - 1.0));
                }
        }

        //! special case for line search
        //! \param  min1     smallest value for first parameter
        //! \param  max1     largest value for first parameter
        //! \param  sections1   total number of values for first parameter
        void configure(double min1, double max1, size_t sections1)
        {
                RANGE_CHECK(min1<=max1);
                RANGE_CHECK(sections1 > 0);

                m_nodeValues.resize(1u);

                if ( sections1 == 1 ) {
                        m_nodeValues[0].push_back(( min1 + max1 ) / 2.0);
                } else {
                        for (size_t section = 0; section < sections1; section++)
                                m_nodeValues[0].push_back(min1 + section * (max1 - min1) / (sections1 - 1.0));
                }
        }


        //! uniform definition of the values to test for all parameters
        //! \param  params  number of model parameters to optimize
        //! \param  values  values used for every coordinate
        void configure(size_t params, const std::vector<double>& values)
        {
                SIZE_CHECK(params > 0);
                SIZE_CHECK(values.size() > 0);
                m_nodeValues.resize(params);
                for(auto& node: m_nodeValues)
                        node=values;
                m_configured=true;
        }

        //! individual definition for every parameter
        //! \param  values  values used. The first dimension is the parameter, the second dimension is the node.
        void configure(const std::vector<std::vector<double> >& values)
        {
                SIZE_CHECK(values.size() > 0);
                m_nodeValues = values;
                m_configured=true;
        }

        //from ISerializable
        /*! If Gridsearch wasn't configured before calling this method, it is default constructed
         *  as a net spanning the range [-1,1] in all dimensions with 5 searchpoints (-1,-0.5,0,0.5,1).
         *  so don't forget to scale the parameter-ranges of the objective function!
         *  The startingPoint can actually be anything, only its dimension has to be correct.
         */
        virtual void init(ObjectiveFunctionType & objectiveFunction, SearchPointType const& startingPoint) {
                checkFeatures(objectiveFunction);

                if(!m_configured)
                        configure(startingPoint.size(),-1,1,5);
                SIZE_CHECK(startingPoint.size() == m_nodeValues.size());
                m_best.point=startingPoint;
        }
        using AbstractSingleObjectiveOptimizer<RealVector >::init;
        /*! Assign linearly progressing grid values to one certain parameter only.
         *  This is especially useful if one parameter needs special treatment
         *  \param index the index of the parameter to which grid values are assigned
         *  \param noOfSections how many grid points should be assigned to that parameter
         *  \param min smallest value for that parameter
         *  \param max largest value for that parameter */
        void assignLinearRange(size_t index, size_t noOfSections, double min, double max)
        {
                SIZE_CHECK( noOfSections >= 1);
                RANGE_CHECK( min <= max );
                SIZE_CHECK( index < m_nodeValues.size() );

                m_nodeValues[index].clear();
                if ( noOfSections == 1 ) {
                        m_nodeValues[index].push_back(( min+max) / 2.0);
                }
                else {
                        m_nodeValues[index].reserve(noOfSections);
                        for (size_t section = 0; section < noOfSections; section++)
                                m_nodeValues[index].push_back(min + section*( max-min ) / ( noOfSections-1.0 ));
                }
        }

        /*! Set exponentially progressing grid values for one certain parameter only.
         *  This is especially useful if one parameter needs special treatment. The
         *  grid points will be filled with values \f$ factor \cdot expbase ^i \f$,
         *  where i does integer steps between min and max.
         *  \param index the index of the parameter that gets new grid values
         *  \param factor the value that the exponential base grid should be multiplied by
         *  \param exp_base the exponential grid will progress on this base (e.g. 2, 10)
         *  \param min the smallest exponent for exp_base
         *  \param max the largest exponent for exp_base  */
        void assignExponentialRange(size_t index, double factor, double exp_base, int min, int max)
        {
                SIZE_CHECK( min <= max );
                RANGE_CHECK( index < m_nodeValues.size() );

                m_nodeValues[index].clear();
                m_nodeValues[index].reserve(max-min);
                for (int section = 0; section <= (max-min); section++)
                        m_nodeValues[index].push_back( factor * std::pow( exp_base, section+min ));
        }

        //! Please note that for the grid search optimizer it does
        //! not make sense to call step more than once, as the
        //! solution does not improve iteratively.
        void step(ObjectiveFunctionType const& objectiveFunction) {
                size_t dimensions = m_nodeValues.size();
                std::vector<size_t> index(dimensions, 0);
                m_best.value = 1e100;
                RealVector point(dimensions);

                // loop through all grid points
                while (true)
                {
                        // define the parameters
                        for (size_t dimension = 0; dimension < dimensions; dimension++)
                                point(dimension) = m_nodeValues[dimension][index[dimension]];

                        // evaluate the model
                        if (objectiveFunction.isFeasible(point))
                        {
                                double error = objectiveFunction.eval(point);

#ifdef SHARK_CV_VERBOSE_1
                                Rcpp::Rcout << "." << std::flush;
#endif
#ifdef SHARK_CV_VERBOSE
                                Rcpp::Rcout << point << "\t" << error << std::endl;
#endif
                                if (error < m_best.value)
                                {
                                        m_best.value = error;
                                        m_best.point = point;           // [TG] swap() solution is out, caused ugly memory bug, I changed this back
                                }
                        }

                        // next index
                        size_t dimension = 0;
                        for (; dimension < dimensions; dimension++)
                        {
                                index[dimension]++;
                                if (index[dimension] < m_nodeValues[dimension].size()) break;
                                index[dimension] = 0;
                        }
                        if (dimension == dimensions) break;
                }
#ifdef SHARK_CV_VERBOSE_1
                Rcpp::Rcout << std::endl;
#endif
        }

protected:

        //! The array columns contain the grid values for the corresponding parameter axis.
        std::vector<std::vector<double> > m_nodeValues;

        bool m_configured;
};


//!
//! \brief Nested grid search
//!
//! \par
//! The NestedGridSearch class is an iterative optimizer,
//! doing one grid search in every iteration. In every
//! iteration, it halves the grid extent doubling the
//! resolution in every coordinate.
//!
//! \par
//! Although nested grid search is much less exhaustive
//! than standard grid search, it still suffers from
//! exponential time and memory complexity in the number
//! of variables optimized. Therefore, if the number of
//! variables is larger than 2 or 3, consider using the
//! CMA instead.
//!
//! \par
//! Nested grid search works as follows: The optimizer
//! defined a 5x5x...x5 equi-distant grid (depending on
//! the search space dimension) on an initially defined
//! search cube. During every grid search iteration,
//! the error is computed for all  grid points.
//!Then the grid is moved
//! to the best grid point found so far and contracted
//! by a factor of two in each dimension. Each call to
//! the optimize() function performs one such step.
//!
//! \par
//! Let N denote the number of parameters to optimize.
//! To compute the error landscape at the current
//! zoom level, the algorithm has to do \f$ 5^N \f$
//! error function evaluations in every iteration.
//!
//! \par
//! The grid is always centered around the best
//! solution currently known. If this solution is
//! located at the boundary, the landscape may exceed
//! the parameter range defined m_minimum and m_maximum.
//! These invalid landscape values are not used.
//!
class NestedGridSearch : public AbstractSingleObjectiveOptimizer<RealVector >
{
public:
        //! Constructor
        NestedGridSearch()
        {
                m_configured=false;
        }

        /// \brief From INameable: return the class name.
        std::string name() const
        {
                return "NestedGridSearch";
        }

        //!
        //! \brief Initialization of the nested grid search.
        //!
        //! \par
        //! The min and max arrays define ranges for every parameter to optimize.
        //! These ranges are strict, that is, the algorithm will not try values
        //! beyond the range, even if is finds a boundary minimum.
        //!
        //! \param  min    lower end of the parameter range
        //! \param  max    upper end of the parameter range
        void configure(const std::vector<double>& min, const std::vector<double>& max)
        {
                size_t dimensions = min.size();
                SIZE_CHECK(max.size() == dimensions);

                m_minimum = min;
                m_maximum = max;

                m_stepsize.resize(dimensions);
                m_best.point.resize(dimensions);
                for (size_t dimension = 0; dimension < dimensions; dimension++)
                {
                        m_best.point(dimension)=0.5 *(min[dimension] + max[dimension]);
                        m_stepsize[dimension] = 0.25 * (max[dimension] - min[dimension]);
                }
                m_configured=true;
        }

        //!
        //! \brief Initialization of the nested grid search.
        //!
        //! \par
        //! The min and max values define ranges for every parameter to optimize.
        //! These ranges are strict, that is, the algorithm will not try values
        //! beyond the range, even if is finds a boundary minimum.
        //!
        //! \param parameters number of parameters to optimize
        //! \param  min    lower end of the parameter range
        //! \param  max    upper end of the parameter range
        void configure(size_t parameters, double min, double max)
        {
                SIZE_CHECK(parameters > 0);

                m_minimum=std::vector<double>(parameters,min);
                m_maximum=std::vector<double>(parameters,max);
                m_stepsize=std::vector<double>(parameters,0.25 * (max - min));

                m_best.point.resize(parameters);

                double start=0.5 *(min + max);
                for(double& value: m_best.point)
                        value=start;
                m_configured=true;
        }

        //from ISerializable
        //! if NestedGridSearch was not configured before this call, it is default initialized ti the range[-1,1] for every parameter
        virtual void init(ObjectiveFunctionType & objectiveFunction, SearchPointType const& startingPoint) {
                objectiveFunction.init();
                checkFeatures(objectiveFunction);

                if(!m_configured)
                        configure(startingPoint.size(),-1,1);
                SIZE_CHECK(m_stepsize.size()==startingPoint.size());

        }
        using AbstractSingleObjectiveOptimizer<RealVector >::init;


        //! Every call of the optimization member computes the
        //! error landscape on the current grid. It picks the
        //! best error value and zooms into the error landscape
        //! by a factor of 2.
        void step(ObjectiveFunctionType const& objectiveFunction) {
                size_t dimensions = m_stepsize.size();
                SIZE_CHECK(m_minimum.size() == dimensions);
                SIZE_CHECK(m_maximum.size() == dimensions);

                m_best.value = 1e99;
                std::vector<size_t> index(dimensions,0);

                RealVector point=m_best.point;

                // loop through the grid
                while (true)
                {
                        // compute the grid point,

                        // set the parameters
                        bool compute=true;
                        for (size_t d = 0; d < dimensions; d++)
                        {
                                point(d) += (index[d] - 2.0) * m_stepsize[d];
                                if (point(d) < m_minimum[d] || point(d) > m_maximum[d])
                                {
                                        compute = false;
                                        break;
                                }
                        }

                        // evaluate the grid point
                        if (compute && objectiveFunction.isFeasible(point))
                        {
                                double error = objectiveFunction.eval(point);

                                // remember the best solution
                                if (error < m_best.value)
                                {
                                        m_best.value = error;
                                        m_best.point=point;
                                }
                        }
                        // move to the next grid point
                        size_t d = 0;
                        for (; d < dimensions; d++)
                        {
                                index[d]++;
                                if (index[d] <= 4) break;
                                index[d] = 0;
                        }
                        if (d == dimensions) break;
                }
                // decrease the step sizes
                for(double& step: m_stepsize)
                        step *= 0.5;
        }

protected:
        //! minimum parameter value to check
        std::vector<double> m_minimum;

        //! maximum parameter value to check
        std::vector<double> m_maximum;

        //! current step size for every parameter
        std::vector<double> m_stepsize;

        bool m_configured;
};


//!
//! \brief Optimize by trying out predefined configurations
//!
//! \par
//! The PointSearch class is similair to the GridSearch class
//! by the property that it optimizes a model in a single pass
//! just trying out a predefined number of parameter configurations.
//! The main difference is that every parameter configuration has
//! to be explicitly defined. It is not possible to define a set
//! of values for every axis; see GridSearch for this purpose.
//! Thus, the PointSearch class allows for more flexibility.
//!
//! If no configure method is called, this class just samples random points.
//! They are uniformly distributed in [-1,1].
//! parameters^2 points but minimum 20 are sampled in this case.
//!
class PointSearch : public AbstractSingleObjectiveOptimizer<RealVector >
{
public:
        //! Constructor
        PointSearch() {
                m_configured=false;
        }

        /// \brief From INameable: return the class name.
        std::string name() const
        {
                return "PointSearch";
        }

        //! Initialization of the search points.
        //! \param  points  array of points to evaluate
        void configure(const std::vector<RealVector>& points) {
                m_points=points;
                m_configured=true;
        }

        //!samples random points in the range [min,max]^parameters
        void configure(size_t parameters,size_t samples, double min,double max) {
                RANGE_CHECK(min<=max);
                m_points.resize(samples);
                for(size_t sample=0; sample!=samples; ++sample)
                {
                        m_points[sample].resize(parameters);
                        for(size_t param=0; param!=parameters; ++param)
                        {
                                m_points[sample](param)=Rng::uni(min,max);
                        }
                }
                m_configured=true;
        }

        //! If the class wasn't configured before, this method samples random uniform distributed points in [-1,1]^n.
        void init(ObjectiveFunctionType& objectiveFunction, SearchPointType const& startingPoint) {
                objectiveFunction.init();
                checkFeatures(objectiveFunction);

                if(!m_configured)
                {
                        size_t parameters=startingPoint.size();
                        size_t samples=std::min(sqr(parameters),(size_t)20);
                        configure(parameters,samples,-1,1);
                }
        }
        using AbstractSingleObjectiveOptimizer<RealVector >::init;
        //! Please note that for the point search optimizer it does
        //! not make sense to call step more than once, as the
        //! solution does not improve iteratively.
        void step(ObjectiveFunctionType const& objectiveFunction) {
                size_t numPoints = m_points.size();
                m_best.value = 1e100;
                size_t bestIndex=0;

                // loop through all points
                for (size_t point = 0; point < numPoints; point++)
                {
                        // evaluate the model
                        if (objectiveFunction.isFeasible(m_points[point]))
                        {
                                double error = objectiveFunction.eval(m_points[point]);
                                if (error < m_best.value)
                                {
                                        m_best.value = error;
                                        bestIndex=point;
                                }
                        }
                }
                m_best.point=m_points[bestIndex];
        }

protected:
        //! The array holds one parameter configuration in every column.
        std::vector<RealVector> m_points;

        //! verbosity level
        bool m_configured;
};


}
#endif

