13097754652

Committed 02 Feb 2025 10:03AM UTC coverage: 73.25% (-0.003%) from 73.253%

Build # 13097754652

Build Type

Pull #2150

github

Committed by

web-flow

Commit Message

Merge 77d357318 into 3a61fcd13

Pull Request Pull Request #2150: added china calendar for the year 2025

Run Details

5 of 5 new or added lines in 1 file covered. (100.0%)

3 existing lines in 2 files now uncovered.

56142 of 76644 relevant lines covered (73.25%)

8680278.19 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

97.17

/ql/experimental/math/fireflyalgorithm.cpp

/* -*- 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.
*/

#include <ql/experimental/math/fireflyalgorithm.hpp>
#include <ql/math/randomnumbers/sobolrsg.hpp>
#include <algorithm>
#include <cmath>
#include <utility>

namespace QuantLib {
    FireflyAlgorithm::FireflyAlgorithm(Size M,
                                       ext::shared_ptr<Intensity> intensity,
                                       ext::shared_ptr<RandomWalk> randomWalk,
                                       Size Mde,
                                       Real mutation,
                                       Real crossover,
                                       unsigned long seed)
    : mutation_(mutation), crossover_(crossover), M_(M), Mde_(Mde), Mfa_(M_ - Mde_),
      intensity_(std::move(intensity)), randomWalk_(std::move(randomWalk)),
      generator_(seed), distribution_(Mfa_, Mde > 0 ? M_ - 1 : M_),
      rng_(seed) {
        QL_REQUIRE(M_ >= Mde_,
            "Differential Evolution subpopulation cannot be larger than total population");
    }

    void FireflyAlgorithm::startState(Problem &P, const EndCriteria &endCriteria) {
        N_ = P.currentValue().size();
        x_.reserve(M_);
        xI_.reserve(M_);
        xRW_.reserve(M_);
        values_.reserve(M_);
        uX_ = P.constraint().upperBound(P.currentValue());
        lX_ = P.constraint().lowerBound(P.currentValue());
        Array bounds = uX_ - lX_;

        //Random initialization is done by Sobol sequence
        SobolRsg sobol(N_);

        //Prepare containers
        for (Size i = 0; i < M_; i++) {
            const SobolRsg::sample_type::value_type &sample = sobol.nextSequence().value;
            x_.emplace_back(N_, 0.0);
            xI_.emplace_back(N_, 0.0);
            xRW_.emplace_back(N_, 0.0);
            Array& x = x_.back();
            for (Size j = 0; j < N_; j++) {
                //Assign X=lb+(ub-lb)*random
                x[j] = lX_[j] + bounds[j] * sample[j];
            }
            //Evaluate point
            values_.emplace_back(P.value(x), i);
        }

        //init intensity & randomWalk
        intensity_->init(this);
        randomWalk_->init(this);
    }

