adaptive_sampling.h

#ifndef SOT_ADAPTIVE_SAMPLING_H
#define SOT_ADAPTIVE_SAMPLING_H

#include "common.h"
#include "utils.h"
#include "problem.h"
#include "genetic_algorithm.h"

namespace sot {
    

    class Sampling {
    public:

        virtual void reset(int budget) = 0;
        

        virtual mat makePoints(const vec &xBest, const mat &points, const vec &sigma, int newPoints) = 0;
    };
    

    template<class MeritFunction = MeritWeightedDistance>
    class SRBF : public Sampling {
    protected:
        std::shared_ptr<Problem> mData; 
        std::shared_ptr<Surrogate> mSurf; 
        int mNumCand; 
        int mDim; 
        vec mxLow; 
        vec mxUp; 
        double mDistTol; 
        int mNumEvals = 0; 
        int mBudget; 
        MeritFunction mMerit; 
    public:

        SRBF(const std::shared_ptr<Problem>& data, const std::shared_ptr<Surrogate>& surf, int numCand, int budget) {
            mData = std::shared_ptr<Problem>(data);
            mSurf = std::shared_ptr<Surrogate>(surf);
            mBudget = budget;
            mNumCand = numCand;
            mDim = data->dim();
            mxLow = data->lBounds();
            mxUp = data->uBounds();
            mDistTol = 1e-3*sqrt(arma::sum(arma::square(mxUp - mxLow)));
        }
        

        void reset(int budget) {
            mBudget = budget;
            mNumEvals = 0;
        }
        

        mat makePoints(const vec &xBest, const mat &points, const vec &sigma, int newPoints) {

            // Generate perturbations
            mat cand = arma::repmat(xBest, 1, mNumCand);
            mat pert = arma::randn(mDim, mNumCand);
            pert.each_col() %= sigma;
            cand += pert;
            
            for(int i=0; i < mNumCand; i++) {
                for(int j=0; j < mDim; j++) {
                    if(cand(j, i) > mxUp(j)) { // 
                        cand(j, i) = fmax(2*mxUp(j) - cand(j, i), mxLow(j)); 
                    }
                    else if(cand(j, i) < mxLow(j)) { 
                        cand(j, i) = fmin(2*mxLow(j) - cand(j, i), mxUp(j)); 
                    }
                }
            }
            
            // Update counter
            mNumEvals += newPoints;
            
            return mMerit.pickPoints(cand, mSurf, points, newPoints, mDistTol);
        }
    };
    

    template<class MeritFunction = MeritWeightedDistance>
    class DYCORS : public Sampling {
    protected:
        std::shared_ptr<Problem> mData; 
        std::shared_ptr<Surrogate> mSurf; 
        int mNumCand; 
        int mDim; 
        vec mxLow; 
        vec mxUp;  
        double mDistTol; 
        int mNumEvals = 0; 
        int mBudget; 
        MeritFunction mMerit; 
    public:

        DYCORS(const std::shared_ptr<Problem>& data, const std::shared_ptr<Surrogate>& surf, 
                int numCand, int budget) {
            mData = std::shared_ptr<Problem>(data);
            mSurf = std::shared_ptr<Surrogate>(surf);
            mBudget = budget;
            mNumCand = numCand;
            mDim = data->dim();
            mxLow = data->lBounds();
            mxUp = data->uBounds();
            mDistTol = 1e-3*sqrt(arma::sum(arma::square(mxUp - mxLow)));
        }
        

        void reset(int budget) {
            mBudget = budget;
            mNumEvals = 0;
        }
        

        mat makePoints(const vec &xBest, const mat &points, const vec &sigma, int newPoints) {                
            double dds_prob = std::min(20.0/mDim, 1.0) * 
                (1.0 - (std::log(mNumEvals + 1.0) / std::log(double(mBudget))));
            mat cand = arma::repmat(xBest, 1, mNumCand);
            
            for(int i=0; i < mNumCand; i++) {

                int count = 0;
                for(int j=0; j < mDim; j++) {
                    if(rand() < dds_prob) {
                        count++;
                        cand(j, i) += sigma(j)*randn();
                    }
                }
                // If no index was perturbed we force one
                if(count == 0) {
                    int ind = randi(mDim-1);
                    cand(ind, i) += sigma(ind)*randn();
                }

                // Make sure we are still in the domain
                for(int j=0; j < mDim; j++) {
                    if(cand(j, i) > mxUp(j)) { 
                        cand(j, i) = fmax(2*mxUp(j) - cand(j, i), mxLow(j)); 
                    }
                    else if(cand(j, i) < mxLow(j)) { 
                        cand(j, i) = fmin(2*mxLow(j) - cand(j, i), mxUp(j)); 
                    }
                }
            }
            
            // Update counter
            mNumEvals += newPoints;
            
            return mMerit.pickPoints(cand, mSurf, points, newPoints, mDistTol);
        }
    };
    

