DifferentialEvolution_MT Class Reference

#include <qle/math/differentialevolution_mt.hpp>

Inheritance diagram for DifferentialEvolution_MT:

Collaboration diagram for DifferentialEvolution_MT:

Public Types
using	Strategy = DifferentialEvolution::Strategy
using	CrossoverType = DifferentialEvolution::CrossoverType
using	Candidate = DifferentialEvolution::Candidate
using	Configuration = DifferentialEvolution::Configuration

Public Member Functions
	DifferentialEvolution_MT (Configuration configuration=Configuration(), std::string maxTime="")
virtual EndCriteria::Type	minimize (Problem_MT &p, const EndCriteria &endCriteria) override
	minimize the optimization problem P More...
const Configuration &	configuration () const
Public Member Functions inherited from OptimizationMethod_MT
virtual	~OptimizationMethod_MT ()
virtual EndCriteria::Type	minimize (Problem_MT &P, const EndCriteria &endCriteria)=0
	minimize the optimization problem P More...
EndCriteria::Type	minimize (Problem &, const EndCriteria &) override

Private Member Functions
void	updateCost (std::vector< Candidate > &population, Problem_MT &p) const
bool	checkMaxTime () const
void	fillInitialPopulation (std::vector< Candidate > &population, const Problem_MT &p) const
void	getCrossoverMask (std::vector< Array > &crossoverMask, std::vector< Array > &invCrossoverMask, const Array &mutationProbabilities) const
Array	getMutationProbabilities (const std::vector< Candidate > &population) const
void	adaptSizeWeights () const
void	adaptCrossover () const
void	calculateNextGeneration (std::vector< Candidate > &population, Problem_MT &p) const
Array	rotateArray (Array inputArray) const
void	crossover (const std::vector< Candidate > &oldPopulation, std::vector< Candidate > &population, const std::vector< Candidate > &mutantPopulation, const std::vector< Candidate > &mirrorPopulation, Problem_MT &p) const

Private Attributes
Configuration	configuration_
std::string	maxTime_
Array	upperBound_
Array	lowerBound_
Array	currGenSizeWeights_
Array	currGenCrossover_
Candidate	bestMemberEver_
MersenneTwisterUniformRng	rng_

Detailed Description

Definition at line 54 of file differentialevolution_mt.hpp.

Member Typedef Documentation

Strategy

using Strategy = DifferentialEvolution::Strategy

Definition at line 56 of file differentialevolution_mt.hpp.

CrossoverType

using CrossoverType = DifferentialEvolution::CrossoverType

Definition at line 57 of file differentialevolution_mt.hpp.

Candidate

using Candidate = DifferentialEvolution::Candidate

Definition at line 58 of file differentialevolution_mt.hpp.

Configuration

using Configuration = DifferentialEvolution::Configuration

Definition at line 59 of file differentialevolution_mt.hpp.

Constructor & Destructor Documentation

DifferentialEvolution_MT()

DifferentialEvolution_MT	(	Configuration	configuration = Configuration(),
		std::string	maxTime = ""
	)

Definition at line 61 of file differentialevolution_mt.hpp.

        : configuration_(configuration), maxTime_(maxTime), rng_(configuration.seed) {}

Member Function Documentation

minimize()

EndCriteria::Type minimize ( Problem_MT &  P, const EndCriteria &  endCriteria  )

overridevirtual

minimize the optimization problem P

Implements OptimizationMethod_MT.

