dds.h

#ifndef SOT_DDS_H
#define SOT_DDS_H

#include <iostream>
#include "common.h"
#include "utils.h"
#include <thread>
#include <mutex>

namespace sot {


    class DDS {
    protected:
        std::shared_ptr<Problem> mData; 
        std::shared_ptr<ExpDesign> mExpDes; 
        int mMaxEvals; 
        int mNumEvals; 
        int mInitPoints; 
        int mDim; 
        vec mxLow; 
        vec mxUp; 
        std::string mName = "DDS"; 
        int mNumThreads; 
        int mEvalCount = 0; 
        std::mutex mMutex; 

        void evalBatch(const mat &batch, vec &funVals) {
            mMutex.lock();
            int myEval = mEvalCount;
            mEvalCount++;
            mMutex.unlock();
            
            while(myEval < batch.n_cols) {
                vec x = batch.col(myEval);
                funVals[myEval] = mData->eval(x);
                
                mMutex.lock();
                myEval = mEvalCount;
                mEvalCount++;
                mMutex.unlock();
            }
        }
    public:

        DDS(std::shared_ptr<Problem>& data, std::shared_ptr<ExpDesign>& expDes, int maxEvals) {
            mData = std::shared_ptr<Problem>(data);
            mExpDes = std::shared_ptr<ExpDesign>(expDes);
            mMaxEvals = maxEvals;
            mNumEvals = 0;
            mInitPoints = expDes->numPoints();
            mDim = data->dim();
            mxLow = data->lBounds();
            mxUp = data->uBounds();
            mName = "DDS";
            mNumThreads = 1;
            if(mMaxEvals < mInitPoints) { 
                throw std::logic_error("Experimental design larger than evaluation budget"); 
            }
        }

        DDS(std::shared_ptr<Problem>& data, std::shared_ptr<ExpDesign>& expDes, int maxEvals, int numThreads) 
            : DDS(data, expDes, maxEvals) {
            mNumThreads = numThreads;
        }
        

        Result run() {
            std::vector<std::thread> threads(mNumThreads);
            Result res(mMaxEvals, mDim);
            mNumEvals = 0;
            
            vec sigma = 0.2*(mxUp - mxLow);
            mat initDes = fromUnitBox(mExpDes->generatePoints(), mxLow, mxUp);
            vec initFunVal = arma::zeros(mInitPoints);

            if(mNumThreads > 1) { // Evaluate in synchronous parallel
                mEvalCount = 0;            
                for(int i=0; i < mNumThreads; i++) {
                    threads[i] = std::thread(&sot::DDS::evalBatch, this, 
                            std::ref(initDes), std::ref(initFunVal));
                }

                for(int i=0; i < mNumThreads; i++) {
                    threads[i].join();
                }
            }
            else { // Evaluate in serial
                for(int i=0; i < mInitPoints; i++) {
                    vec x = initDes.col(i);
                    initFunVal(i) = mData->eval(x);
                }
            }
            
            res.addEvals(initDes, initFunVal);
            mNumEvals += mInitPoints;
            
            while (mNumEvals < mMaxEvals) {
                
                double ddsProb = 1 - std::log((double)mNumEvals - mInitPoints)/
                    std::log((double)mMaxEvals - mInitPoints);
                ddsProb = fmax(ddsProb, 1.0/mDim);
                
                int newEvals = std::min<int>(mNumThreads, mMaxEvals - mNumEvals);
                mat batch = arma::zeros<mat>(mDim, newEvals);
                vec batchVals = arma::zeros(newEvals);  

                for(int i=0; i < newEvals; i++) {
                    vec cand = res.xBest();
                    int count = 0;
                    for(int j=0; j < mDim; j++) {
                        if(rand() < ddsProb) {
                            count++;
                            cand(j) += sigma(j) * randn();
                            if(cand(j) > mxUp(j)) { 
                                cand(j) = fmax(2*mxUp(j) - cand(j), mxLow(j)); 
                            }
                            else if(cand(j) < mxLow(j)) { 
                                cand(j) = fmin(2*mxLow(j) - cand(j), mxUp(j)); 
                            }
                        }
                    }
                    // If no index was perturbed we force one
                    if(count == 0) {
                        int ind = randi(mDim-1);
                        cand(ind) += sigma(ind) * randn();
                        if(cand(ind) > mxUp(ind)) { 
                            cand(ind) = fmax(2*mxUp(ind) - cand(ind), mxLow(ind)); 
                        }
                        else if(cand(ind) < mxLow(ind)) { 
                            cand(ind) = fmin(2*mxLow(ind) - cand(ind), mxUp(ind)); 
                        }
                    }
                    batch.col(i) = cand;
                }
                
                               
                if(newEvals > 1) { // Evaluate in synchronous parallel
                    mEvalCount = 0;
                    for(int i=0; i < newEvals; i++) {
                        threads[i] = std::thread(&sot::DDS::evalBatch, this, 
                                std::ref(batch), std::ref(batchVals));
                    }

                    for(int i=0; i < newEvals; i++) {
                        threads[i].join();
                    }
                }
                else { // Evaluate in serial
                    batchVals(0) = mData->eval((vec)batch);
                }
                      
                // Update evaluation counter
                mNumEvals += newEvals;                
                res.addEvals(batch, batchVals);                
            }
                                
            return res;
        }
    };
}

#endif