    template<class MeritFunction = MeritWeightedDistance>
    class Uniform : public Sampling {
    protected:
        std::shared_ptr<Problem> mData; 
        std::shared_ptr<Surrogate> mSurf; 
        int mNumCand; 
        int mDim; 
        vec mxLow; 
        vec mxUp;  
        double mDistTol; 
        int mNumEvals = 0; 
        int mBudget; 
        MeritFunction mMerit; 
    public:

        Uniform(const std::shared_ptr<Problem>& data, const std::shared_ptr<Surrogate>& surf, int numCand, int budget) {
            mData = std::shared_ptr<Problem>(data);
            mSurf = std::shared_ptr<Surrogate>(surf);
            mBudget = budget;
            mNumCand = numCand;
            mDim = data->dim();
            mxLow = data->lBounds();
            mxUp = data->uBounds();
            mDistTol = 1e-3*sqrt(arma::sum(arma::square(mxUp - mxLow)));
        }
        

        void reset(int budget) {
            mBudget = budget;
            mNumEvals = 0;
        }
        

        mat makePoints(const vec &xBest, const mat &points, const vec &sigma, int newPoints) {
         
            mat cand = arma::randu<mat>(mDim, mNumCand);
            for(int j=0; j < mDim; j++) {
                cand.row(j) = mxLow(j) + (mxUp(j) - mxLow(j)) * cand.row(j);
            }
        
            // Update counter
            mNumEvals += newPoints;
            
            return mMerit.pickPoints(cand, mSurf, points, newPoints, mDistTol);
        }
    };


    class GAWrapper : public Problem {
    protected:
        std::shared_ptr<Surrogate> mSurf; 
        int mDim; 
        mat mPoints; 
        vec mxLow; 
        vec mxUp; 
        vec mOptimum; 
        double mMinimum = 0; 
        std::string mName = "GA surrogate wrapper"; 
        double mDistTol; 
    public:

        GAWrapper(const std::shared_ptr<Problem> &data, const std::shared_ptr<Surrogate> &surf,
                const mat &points, double distTol) {
            mDim = data->dim();
            mxLow = data->lBounds();
            mxUp = data->uBounds();
            mSurf = std::shared_ptr<Surrogate>(surf);
            mPoints = points;
            mDistTol = distTol;
        }
        vec lBounds() const { return mxLow; }
        vec uBounds() const { return mxUp; }
        int dim() const { return mDim; }
        double min() const { return mMinimum; }
        vec optimum() const { return mOptimum; }
        std::string name() const { return mName; }
        double eval(const vec &x) const { return 0; }

        vec evals(const mat &X) const {

            mat dists = arma::sqrt(squaredPairwiseDistance<mat>(mPoints, X));
            // Evaluate the Surrogate at the points
            vec surfVals;

            if(mSurf->numPoints() == X.n_cols && mPoints.n_cols == X.n_cols) {
                surfVals = mSurf->evals(X, dists);
            }
            else {
                surfVals = mSurf->evals(X);
            }

            vec minDists = arma::min(dists).t();
            // Set the points that are too close to something large
            surfVals.elem(arma::find(minDists < mDistTol)).fill(arma::datum::inf);
            return surfVals;
        }
    };


    class GASampling : public Sampling {
    protected:
        std::shared_ptr<Problem> mData; 
        std::shared_ptr<Surrogate> mSurf; 
        int mDim; 
        vec mxLow; 
        vec mxUp;  
        int mNumIndividuals; 
        int mNumGenerations; 
        double mDistTol; 
    public:

        GASampling(const std::shared_ptr<Problem>& data, const std::shared_ptr<Surrogate>& surf,
                   int numIndividuals, int numGenerations) {
            mData = std::shared_ptr<Problem>(data);
            mSurf = std::shared_ptr<Surrogate>(surf);
            mDim = data->dim();
            mxLow = data->lBounds();
            mxUp = data->uBounds();
            mNumIndividuals = numIndividuals;
            mNumGenerations = numGenerations;
            mDistTol = 1e-3*sqrt(arma::sum(arma::square(mxUp - mxLow)));
        }
        
        void reset(int budget) {}

        mat makePoints(const vec &xBest, const mat &points, const vec &sigma, int newPoints) {
            if (newPoints > 1) {
                mat points2 = points;
                mat xNew = arma::zeros(mDim, newPoints);
                for(int i=0; i < newPoints; i++) {
                    std::shared_ptr<Problem> gaWrapper(
                            std::make_shared<GAWrapper>(mData, mSurf, points2, arma::norm(sigma)));
                    GeneticAlgorithm ga(gaWrapper, mNumIndividuals, mNumGenerations);
                    Result res = ga.run();
                    xNew.col(i) = res.xBest();
                    points2 = arma::join_horiz(points2, res.xBest());
                }
                return xNew;
            }
            else {
                std::shared_ptr<Problem> gaWrapper(
                        std::make_shared<GAWrapper>(mData, mSurf, points, arma::norm(sigma)));
                GeneticAlgorithm ga(gaWrapper, mNumIndividuals, mNumGenerations);
                Result res = ga.run();
                return res.xBest();
            }
        }
    };
}

#endif