Definition at line 72 of file differentialevolution_mt.cpp.

                                                                                                {
    EndCriteria::Type ecType;
 
    if (configuration().upperBound.empty()) {
        upperBound_ = p.constraint().upperBound(p.currentValue());
    } else {
        QL_REQUIRE(configuration().upperBound.size() == p.currentValue().size(),
                   "wrong upper bound size in differential evolution configuration");
        upperBound_ = configuration().upperBound;
    }
    if (configuration().lowerBound.empty()) {
        lowerBound_ = p.constraint().lowerBound(p.currentValue());
    } else {
        QL_REQUIRE(configuration().lowerBound.size() == p.currentValue().size(),
                   "wrong lower bound size in differential evolution configuration");
        lowerBound_ = configuration().lowerBound;
    }
 
    currGenSizeWeights_ = Array(configuration().populationMembers, configuration().stepsizeWeight);
    currGenCrossover_ = Array(configuration().populationMembers, configuration().crossoverProbability);
 
    std::vector<Candidate> population;
    if (!configuration().initialPopulation.empty()) {
        population.resize(configuration().initialPopulation.size());
        for (Size i = 0; i < population.size(); ++i) {
            population[i].values = configuration().initialPopulation[i];
            QL_REQUIRE(population[i].values.size() == p.currentValue().size(),
                       "wrong values size in initial population");
        }
    } else {
        population = std::vector<Candidate>(configuration().populationMembers, Candidate(p.currentValue().size()));
        fillInitialPopulation(population, p);
    }
 
    updateCost(population, p);
 
    std::partial_sort(population.begin(), population.begin() + 1, population.end(), sort_by_cost());
    bestMemberEver_ = population.front();
    Real fxOld = population.front().cost;
    Size iteration = 0, stationaryPointIteration = 0;
 
    // main loop - calculate consecutive emerging populations
    while (!endCriteria.checkMaxIterations(iteration++, ecType) && !checkMaxTime()) {
        calculateNextGeneration(population, p);
        std::partial_sort(population.begin(), population.begin() + 1, population.end(), sort_by_cost());
        if (population.front().cost < bestMemberEver_.cost)
            bestMemberEver_ = population.front();
        Real fxNew = population.front().cost;
        if (endCriteria.checkStationaryFunctionValue(fxOld, fxNew, stationaryPointIteration, ecType))
            break;
        fxOld = fxNew;
    };
    p.setCurrentValue(bestMemberEver_.values);
    p.setFunctionValue(bestMemberEver_.cost);
    if (checkMaxTime())
        ecType = EndCriteria::Type::Unknown;
    return ecType;
}

Here is the call graph for this function:

configuration()

const Configuration & configuration

(

)

const

Definition at line 66 of file differentialevolution_mt.hpp.

{ return configuration_; }

Here is the caller graph for this function:

updateCost()

void updateCost ( std::vector< Candidate > &  population, Problem_MT &  p  ) const

private

Definition at line 304 of file differentialevolution_mt.cpp.

                                                                                               {
    struct Worker {
        Worker(std::vector<Candidate>& p, const Size s, const Size e, const QuantLib::ext::shared_ptr<CostFunction> c)
            : population(p), start(s), end(e), costFunction(c) {}
        void operator()() {
            for (Size popIter = start; popIter < end; popIter++) {
                // evaluate objective function as soon as possible to avoid unnecessary loops
                try {
                    population[popIter].cost = costFunction->value(population[popIter].values);
                } catch (Error&) {
                    population[popIter].cost = QL_MAX_REAL;
                }
                if (!std::isfinite(population[popIter].cost))
                    population[popIter].cost = QL_MAX_REAL;
            }
        }
        std::vector<Candidate>& population;
        const Size start, end;
        const QuantLib::ext::shared_ptr<CostFunction> costFunction;
    };
 
    Size threads = p.costFunctions().size();
    QL_REQUIRE(threads > 0, "DifferentialEvolution_MT: number of available threads is zero");
 
    std::vector<Size> chunks(threads, std::max<Size>(population.size() / threads, 1));
    int rest = population.size() - threads * chunks[0];
    do {
        Size i = 0;
        while (i < chunks.size() && rest > 0) {
            chunks[i]++;
            rest--;
            ++i;
        }
    } while (rest > 0);
 
    std::vector<QuantLib::ext::shared_ptr<std::thread>> workers(threads);
 
    Size start = 0, end;
    for (Size thread = 0; thread < threads; ++thread) {
        end = std::min(start + chunks[thread], population.size());
        Worker worker(population, start, end, p.costFunctions()[thread]);
        workers[thread] = QuantLib::ext::make_shared<std::thread>(worker);
        start = end;
    }
 
    for (Size thread = 0; thread < threads; ++thread) {
        workers[thread]->join();
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

checkMaxTime()

bool checkMaxTime ( ) const

private

Definition at line 64 of file differentialevolution_mt.cpp.

                                                  {
    if (maxTime_.empty())
        return false;
    QL_REQUIRE(maxTime_.length() == 15, "maxTime (" << maxTime_ << ") must have format YYYYMMDDTHHMMSS");
    std::string ts = to_iso_string(boost::posix_time::microsec_clock::local_time()).substr(0, 15);
    return ts > maxTime_;
}

Here is the caller graph for this function:

fillInitialPopulation()

void fillInitialPopulation ( std::vector< Candidate > &  population, const Problem_MT &  p  ) const

private

Definition at line 417 of file differentialevolution_mt.cpp.

                                                                                                                {
 
    // use initial values provided by the user
    population.front().values = p.currentValue();
    // rest of the initial population is random
    for (Size j = 1; j < population.size(); ++j) {
        for (Size i = 0; i < p.currentValue().size(); ++i) {
            Real l = lowerBound_[i], u = upperBound_[i];
            population[j].values[i] = l + (u - l) * rng_.nextReal();
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

getCrossoverMask()

void getCrossoverMask ( std::vector< Array > &  crossoverMask, std::vector< Array > &  invCrossoverMask, const Array &  mutationProbabilities  ) const

private

Definition at line 354 of file differentialevolution_mt.cpp.

                                                                                          {
    for (Size cmIter = 0; cmIter < crossoverMask.size(); cmIter++) {
        for (Size memIter = 0; memIter < crossoverMask[cmIter].size(); memIter++) {
            if (rng_.nextReal() < mutationProbabilities[cmIter]) {
                invCrossoverMask[cmIter][memIter] = 0.0;
            } else {
                crossoverMask[cmIter][memIter] = 0.0;
            }
        }
    }
}

Here is the caller graph for this function:

getMutationProbabilities()

Array getMutationProbabilities ( const std::vector< Candidate > &  population ) const

private

Definition at line 367 of file differentialevolution_mt.cpp.

                                                                                                     {
    Array mutationProbabilities = currGenCrossover_;
    switch (configuration().crossoverType) {
    case DifferentialEvolution::Normal:
        break;
    case DifferentialEvolution::Binomial:
        mutationProbabilities =
            currGenCrossover_ * (1.0 - 1.0 / population.front().values.size()) + 1.0 / population.front().values.size();
        break;
    case DifferentialEvolution::Exponential:
        for (Size coIter = 0; coIter < currGenCrossover_.size(); coIter++) {
            mutationProbabilities[coIter] =
                (1.0 - std::pow(currGenCrossover_[coIter], (int)population.front().values.size())) /
                (population.front().values.size() * (1.0 - currGenCrossover_[coIter]));
        }
        break;
    default:
        QL_FAIL("Unknown crossover type (" << Integer(configuration().crossoverType) << ")");
        break;
    }
    return mutationProbabilities;
}

Here is the call graph for this function:

Here is the caller graph for this function:

adaptSizeWeights()

void adaptSizeWeights ( ) const

private

Definition at line 395 of file differentialevolution_mt.cpp.

                                                      {
    // [=Fl & =Fu] respectively see Brest, J. et al., 2006,
    // "Self-Adapting Control Parameters in Differential
    // Evolution"
    Real sizeWeightLowerBound = 0.1, sizeWeightUpperBound = 0.9;
    // [=tau1] A Comparative Study on Numerical Benchmark
    // Problems." page 649 for reference
    Real sizeWeightChangeProb = 0.1;
    for (Size coIter = 0; coIter < currGenSizeWeights_.size(); coIter++) {
        if (rng_.nextReal() < sizeWeightChangeProb)
            currGenSizeWeights_[coIter] = sizeWeightLowerBound + rng_.nextReal() * sizeWeightUpperBound;
    }
}

Here is the caller graph for this function:

adaptCrossover()

void adaptCrossover ( ) const

private

Definition at line 409 of file differentialevolution_mt.cpp.

                                                    {
    Real crossoverChangeProb = 0.1; // [=tau2]
    for (Size coIter = 0; coIter < currGenCrossover_.size(); coIter++) {
        if (rng_.nextReal() < crossoverChangeProb)
            currGenCrossover_[coIter] = rng_.nextReal();
    }
}

Here is the caller graph for this function:

calculateNextGeneration()

void calculateNextGeneration ( std::vector< Candidate > &  population, Problem_MT &  p  ) const

private

Definition at line 131 of file differentialevolution_mt.cpp.

                                                                                                            {
 
    std::vector<Candidate> mirrorPopulation;
    std::vector<Candidate> oldPopulation = population;
 
    switch (configuration().strategy) {
 
    case QuantLib::DifferentialEvolution::Rand1Standard: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop2 = population;
        randomize(population.begin(), population.end(), rng_);
        mirrorPopulation = shuffledPop1;
 
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            population[popIter].values =
                population[popIter].values +
                configuration().stepsizeWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
        }
    } break;
 
    case QuantLib::DifferentialEvolution::BestMemberWithJitter: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
        Array jitter(population[0].values.size(), 0.0);
 
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            for (double& jitterIter : jitter) {
                jitterIter = rng_.nextReal();
            }
            population[popIter].values =
                bestMemberEver_.values + (shuffledPop1[popIter].values - population[popIter].values) *
                                             (0.0001 * jitter + configuration().stepsizeWeight);
        }
        mirrorPopulation = std::vector<Candidate>(population.size(), bestMemberEver_);
    } break;
 
    case QuantLib::DifferentialEvolution::CurrentToBest2Diffs: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
 
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            population[popIter].values =
                oldPopulation[popIter].values +
                configuration().stepsizeWeight * (bestMemberEver_.values - oldPopulation[popIter].values) +
                configuration().stepsizeWeight * (population[popIter].values - shuffledPop1[popIter].values);
        }
        mirrorPopulation = shuffledPop1;
    } break;
 
    case QuantLib::DifferentialEvolution::Rand1DiffWithPerVectorDither: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop2 = population;
        randomize(population.begin(), population.end(), rng_);
        mirrorPopulation = shuffledPop1;
        Array FWeight = Array(population.front().values.size(), 0.0);
        for (double& fwIter : FWeight)
            fwIter = (1.0 - configuration().stepsizeWeight) * rng_.nextReal() + configuration().stepsizeWeight;
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            population[popIter].values =
                population[popIter].values + FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
        }
    } break;
 
    case QuantLib::DifferentialEvolution::Rand1DiffWithDither: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop2 = population;
        randomize(population.begin(), population.end(), rng_);
        mirrorPopulation = shuffledPop1;
        Real FWeight = (1.0 - configuration().stepsizeWeight) * rng_.nextReal() + configuration().stepsizeWeight;
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            population[popIter].values =
                population[popIter].values + FWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
        }
    } break;
 
    case QuantLib::DifferentialEvolution::EitherOrWithOptimalRecombination: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop2 = population;
        randomize(population.begin(), population.end(), rng_);
        mirrorPopulation = shuffledPop1;
        Real probFWeight = 0.5;
        if (rng_.nextReal() < probFWeight) {
            for (Size popIter = 0; popIter < population.size(); popIter++) {
                population[popIter].values =
                    oldPopulation[popIter].values +
                    configuration().stepsizeWeight * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
            }
        } else {
            Real K = 0.5 * (configuration().stepsizeWeight + 1); // invariant with respect to probFWeight used
            for (Size popIter = 0; popIter < population.size(); popIter++) {
                population[popIter].values =
                    oldPopulation[popIter].values + K * (shuffledPop1[popIter].values - shuffledPop2[popIter].values -
                                                         2.0 * population[popIter].values);
            }
        }
    } break;
 
    case QuantLib::DifferentialEvolution::Rand1SelfadaptiveWithRotation: {
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop1 = population;
        randomize(population.begin(), population.end(), rng_);
        std::vector<Candidate> shuffledPop2 = population;
        randomize(population.begin(), population.end(), rng_);
        mirrorPopulation = shuffledPop1;
 
        adaptSizeWeights();
 
        for (Size popIter = 0; popIter < population.size(); popIter++) {
            if (rng_.nextReal() < 0.1) {
                population[popIter].values = rotateArray(bestMemberEver_.values);
            } else {
                population[popIter].values =
                    bestMemberEver_.values +
                    currGenSizeWeights_[popIter] * (shuffledPop1[popIter].values - shuffledPop2[popIter].values);
            }
        }
    } break;
 
    default:
        QL_FAIL("Unknown strategy (" << Integer(configuration().strategy) << ")");
    }
    // in order to avoid unnecessary copying we use the same population object for mutants
    crossover(oldPopulation, population, population, mirrorPopulation, p);
}

