fireflyalgorithm.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_fireflyalgorithm_hpp
#define quantlib_optimization_fireflyalgorithm_hpp
 
#include <ql/math/optimization/problem.hpp>
#include <ql/math/optimization/constraint.hpp>
#include <ql/experimental/math/isotropicrandomwalk.hpp>
#include <ql/experimental/math/levyflightdistribution.hpp>
#include <ql/math/randomnumbers/mt19937uniformrng.hpp>
#include <ql/math/randomnumbers/seedgenerator.hpp>
 
#include <cmath>
#include <random>
 
namespace QuantLib {
 
    class FireflyAlgorithm : public OptimizationMethod {
      public:
        class RandomWalk;
        class Intensity;
        FireflyAlgorithm(Size M,
                         ext::shared_ptr<Intensity> intensity,
                         ext::shared_ptr<RandomWalk> randomWalk,
                         Size Mde = 0,
                         Real mutationFactor = 1.0,
                         Real crossoverFactor = 0.5,
                         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_, xI_, xRW_; 
        std::vector<std::pair<Real, Size> > values_;
        Array lX_, uX_;
        Real mutation_, crossover_;
        Size M_, N_, Mde_, Mfa_;
        ext::shared_ptr<Intensity> intensity_;
        ext::shared_ptr<RandomWalk> randomWalk_;
        std::mt19937 generator_;
        std::uniform_int_distribution<QuantLib::Size> distribution_;
        MersenneTwisterUniformRng rng_;
    };
 
 
    class FireflyAlgorithm::Intensity {
        friend class FireflyAlgorithm;
    public:
      virtual ~Intensity() = default;
      void findBrightest();
    protected:
        Size Mfa_, N_;
        const std::vector<Array> *x_;
        const std::vector<std::pair<Real, Size> > *values_;
        std::vector<Array> *xI_;
 
        virtual Real intensityImpl(Real valueX, Real valueY, Real distance) = 0;
        inline Real distance(const Array& x, const Array& y) const {
            Real d = 0.0;
            for (Size i = 0; i < N_; i++) {
                Real diff = x[i] - y[i];
                d += diff*diff;
            }
            return d;
        }
 
    private:
        void init(FireflyAlgorithm *fa) {
            x_ = &fa->x_;
            xI_ = &fa->xI_;
            values_ = &fa->values_;
            Mfa_ = fa->Mfa_;
            N_ = fa->N_;
        }
    };
 
    /*  Exponentially decreasing intensity
    */
    class ExponentialIntensity : public FireflyAlgorithm::Intensity {
      public:
          ExponentialIntensity(Real beta0, Real betaMin, Real gamma)
              : beta0_(beta0), betaMin_(betaMin), gamma_(gamma) {}
      protected:
        Real intensityImpl(Real valueX, Real valueY, Real d) override {
            return (beta0_ - betaMin_) * std::exp(-gamma_ * d) + betaMin_;
        }
          Real beta0_, betaMin_, gamma_;
    };
 
    /*  Inverse law square
    */
    class InverseLawSquareIntensity : public FireflyAlgorithm::Intensity {
    public:
        InverseLawSquareIntensity(Real beta0, Real betaMin)
            : beta0_(beta0), betaMin_(betaMin) {}
    protected:
      Real intensityImpl(Real valueX, Real valueY, Real d) override {
          return (beta0_ - betaMin_) / (d + QL_EPSILON) + betaMin_;
      }
        Real beta0_, betaMin_;
    };
 
 
    class FireflyAlgorithm::RandomWalk {
        friend class FireflyAlgorithm;
    public:
      virtual ~RandomWalk() = default;
      void walk() {
          for (Size i = 0; i < Mfa_; i++) {
              walkImpl((*xRW_)[(*values_)[i].second]);
          }
        }
    protected:
        Size Mfa_, N_;
        const std::vector<Array> *x_;
        const std::vector<std::pair<Real, Size> > *values_;
        std::vector<Array> *xRW_;
        Array *lX_, *uX_;
 
        virtual void walkImpl(Array & xRW) = 0;
        virtual void init(FireflyAlgorithm *fa) {
            x_ = &fa->x_;
            xRW_ = &fa->xRW_;
            values_ = &fa->values_;
            Mfa_ = fa->Mfa_;
            N_ = fa->N_;
            lX_ = &fa->lX_;
            uX_ = &fa->uX_;
        }
    };
 
    /*  Random walk given by distribution template parameter. The
        distribution must be compatible with std::mt19937.
    */
    template <class Distribution>
    class DistributionRandomWalk : public FireflyAlgorithm::RandomWalk {
      public:
        explicit DistributionRandomWalk(Distribution dist, 
                                        Real delta = 0.9, 
                                        unsigned long seed = SeedGenerator::instance().get())
        : walkRandom_(std::mt19937(seed), std::move(dist), 1, Array(1, 1.0), seed),
          delta_(delta) {}
      protected:
        void walkImpl(Array& xRW) override {
            walkRandom_.nextReal(&xRW[0]);
            xRW *= delta_;
        }
        void init(FireflyAlgorithm* fa) override {
            FireflyAlgorithm::RandomWalk::init(fa);
            walkRandom_.setDimension(N_, *lX_, *uX_);
        }
        IsotropicRandomWalk<Distribution, std::mt19937> walkRandom_;
        Real delta_;
    };
    
    /*  Gaussian random walk
    */
    class GaussianWalk : public DistributionRandomWalk<std::normal_distribution<QuantLib::Real>> {
      public:
        explicit GaussianWalk(Real sigma, 
                              Real delta = 0.9, 
                              unsigned long seed = SeedGenerator::instance().get())
        : DistributionRandomWalk<std::normal_distribution<QuantLib::Real>>(
                           std::normal_distribution<QuantLib::Real>(0.0, sigma), delta, seed){}
    };
 
    /*  Levy flight random walk
    */
    class LevyFlightWalk : public DistributionRandomWalk<LevyFlightDistribution> {
      public:
        explicit LevyFlightWalk(Real alpha, 
                                Real xm = 0.5, 
                                Real delta = 0.9,
                                unsigned long seed = SeedGenerator::instance().get())
        : DistributionRandomWalk<LevyFlightDistribution>(
                            LevyFlightDistribution(xm, alpha), delta, seed) {}
    };
 
    /*  Gaussian random walk, but size of step decreases with each iteration step
    */
    class DecreasingGaussianWalk : public GaussianWalk {
      public:
        explicit DecreasingGaussianWalk(
            Real sigma,
            Real delta = 0.9,
            unsigned long seed = SeedGenerator::instance().get())
        : GaussianWalk(sigma, delta, seed), delta0_(delta) {}
      protected:
        void walkImpl(Array& xRW) override {
            iteration_++;
            if (iteration_ > Mfa_) {
                //Every time all the fireflies have been processed
                //multiply delta by itself
                iteration_ = 0;
                delta_ *= delta_;
            }
            GaussianWalk::walkImpl(xRW);
        }
        void init(FireflyAlgorithm* fa) override {
            GaussianWalk::init(fa);
            iteration_ = 0;
            delta_ = delta0_;
        }
 
      private:
        Real delta0_;
        Size iteration_;
    };
}
 
#endif