particleswarmoptimization.hpp

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
 
/*
Copyright (C) 2015 Andres Hernandez
 
This file is part of QuantLib, a free-software/open-source library
for financial quantitative analysts and developers - http://quantlib.org/
 
QuantLib is free software: you can redistribute it and/or modify it
under the terms of the QuantLib license.  You should have received a
copy of the license along with this program; if not, please email
<quantlib-dev@lists.sf.net>. The license is also available online at
<http://quantlib.org/license.shtml>.
 
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 license for more details.
*/
 
#ifndef quantlib_optimization_particleswarmoptimization_hpp
#define quantlib_optimization_particleswarmoptimization_hpp
 
#include <ql/math/optimization/problem.hpp>
#include <ql/math/optimization/constraint.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
#include <ql/experimental/math/isotropicrandomwalk.hpp>
#include <ql/experimental/math/levyflightdistribution.hpp>
#include <ql/math/randomnumbers/seedgenerator.hpp>
 
#include <random>
 
namespace QuantLib {
 
    class ParticleSwarmOptimization : public OptimizationMethod {
      public:
        class Inertia;
        class Topology;
        ParticleSwarmOptimization(Size M,
                                  ext::shared_ptr<Topology> topology,
                                  ext::shared_ptr<Inertia> inertia,
                                  Real c1 = 2.05,
                                  Real c2 = 2.05,
                                  unsigned long seed = SeedGenerator::instance().get());
        explicit ParticleSwarmOptimization(Size M,
                                           ext::shared_ptr<Topology> topology,
                                           ext::shared_ptr<Inertia> inertia,
                                           Real omega,
                                           Real c1,
                                           Real c2,
                                           unsigned long seed = SeedGenerator::instance().get());
        void startState(Problem &P, const EndCriteria &endCriteria);
        EndCriteria::Type minimize(Problem& P, const EndCriteria& endCriteria) override;
 
      protected:
        std::vector<Array> X_, V_, pBX_, gBX_;
        Array pBF_, gBF_;
        Array lX_, uX_;
        Size M_, N_;
        Real c0_, c1_, c2_;
        MersenneTwisterUniformRng rng_;
        ext::shared_ptr<Topology> topology_;
        ext::shared_ptr<Inertia> inertia_;
    };
 
 
    class ParticleSwarmOptimization::Inertia {
        friend class ParticleSwarmOptimization;
      public:
        virtual ~Inertia() = default;
        virtual void setSize(Size M, Size N, Real c0, const EndCriteria &endCriteria) = 0;
        virtual void setValues() = 0;
      protected:
        ParticleSwarmOptimization *pso_;
        std::vector<Array> *X_, *V_, *pBX_, *gBX_;
        Array *pBF_, *gBF_;
        Array *lX_, *uX_;
 
        virtual void init(ParticleSwarmOptimization *pso) {
            pso_ = pso;
            X_ = &pso_->X_;
            V_ = &pso_->V_;
            pBX_ = &pso_->pBX_;
            gBX_ = &pso_->gBX_;
            pBF_ = &pso_->pBF_;
            gBF_ = &pso_->gBF_;
            lX_ = &pso_->lX_;
            uX_ = &pso_->uX_;
        }
    };
 
    /*     Inertia is a static value
    */
    class TrivialInertia : public ParticleSwarmOptimization::Inertia {
      public:
        inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
            c0_ = c0;
            M_ = M;
        }
        inline void setValues() override {
            for (Size i = 0; i < M_; i++) {
                (*V_)[i] *= c0_;
            }
        }
 