    EndCriteria::Type FireflyAlgorithm::minimize(Problem &P, const EndCriteria &endCriteria) {
        QL_REQUIRE(!P.constraint().empty(), "Firefly Algorithm is a constrained optimizer");
        EndCriteria::Type ecType = EndCriteria::None;
        P.reset();
        Size iteration = 0;
        Size iterationStat = 0;
        Size maxIteration = endCriteria.maxIterations();
        Size maxIStationary = endCriteria.maxStationaryStateIterations();
        
        startState(P, endCriteria);

        bool isFA = Mfa_ > 0;
        //Variables for DE
        Array z(N_, 0.0);
        Size indexR1, indexR2;
        decltype(distribution_)::param_type nParam(0, N_ - 1);

        //Set best value & position
        Real bestValue = values_[0].first;
        Size bestPosition = 0;
        for (Size i = 1; i < M_; i++) {
            if (values_[i].first < bestValue) {
                bestPosition = i;
                bestValue = values_[i].first;
            }
        }
        Array bestX = x_[bestPosition];

        //Run optimization
        do {
            iteration++;
            iterationStat++;
            //Check if stopping criteria is met
            if (iteration > maxIteration || iterationStat > maxIStationary)
                break;

            //Divide into two subpopulations
            //First sort values
            std::sort(values_.begin(), values_.end());

            //Differential evolution
            if(Mfa_ < M_){
                Size indexBest = values_[0].second;
                Array& xBest = x_[indexBest];
                for (Size i = Mfa_; i < M_; i++) { 
                    if (!isFA) {
                        //Pure DE requires random index
                        indexBest = distribution_(generator_);
                        xBest = x_[indexBest];
                    }
                    do { 
                        indexR1 = distribution_(generator_);
                    } while(indexR1 == indexBest);
                    do { 
                        indexR2 = distribution_(generator_);
                    } while(indexR2 == indexBest || indexR2 == indexR1);
                    
                    Size index = values_[i].second;
                    Array& x   = x_[index];
                    Array& xR1 = x_[indexR1];
                    Array& xR2 = x_[indexR2];
                                        Size rIndex = distribution_(generator_, nParam);
                    for (Size j = 0; j < N_; j++) {
                        if (j == rIndex || rng_.nextReal() <= crossover_) {
                            //Change x[j] according to crossover
                            z[j] = xBest[j] + mutation_*(xR1[j] - xR2[j]);
                        } else {
                            z[j] = x[j];
                        }
                        //Enforce bounds on positions
                        if (z[j] < lX_[j]) {
                            z[j] = lX_[j];
                        }
                        else if (z[j] > uX_[j]) {
                            z[j] = uX_[j];
                        }
                    }
                    Real val = P.value(z);
                    if (val < values_[index].first) {
                        //Accept new point
                        x = z;
                        values_[index].first = val;
                        //mark best
                        if (val < bestValue) {
                            bestValue = val;
                            bestX = x;
                            iterationStat = 0;
                        }
                    }
                }
            }
                
            //Firefly algorithm
            if(isFA){
                //According to the intensity, determine best global position
                intensity_->findBrightest();

                //Prepare random walk
                randomWalk_->walk();

                //Loop over particles
                for (Size i = 0; i < Mfa_; i++) {
                    Size index = values_[i].second;
                    Array& x   = x_[index];
                    Array& xI  = xI_[index];
                    Array& xRW = xRW_[index];

                    //Loop over dimensions
                    for (Size j = 0; j < N_; j++) {
                        //Update position
                        z[j] = x[j] + xI[j] + xRW[j];
                        //Enforce bounds on positions
                        if (z[j] < lX_[j]) {
                            z[j] = lX_[j];
                        }
                        else if (z[j] > uX_[j]) {
                            z[j] = uX_[j];
                        }
                    }
                    Real val = P.value(z);
                    if(!std::isnan(val))
                                        {
                                                //Accept new point
                        x = z;
                        values_[index].first = val;
                        //mark best
                        if (val < bestValue) {
                            bestValue = val;
                            bestX = x;
                            iterationStat = 0;
                        }
                                        }
                }
            }
        } while (true);
        if (iteration > maxIteration)
            ecType = EndCriteria::MaxIterations;
        else
            ecType = EndCriteria::StationaryPoint;

        //Set result to best point
        P.setCurrentValue(bestX);
        P.setFunctionValue(bestValue);
        return ecType;
    }

    void FireflyAlgorithm::Intensity::findBrightest() {
        //Brightest ignores all others
        Array& xI = (*xI_)[(*values_)[0].second];
        for (Size j = 0; j < N_; j++) {
            xI[j] = 0.0;
        }

        for (Size i = 1; i < Mfa_; i++) {
            //values_ is already sorted
            Size index = (*values_)[i].second;
            const Array& x = (*x_)[index];
            Array& xI = (*xI_)[index];
            for (Size j = 0; j < N_; j++) {
                xI[j] = 0.0;
            }
            Real valueX = (*values_)[i].first;
            for (Size k = 0; k < i - 1; k++){
                const Array& y = (*x_)[(*values_)[k].second];
                Real valueY = (*values_)[k].first;
                Real intensity = intensityImpl(valueX, valueY, distance(x, y));
                xI += intensity*(y - x);
            }
        }
    }
}