Here is the call graph for this function:

Here is the caller graph for this function:

rotateArray()

Array rotateArray ( Array  inputArray ) const

private

Definition at line 390 of file differentialevolution_mt.cpp.

                                                         {
    randomize(a.begin(), a.end(), rng_);
    return a;
}

Here is the caller graph for this function:

crossover()

void crossover ( const std::vector< Candidate > &  oldPopulation, std::vector< Candidate > &  population, const std::vector< Candidate > &  mutantPopulation, const std::vector< Candidate > &  mirrorPopulation, Problem_MT &  p  ) const

private

Definition at line 266 of file differentialevolution_mt.cpp.

                                                                                                            {
 
    if (configuration().crossoverIsAdaptive) {
        adaptCrossover();
    }
 
    Array mutationProbabilities = getMutationProbabilities(population);
 
    std::vector<Array> crossoverMask(population.size(), Array(population.front().values.size(), 1.0));
    std::vector<Array> invCrossoverMask = crossoverMask;
    getCrossoverMask(crossoverMask, invCrossoverMask, mutationProbabilities);
 
    // crossover of the old and mutant population
 
    for (Size popIter = 0; popIter < population.size(); popIter++) {
        population[popIter].values = oldPopulation[popIter].values * invCrossoverMask[popIter] +
                                     mutantPopulation[popIter].values * crossoverMask[popIter];
        // immediately apply bounds if specified
        if (configuration().applyBounds) {
            for (Size memIter = 0; memIter < population[popIter].values.size(); memIter++) {
                if (population[popIter].values[memIter] > upperBound_[memIter])
                    population[popIter].values[memIter] =
                        upperBound_[memIter] +
                        rng_.nextReal() * (mirrorPopulation[popIter].values[memIter] - upperBound_[memIter]);
                if (population[popIter].values[memIter] < lowerBound_[memIter])
                    population[popIter].values[memIter] =
                        lowerBound_[memIter] +
                        rng_.nextReal() * (mirrorPopulation[popIter].values[memIter] - lowerBound_[memIter]);
            }
        }
    }
 
    updateCost(population, p);
}

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

configuration_

Configuration configuration_

private

Definition at line 71 of file differentialevolution_mt.hpp.

maxTime_

std::string maxTime_

private

Definition at line 72 of file differentialevolution_mt.hpp.

upperBound_

Array upperBound_

private

Definition at line 73 of file differentialevolution_mt.hpp.

lowerBound_

Array lowerBound_

private

Definition at line 73 of file differentialevolution_mt.hpp.

currGenSizeWeights_

Array currGenSizeWeights_

mutableprivate

Definition at line 74 of file differentialevolution_mt.hpp.

currGenCrossover_

Array currGenCrossover_

private

Definition at line 74 of file differentialevolution_mt.hpp.

bestMemberEver_

Candidate bestMemberEver_

private

Definition at line 75 of file differentialevolution_mt.hpp.

rng_

MersenneTwisterUniformRng rng_

private

Definition at line 76 of file differentialevolution_mt.hpp.