      private:
        Real c0_;
        Size M_;
    };
 
    /*     Inertia value gets multiplied with a random number
    between (threshold, 1)
    */
    class SimpleRandomInertia : public ParticleSwarmOptimization::Inertia {
      public:
        SimpleRandomInertia(Real threshold = 0.5, unsigned long seed = SeedGenerator::instance().get())
            : threshold_(threshold), rng_(seed) {
            QL_REQUIRE(threshold_ >= 0.0 && threshold_ < 1.0, "Threshold must be a Real in [0, 1)");
        }
        inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
            M_ = M;
            c0_ = c0;
        }
        inline void setValues() override {
            for (Size i = 0; i < M_; i++) {
                Real val = c0_*(threshold_ + (1.0 - threshold_)*rng_.nextReal());
                (*V_)[i] *= val;
            }
        }
 
      private:
        Real c0_, threshold_;
        Size M_;
        MersenneTwisterUniformRng rng_;
    };
 
    /*     Inertia value gets decreased every iteration until it reaches
    a value of threshold when iteration reaches the maximum level
    */
    class DecreasingInertia : public ParticleSwarmOptimization::Inertia {
      public:
        DecreasingInertia(Real threshold = 0.5)
            : threshold_(threshold) {
            QL_REQUIRE(threshold_ >= 0.0 && threshold_ < 1.0, "Threshold must be a Real in [0, 1)");
        }
        inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
            N_ = N;
            c0_ = c0;
            iteration_ = 0;
            maxIterations_ = endCriteria.maxIterations();
        }
        inline void setValues() override {
            Real c0 = c0_*(threshold_ + (1.0 - threshold_)*(maxIterations_ - iteration_) / maxIterations_);
            for (Size i = 0; i < M_; i++) {
                (*V_)[i] *= c0;
            }
        }
 
      private:
        Real c0_, threshold_;
        Size M_, N_, maxIterations_, iteration_;
    };
 
    /*    Alen Lukic, Approximating Kinetic Parameters Using Particle
    Swarm Optimization.
    */
    class AdaptiveInertia : public ParticleSwarmOptimization::Inertia {
      public:
        AdaptiveInertia(Real minInertia, Real maxInertia, Size sh = 5, Size sl = 2)
            :minInertia_(minInertia), maxInertia_(maxInertia),
            sh_(sh), sl_(sl) {};
        inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
            M_ = M;
            c0_ = c0;
            adaptiveCounter = 0;
            best_ = QL_MAX_REAL;
            started_ = false;
        }
        void setValues() override;
 
      private:
        Real c0_, best_;
        Real minInertia_, maxInertia_;
        Size M_;
        Size sh_, sl_;
        Size adaptiveCounter;
        bool started_;
    };
 
    /*    As long as the particle keeps getting frequent updates to its
    personal best value, the inertia behaves like a SimpleRandomInertia,
    but after a number of iterations without improvement, the behaviour
    changes to that of a Levy flight ~ u^{-1/\alpha}
    */
    class LevyFlightInertia : public ParticleSwarmOptimization::Inertia {
      public:
        LevyFlightInertia(Real alpha, Size threshold,
                          unsigned long seed = SeedGenerator::instance().get())
            :rng_(seed), generator_(seed), flight_(generator_, LevyFlightDistribution(1.0, alpha),
                1, Array(1, 1.0), seed),
            threshold_(threshold) {};
        inline void setSize(Size M, Size N, Real c0, const EndCriteria& endCriteria) override {
            M_ = M;
            N_ = N;
            c0_ = c0;
            adaptiveCounter_ = std::vector<Size>(M_, 0);
        }
        inline void setValues() override {
            for (Size i = 0; i < M_; i++) {
                if ((*pBF_)[i] < personalBestF_[i]) {
                    personalBestF_[i] = (*pBF_)[i];
                    adaptiveCounter_[i] = 0;
                }
                else {
                    adaptiveCounter_[i]++;
                }
                if (adaptiveCounter_[i] <= threshold_) {
                    //Simple Random Inertia
                    (*V_)[i] *= c0_*(0.5 + 0.5*rng_.nextReal());
                }
                else {
                    //If particle has not found a new personal best after threshold_ iterations
                    //then trigger a Levy flight pattern for the speed
                    flight_.nextReal<Real *>(&(*V_)[i][0]);
                }
            }
        }
 
      protected:
        void init(ParticleSwarmOptimization* pso) override {
            ParticleSwarmOptimization::Inertia::init(pso);
            personalBestF_ = *pBF_;
            flight_.setDimension(N_, *lX_, *uX_);
        }
 
      private:
        MersenneTwisterUniformRng rng_;
        std::mt19937 generator_;
        IsotropicRandomWalk<LevyFlightDistribution, std::mt19937> flight_;
        Array personalBestF_;
        std::vector<Size> adaptiveCounter_;
        Real c0_;
        Size M_, N_;
        Size threshold_;
    };
 
 
    class ParticleSwarmOptimization::Topology {
        friend class ParticleSwarmOptimization;
      public:
        virtual ~Topology() = default;
        virtual void setSize(Size M) = 0;
        virtual void findSocialBest() = 0;
      protected:
        ParticleSwarmOptimization *pso_;
        std::vector<Array> *X_, *V_, *pBX_, *gBX_;
        Array *pBF_, *gBF_;
      private:
        void init(ParticleSwarmOptimization *pso) {
            pso_ = pso;
            X_ = &pso_->X_;
            V_ = &pso_->V_;
            pBX_ = &pso_->pBX_;
            gBX_ = &pso_->gBX_;
            pBF_ = &pso_->pBF_;
            gBF_ = &pso_->gBF_;
        }
    };
 
    /*  The global best as seen by each particle is the best from amongst
    all particles
    */
    class GlobalTopology : public ParticleSwarmOptimization::Topology {
      public:
        inline void setSize(Size M) override { M_ = M; }
        inline void findSocialBest() override {
            Real bestF = (*pBF_)[0];
            Size bestP = 0;
            for (Size i = 1; i < M_; i++) {
                if (bestF < (*pBF_)[i]) {
                    bestF = (*pBF_)[i];
                    bestP = i;
                }
            }
            Array& x = (*pBX_)[bestP];
            for (Size i = 0; i < M_; i++) {
                if (i != bestP) {
                    (*gBX_)[i] = x;
                    (*gBF_)[i] = bestF;
                }
            }
        }
 
      private:
        Size M_;
    };
 
    /*  The global best as seen by each particle is the best from amongst
    the previous K and next K neighbors. For particle I, the best is
    then taken from amongst the [I - K, I + K] particles.
    */
    class KNeighbors : public ParticleSwarmOptimization::Topology {
      public:
        KNeighbors(Size K = 1) :K_(K) {
            QL_REQUIRE(K > 0, "Neighbors need to be larger than 0");
        }
        inline void setSize(Size M) override {
            M_ = M;
            QL_ENSURE(K_ < M, "Number of neighbors need to be smaller than total particles in swarm");
        }
        void findSocialBest() override;
 
      private:
        Size K_, M_;
    };
 
    /*  H.M. Emara,  Adaptive Clubs-based Particle Swarm Optimization
    Each particle is originally assigned to a default number of clubs
    from among the total set. The best as seen by each particle is the
    best from amongst the clubs to which the particle belongs.
    Underperforming particles join more clubs randomly (up to a maximum
    number) to widen the particles that influence them, while
    overperforming particles leave clubs randomly (down to a minimum
    number) to avoid early convergence to local minima.
    */
    class ClubsTopology : public ParticleSwarmOptimization::Topology {
      public:
        ClubsTopology(Size defaultClubs, Size totalClubs,
            Size maxClubs, Size minClubs,
            Size resetIteration, unsigned long seed = SeedGenerator::instance().get());
        void setSize(Size M) override;
        void findSocialBest() override;
 
      private:
        Size totalClubs_, maxClubs_, minClubs_, defaultClubs_;
        Size iteration_ = 0, resetIteration_;
        Size M_;
        std::vector<std::vector<bool> > clubs4particles_;
        std::vector<std::vector<bool> > particles4clubs_;
        std::vector<Size> bestByClub_;
        std::vector<Size> worstByClub_;
        std::mt19937 generator_;
        std::uniform_int_distribution<QuantLib::Size> distribution_;
        using param_type = decltype(distribution_)::param_type;
 
        void leaveRandomClub(Size particle, Size currentClubs);
        void joinRandomClub(Size particle, Size currentClubs);
    };
 
}
 
#endif