1	/* -- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -- */
2
3	/*
4	Copyright (C) 2015 Andres Hernandez
5
6	This file is part of QuantLib, a free-software/open-source library
7	for financial quantitative analysts and developers - http://quantlib.org/
8
9	QuantLib is free software: you can redistribute it and/or modify it
10	under the terms of the QuantLib license. You should have received a
11	copy of the license along with this program; if not, please email
12	<quantlib-dev@lists.sf.net>. The license is also available online at
13	<http://quantlib.org/license.shtml>.
14
15	This program is distributed in the hope that it will be useful, but WITHOUT
16	ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
17	FOR A PARTICULAR PURPOSE. See the license for more details.
18	*/
19
20	#include <ql/experimental/math/fireflyalgorithm.hpp>
21	#include <ql/math/randomnumbers/sobolrsg.hpp>
22	#include <algorithm>
23	#include <cmath>
24	#include <utility>
25
26	namespace QuantLib {
27	FireflyAlgorithm::FireflyAlgorithm(Size M,	1✔
28	ext::shared_ptr<Intensity> intensity,
29	ext::shared_ptr<RandomWalk> randomWalk,
30	Size Mde,
31	Real mutation,
32	Real crossover,
33	unsigned long seed)	1✔
34	: mutation_(mutation), crossover_(crossover), M_(M), Mde_(Mde), Mfa_(M_ - Mde_),	1✔
35	intensity_(std::move(intensity)), randomWalk_(std::move(randomWalk)),
36	generator_(seed), distribution_(Mfa_, Mde > 0 ? M_ - 1 : M_),	1✔
37	rng_(seed) {	2✔
38	QL_REQUIRE(M_ >= Mde_,	1✔
39	"Differential Evolution subpopulation cannot be larger than total population");
40	}	1✔
41
42	void FireflyAlgorithm::startState(Problem &P, const EndCriteria &endCriteria) {	1✔
43	N_ = P.currentValue().size();	1✔
44	x_.reserve(M_);	1✔
45	xI_.reserve(M_);	1✔
46	xRW_.reserve(M_);	1✔
47	values_.reserve(M_);	1✔
48	uX_ = P.constraint().upperBound(P.currentValue());	1✔
49	lX_ = P.constraint().lowerBound(P.currentValue());	1✔
50	Array bounds = uX_ - lX_;	1✔
51
52	//Random initialization is done by Sobol sequence
53	SobolRsg sobol(N_);	1✔
54
55	//Prepare containers
56	for (Size i = 0; i < M_; i++) {	151✔
57	const SobolRsg::sample_type::value_type &sample = sobol.nextSequence().value;	150✔
58	x_.emplace_back(N_, 0.0);	150✔
59	xI_.emplace_back(N_, 0.0);	150✔
60	xRW_.emplace_back(N_, 0.0);	150✔
61	Array& x = x_.back();
62	for (Size j = 0; j < N_; j++) {	450✔
63	//Assign X=lb+(ub-lb)*random
64	x[j] = lX_[j] + bounds[j] * sample[j];	300✔
65	}
66	//Evaluate point
67	values_.emplace_back(P.value(x), i);	150✔
68	}
69
70	//init intensity & randomWalk
71	intensity_->init(this);	1✔
72	randomWalk_->init(this);	1✔
73	}	2✔
74
75	EndCriteria::Type FireflyAlgorithm::minimize(Problem &P, const EndCriteria &endCriteria) {	1✔
76	QL_REQUIRE(!P.constraint().empty(), "Firefly Algorithm is a constrained optimizer");	1✔
77	EndCriteria::Type ecType = EndCriteria::None;
78	P.reset();
79	Size iteration = 0;
80	Size iterationStat = 0;
81	Size maxIteration = endCriteria.maxIterations();	1✔
82	Size maxIStationary = endCriteria.maxStationaryStateIterations();	1✔
83
84	startState(P, endCriteria);	1✔
85
86	bool isFA = Mfa_ > 0;	1✔
87	//Variables for DE
88	Array z(N_, 0.0);	1✔
89	Size indexR1, indexR2;
90	decltype(distribution_)::param_type nParam(0, N_ - 1);	1✔
91
92	//Set best value & position
93	Real bestValue = values_[0].first;	1✔
94	Size bestPosition = 0;
95	for (Size i = 1; i < M_; i++) {	150✔
96	if (values_[i].first < bestValue) {	149✔
97	bestPosition = i;
98	bestValue = values_[i].first;
99	}
100	}
101	Array bestX = x_[bestPosition];	1✔
102
103	//Run optimization
104	do {
105	iteration++;	1,011✔
106	iterationStat++;	1,011✔
107	//Check if stopping criteria is met
108	if (iteration > maxIteration \|\| iterationStat > maxIStationary)	1,011✔
109	break;
110
111	//Divide into two subpopulations
112	//First sort values
113	std::sort(values_.begin(), values_.end());	1,010✔
114
115	//Differential evolution
116	if(Mfa_ < M_){	1,010✔
117	Size indexBest = values_[0].second;	1,010✔
118	Array& xBest = x_[indexBest];
119	for (Size i = Mfa_; i < M_; i++) {	41,410✔
120	if (!isFA) {	40,400✔
121	//Pure DE requires random index
122	indexBest = distribution_(generator_);	×
123	xBest = x_[indexBest];	×
124	}
125	do {
126	indexR1 = distribution_(generator_);	40,677✔
127	} while(indexR1 == indexBest);	40,677✔
128	do {
129	indexR2 = distribution_(generator_);
130	} while(indexR2 == indexBest \|\| indexR2 == indexR1);	41,730✔
131
132	Size index = values_[i].second;	40,400✔
133	Array& x = x_[index];
134	Array& xR1 = x_[indexR1];
135	Array& xR2 = x_[indexR2];
136	Size rIndex = distribution_(generator_, nParam);	40,400✔
137	for (Size j = 0; j < N_; j++) {	121,200✔
138	if (j == rIndex \|\| rng_.nextReal() <= crossover_) {	80,800✔
139	//Change x[j] according to crossover
140	z[j] = xBest[j] + mutation_*(xR1[j] - xR2[j]);	60,612✔
141	} else {
142	z[j] = x[j];	20,188✔
143	}
144	//Enforce bounds on positions
145	if (z[j] < lX_[j]) {	80,800✔
146	z[j] = lX_[j];	6,947✔
147	}
148	else if (z[j] > uX_[j]) {	73,853✔
149	z[j] = uX_[j];	17,102✔
150	}
151	}
152	Real val = P.value(z);	40,400✔
153	if (val < values_[index].first) {	40,400✔
154	//Accept new point
155	x = z;	8,171✔
156	values_[index].first = val;	8,171✔
157	//mark best
158	if (val < bestValue) {	8,171✔
159	bestValue = val;
160	bestX = x;	2✔
161	iterationStat = 0;
162	}
163	}
164	}
165	}
166
167	//Firefly algorithm
168	if(isFA){	1,010✔
169	//According to the intensity, determine best global position
170	intensity_->findBrightest();	1,010✔
171
172	//Prepare random walk
173	randomWalk_->walk();	1,010✔
174
175	//Loop over particles
176	for (Size i = 0; i < Mfa_; i++) {	112,110✔
177	Size index = values_[i].second;	111,100✔
178	Array& x = x_[index];
179	Array& xI = xI_[index];
180	Array& xRW = xRW_[index];
181
182	//Loop over dimensions
183	for (Size j = 0; j < N_; j++) {	333,300✔
184	//Update position
185	z[j] = x[j] + xI[j] + xRW[j];	222,200✔
186	//Enforce bounds on positions
187	if (z[j] < lX_[j]) {	222,200✔
188	z[j] = lX_[j];	2,145✔
189	}
190	else if (z[j] > uX_[j]) {	220,055✔
191	z[j] = uX_[j];	10,455✔
192	}
193	}
194	Real val = P.value(z);	111,100✔
195	if(!std::isnan(val))	111,100✔
196	{
197	//Accept new point
198	x = z;	111,100✔
199	values_[index].first = val;	111,100✔
200	//mark best
201	if (val < bestValue) {	111,100✔
202	bestValue = val;
UNCOV 203	bestX = x;	×
204	iterationStat = 0;
205	}
206	}
207	}
208	}
209	} while (true);
210	if (iteration > maxIteration)	1✔
211	ecType = EndCriteria::MaxIterations;
212	else
213	ecType = EndCriteria::StationaryPoint;
214
215	//Set result to best point
216	P.setCurrentValue(bestX);	2✔
217	P.setFunctionValue(bestValue);
218	return ecType;	1✔
219	}
220
221	void FireflyAlgorithm::Intensity::findBrightest() {	1,010✔
222	//Brightest ignores all others
223	Array& xI = (xI_)[(values_)[0].second];	1,010✔
224	for (Size j = 0; j < N_; j++) {	3,030✔
225	xI[j] = 0.0;	2,020✔
226	}
227
228	for (Size i = 1; i < Mfa_; i++) {	111,100✔
229	//values_ is already sorted
230	Size index = (*values_)[i].second;	110,090✔
231	const Array& x = (*x_)[index];	110,090✔
232	Array& xI = (*xI_)[index];	110,090✔
233	for (Size j = 0; j < N_; j++) {	330,270✔
234	xI[j] = 0.0;	220,180✔
235	}
236	Real valueX = (*values_)[i].first;	110,090✔
237	for (Size k = 0; k < i - 1; k++){	6,054,950✔
238	const Array& y = (x_)[(values_)[k].second];	5,944,860✔
239	Real valueY = (*values_)[k].first;	5,944,860✔
240	Real intensity = intensityImpl(valueX, valueY, distance(x, y));	5,944,860✔
241	xI += intensity*(y - x);	11,889,720✔
242	}
243	}
244	}	1,010✔
245	}
246

lballabio / QuantLib / 13097754652

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous