/src/quantlib/ql/experimental/math/particleswarmoptimization.cpp

Source
/* -*- 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
<https://www.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/particleswarmoptimization.hpp>
#include <ql/math/randomnumbers/sobolrsg.hpp>
#include <cmath>
#include <utility>

using std::sqrt;

namespace QuantLib {
    ParticleSwarmOptimization::ParticleSwarmOptimization(Size M,
                                                         ext::shared_ptr<Topology> topology,
                                                         ext::shared_ptr<Inertia> inertia,
                                                         Real c1,
                                                         Real c2,
                                                         unsigned long seed)
    : M_(M), rng_(seed), topology_(std::move(topology)), inertia_(std::move(inertia)) {
        Real phi = c1 + c2;
        QL_ENSURE(phi*phi - 4 * phi != 0.0, "Invalid phi");
        c0_ = 2.0 / std::abs(2.0 - phi - sqrt(phi*phi - 4 * phi));
        c1_ = c0_*c1;
        c2_ = c0_*c2;
    }

    ParticleSwarmOptimization::ParticleSwarmOptimization(Size M,
                                                         ext::shared_ptr<Topology> topology,
                                                         ext::shared_ptr<Inertia> inertia,
                                                         Real omega,
                                                         Real c1,
                                                         Real c2,
                                                         unsigned long seed)
    : M_(M), c0_(omega), c1_(c1), c2_(c2), rng_(seed), topology_(std::move(topology)),
      inertia_(std::move(inertia)) {}

    void ParticleSwarmOptimization::startState(Problem &P, const EndCriteria &endCriteria) {
        QL_REQUIRE(topology_, "Invalid topology");
        QL_REQUIRE(inertia_, "Invalid inertia");
        N_ = P.currentValue().size();
        topology_->setSize(M_);
        inertia_->setSize(M_, N_, c0_, endCriteria);
        X_.reserve(M_);
        V_.reserve(M_);
        pBX_.reserve(M_);
        pBF_ = Array(M_);
        gBX_.reserve(M_);
        gBF_ = Array(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_ * 2);

        //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);
            Array& x = X_.back();
            V_.emplace_back(N_, 0.0);
            Array& v = V_.back();
            gBX_.emplace_back(N_, 0.0);
            for (Size j = 0; j < N_; j++) {
                //Assign X=lb+(ub-lb)*random
                x[j] = lX_[j] + bounds[j] * sample[2 * j];
                //Assign V=(ub-lb)*2*random-(ub-lb) -> between (lb-ub) and (ub-lb)
                v[j] = bounds[j] * (2.0*sample[2 * j + 1] - 1.0);
            }
            //Evaluate X and assign as personal best
            pBX_.push_back(X_.back());
            pBF_[i] = P.value(X_.back());
        }

        //init topology & inertia
        topology_->init(this);
        inertia_->init(this);
    }

    EndCriteria::Type ParticleSwarmOptimization::minimize(Problem &P, const EndCriteria &endCriteria) {
        QL_REQUIRE(!P.constraint().empty(), "PSO 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();
        Real bestValue = QL_MAX_REAL;
        Size bestPosition = 0;

        startState(P, endCriteria);
        //Set best value & position
        for (Size i = 0; i < M_; i++) {
            if (pBF_[i] < bestValue) {
                bestValue = pBF_[i];
                bestPosition = i;
            }
        }

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

            //According to the topology, determine best global position
            topology_->findSocialBest();

            //Call inertia to change internal state
            inertia_->setValues();

            //Loop over particles
            for (Size i = 0; i < M_; i++) {
                Array& x = X_[i];
                Array& pB = pBX_[i];
                const Array& gB = gBX_[i];
                Array& v = V_[i];

                //Loop over dimensions
                for (Size j = 0; j < N_; j++) {
                    //Update velocity
                    v[j] += c1_*rng_.nextReal()*(pB[j] - x[j]) + c2_*rng_.nextReal()*(gB[j] - x[j]);
                    //Update position
                    x[j] += v[j];
                    //Enforce bounds on positions
                    if (x[j] < lX_[j]) {
                        x[j] = lX_[j];
                        v[j] = 0.0;
                    }
                    else if (x[j] > uX_[j]) {
                        x[j] = uX_[j];
                        v[j] = 0.0;
                    }
                }
                //Evaluate x
                Real f = P.value(x);
                if (f < pBF_[i]) {
                    //Update personal best
                    pBF_[i] = f;
                    pB = x;
                    //Check stationary condition
                    if (f < bestValue) {
                        bestValue = f;
                        bestPosition = i;
                        iterationStat = 0;
                    }
                }
            }
        } while (true);
        if (iteration > maxIteration)
            ecType = EndCriteria::MaxIterations;
        else
            ecType = EndCriteria::StationaryPoint;

        //Set result to best point
        P.setCurrentValue(pBX_[bestPosition]);
        P.setFunctionValue(bestValue);
        return ecType;
    }

    void AdaptiveInertia::setValues() {
        Real currBest = (*pBF_)[0];
        for (Size i = 1; i < M_; i++) {
            if (currBest >(*pBF_)[i]) currBest = (*pBF_)[i];
        }
        if (started_) { //First iteration leaves inertia unchanged
            if (currBest < best_) {
                best_ = currBest;
                adaptiveCounter--;
            }
            else {
                adaptiveCounter++;
            }
            if (adaptiveCounter > sh_) {
                c0_ = std::max(minInertia_, std::min(maxInertia_, c0_*0.5));
            }
            else if (adaptiveCounter < sl_) {
                c0_ = std::max(minInertia_, std::min(maxInertia_, c0_*2.0));
            }
        }
        else {
            best_ = currBest;
            started_ = true;
        }
        for (Size i = 0; i < M_; i++) {
            (*V_)[i] *= c0_;
        }
    }

    void KNeighbors::findSocialBest() {
        for (Size i = 0; i < M_; i++) {
            Real bestF = (*pBF_)[i];
            Size bestX = 0;
            //Search K_ neightbors upwards
            Size upper = std::min(i + K_, M_);
            //Search K_ neighbors downwards
            Size lower = std::max(i, K_ + 1) - K_ - 1;
            for (Size j = lower; j < upper; j++) {
                if ((*pBF_)[j] < bestF) {
                    bestF = (*pBF_)[j];
                    bestX = j;
                }
            }
            if (i + K_ >= M_) { //loop around if i+K >= M_
                for (Size j = 0; j < i + K_ - M_; j++) {
                    if ((*pBF_)[j] < bestF) {
                        bestF = (*pBF_)[j];
                        bestX = j;
                    }
                }
            }
            else if (i < K_) {//loop around from above
                for (Size j = M_ - (K_ - i) - 1; j < M_; j++) {
                    if ((*pBF_)[j] < bestF) {
                        bestF = (*pBF_)[j];
                        bestX = j;
                    }
                }
            }
            (*gBX_)[i] = (*pBX_)[bestX];
            (*gBF_)[i] = bestF;
        }
    }

    ClubsTopology::ClubsTopology(Size defaultClubs,
                                 Size totalClubs,
                                 Size maxClubs,
                                 Size minClubs,
                                 Size resetIteration,
                                 unsigned long seed)
    : totalClubs_(totalClubs), maxClubs_(maxClubs), minClubs_(minClubs),
      defaultClubs_(defaultClubs), resetIteration_(resetIteration), bestByClub_(totalClubs, 0),
      worstByClub_(totalClubs, 0), generator_(seed), distribution_(1, totalClubs_) {
        QL_REQUIRE(totalClubs_ >= defaultClubs_,
            "Total number of clubs must be larger or equal than default clubs");
        QL_REQUIRE(defaultClubs_ >= minClubs_,
            "Number of default clubs must be larger or equal than minimum clubs");
        QL_REQUIRE(maxClubs_ >= defaultClubs_,
            "Number of maximum clubs must be larger or equal than default clubs");
        QL_REQUIRE(totalClubs_ >= maxClubs_,
            "Total number of clubs must be larger or equal than maximum clubs");
    }

    void ClubsTopology::setSize(Size M) {
        M_ = M;

        if (defaultClubs_ < totalClubs_) {
            clubs4particles_ = std::vector<std::vector<bool> >(M_, std::vector<bool>(totalClubs_, false));
            particles4clubs_ = std::vector<std::vector<bool> >(totalClubs_, std::vector<bool>(M_, false));
            //Assign particles to clubs randomly
            for (Size i = 0; i < M_; i++) {
                std::vector<bool> &clubSet = clubs4particles_[i];
                for (Size j = 0; j < defaultClubs_; j++) {
                    Size index = distribution_(generator_);
                    while (clubSet[index]) { index = distribution_(generator_); }
                    clubSet[index] = true;
                    particles4clubs_[index][i] = true;
                }
            }
        }
        else {
            //Since totalClubs_ == defaultClubs_, then just initialize to true
            clubs4particles_ = std::vector<std::vector<bool> >(M_, std::vector<bool>(totalClubs_, true));
            particles4clubs_ = std::vector<std::vector<bool> >(totalClubs_, std::vector<bool>(M_, true));
        }
    }

    void ClubsTopology::findSocialBest() {
        //Update iteration
        iteration_++;
        bool reset = false;
        if (iteration_ == resetIteration_) {
            iteration_ = 0;
            reset = true;
        }

        //Find best by current club
        for (Size i = 0; i < totalClubs_; i++) {
            Real bestByClub = QL_MAX_REAL;
            Real worstByClub = -QL_MAX_REAL;
            Size bestP = 0;
            Size worstP = 0;
            const std::vector<bool> &particlesSet = particles4clubs_[i];
            for (Size j = 0; j < M_; j++) {
                if (particlesSet[j]) {
                    if (bestByClub >(*pBF_)[j]) {
                        bestByClub = (*pBF_)[j];
                        bestP = j;
                    }
                    else if (worstByClub < (*pBF_)[j]) {
                        worstByClub = (*pBF_)[j];
                        worstP = j;
                    }
                }
            }
            bestByClub_[i] = bestP;
            worstByClub_[i] = worstP;
        }

        //Update clubs && global best
        for (Size i = 0; i < M_; i++) {
            std::vector<bool> &clubSet = clubs4particles_[i];
            bool best = true;
            bool worst = true;
            Size currentClubs = 0;
            for (Size j = 0; j < totalClubs_; j++) {
                if (clubSet[j]) {
                    //If still thought of the best, check if best in club j
                    if (best && i != bestByClub_[j]) best = false;
                    //If still thought of the worst, check if worst in club j
                    if (worst && i != worstByClub_[j]) worst = false;
                    //Update currentClubs
                    currentClubs++;
                }
            }
            //Update clubs
            if (best) {
                //Leave random club
                leaveRandomClub(i, currentClubs);
            }
            else if (worst) {
                //Join random club
                joinRandomClub(i, currentClubs);
            }
            else if (reset && currentClubs != defaultClubs_) {
                //If membership != defaultClubs_, then leave or join accordingly
                if (currentClubs < defaultClubs_) {
                    //Join random club
                    joinRandomClub(i, currentClubs);
                }
                else {
                    //Leave random club
                    leaveRandomClub(i, currentClubs);
                }
            }

            //Update global best
            Real bestNeighborF = QL_MAX_REAL;
            Size bestNeighborX = 0;
            for (Size j = 0; j < totalClubs_; j++) {
                if (clubSet[j] && bestNeighborF >(*pBF_)[bestByClub_[j]]) {
                    bestNeighborF = (*pBF_)[bestByClub_[j]];
                    bestNeighborX = j;
                }
            }
            (*gBX_)[i] = (*pBX_)[bestNeighborX];
            (*gBF_)[i] = bestNeighborF;
        }
    }

    void ClubsTopology::leaveRandomClub(Size particle, Size currentClubs) {
        Size randIndex = distribution_(generator_, param_type(1, currentClubs));
        Size index = 1;
        std::vector<bool> &clubSet = clubs4particles_[particle];
        for (Size j = 0; j < totalClubs_; j++) {
            if (clubSet[j]) {
                if (index == randIndex) {
                    clubSet[j] = false;
                    particles4clubs_[j][particle] = false;
                    break;
                }
                index++;
            }
        }
    }

    void ClubsTopology::joinRandomClub(Size particle, Size currentClubs) {
        Size randIndex = totalClubs_ == currentClubs ? 1 :
            distribution_(generator_, param_type(1, totalClubs_ - currentClubs));
        Size index = 1;
        std::vector<bool> &clubSet = clubs4particles_[particle];
        for (Size j = 0; j < totalClubs_; j++) {
            if (!clubSet[j]) {
                if (index == randIndex) {
                    clubSet[j] = true;
                    particles4clubs_[j][particle] = true;
                    break;
                }
                index++;
            }
        }
    }
}


Coverage Report

Created: 2026-02-03 07:02

Line	Count	Source
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		<https://www.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/particleswarmoptimization.hpp>
21		#include <ql/math/randomnumbers/sobolrsg.hpp>
22		#include <cmath>
23		#include <utility>
24
25		using std::sqrt;
26
27		namespace QuantLib {
28		ParticleSwarmOptimization::ParticleSwarmOptimization(Size M,
29		ext::shared_ptr<Topology> topology,
30		ext::shared_ptr<Inertia> inertia,
31		Real c1,
32		Real c2,
33		unsigned long seed)
34	0	: M_(M), rng_(seed), topology_(std::move(topology)), inertia_(std::move(inertia)) {
35	0	Real phi = c1 + c2;
36	0	QL_ENSURE(phiphi - 4 phi != 0.0, "Invalid phi");
37	0	c0_ = 2.0 / std::abs(2.0 - phi - sqrt(phiphi - 4 phi));
38	0	c1_ = c0_*c1;
39	0	c2_ = c0_*c2;
40	0	}
41
42		ParticleSwarmOptimization::ParticleSwarmOptimization(Size M,
43		ext::shared_ptr<Topology> topology,
44		ext::shared_ptr<Inertia> inertia,
45		Real omega,
46		Real c1,
47		Real c2,
48		unsigned long seed)
49	0	: M_(M), c0_(omega), c1_(c1), c2_(c2), rng_(seed), topology_(std::move(topology)),
50	0	inertia_(std::move(inertia)) {}
51
52	0	void ParticleSwarmOptimization::startState(Problem &P, const EndCriteria &endCriteria) {
53	0	QL_REQUIRE(topology_, "Invalid topology");
54	0	QL_REQUIRE(inertia_, "Invalid inertia");
55	0	N_ = P.currentValue().size();
56	0	topology_->setSize(M_);
57	0	inertia_->setSize(M_, N_, c0_, endCriteria);
58	0	X_.reserve(M_);
59	0	V_.reserve(M_);
60	0	pBX_.reserve(M_);
61	0	pBF_ = Array(M_);
62	0	gBX_.reserve(M_);
63	0	gBF_ = Array(M_);
64	0	uX_ = P.constraint().upperBound(P.currentValue());
65	0	lX_ = P.constraint().lowerBound(P.currentValue());
66	0	Array bounds = uX_ - lX_;
67
68		//Random initialization is done by Sobol sequence
69	0	SobolRsg sobol(N_ * 2);
70
71		//Prepare containers
72	0	for (Size i = 0; i < M_; i++) {
73	0	const SobolRsg::sample_type::value_type &sample = sobol.nextSequence().value;
74	0	X_.emplace_back(N_, 0.0);
75	0	Array& x = X_.back();
76	0	V_.emplace_back(N_, 0.0);
77	0	Array& v = V_.back();
78	0	gBX_.emplace_back(N_, 0.0);
79	0	for (Size j = 0; j < N_; j++) {
80		//Assign X=lb+(ub-lb)*random
81	0	x[j] = lX_[j] + bounds[j] * sample[2 * j];
82		//Assign V=(ub-lb)2random-(ub-lb) -> between (lb-ub) and (ub-lb)
83	0	v[j] = bounds[j] * (2.0sample[2 j + 1] - 1.0);
84	0	}
85		//Evaluate X and assign as personal best
86	0	pBX_.push_back(X_.back());
87	0	pBF_[i] = P.value(X_.back());
88	0	}
89
90		//init topology & inertia
91	0	topology_->init(this);
92	0	inertia_->init(this);
93	0	}
94
95	0	EndCriteria::Type ParticleSwarmOptimization::minimize(Problem &P, const EndCriteria &endCriteria) {
96	0	QL_REQUIRE(!P.constraint().empty(), "PSO is a constrained optimizer");
97
98	0	EndCriteria::Type ecType = EndCriteria::None;
99	0	P.reset();
100	0	Size iteration = 0;
101	0	Size iterationStat = 0;
102	0	Size maxIteration = endCriteria.maxIterations();
103	0	Size maxIStationary = endCriteria.maxStationaryStateIterations();
104	0	Real bestValue = QL_MAX_REAL;
105	0	Size bestPosition = 0;
106
107	0	startState(P, endCriteria);
108		//Set best value & position
109	0	for (Size i = 0; i < M_; i++) {
110	0	if (pBF_[i] < bestValue) {
111	0	bestValue = pBF_[i];
112	0	bestPosition = i;
113	0	}
114	0	}
115
116		//Run optimization
117	0	do {
118	0	iteration++;
119	0	iterationStat++;
120		//Check if stopping criteria is met
121	0	if (iteration > maxIteration \|\| iterationStat > maxIStationary)
122	0	break;
123
124		//According to the topology, determine best global position
125	0	topology_->findSocialBest();
126
127		//Call inertia to change internal state
128	0	inertia_->setValues();
129
130		//Loop over particles
131	0	for (Size i = 0; i < M_; i++) {
132	0	Array& x = X_[i];
133	0	Array& pB = pBX_[i];
134	0	const Array& gB = gBX_[i];
135	0	Array& v = V_[i];
136
137		//Loop over dimensions
138	0	for (Size j = 0; j < N_; j++) {
139		//Update velocity
140	0	v[j] += c1_rng_.nextReal()(pB[j] - x[j]) + c2_rng_.nextReal()(gB[j] - x[j]);
141		//Update position
142	0	x[j] += v[j];
143		//Enforce bounds on positions
144	0	if (x[j] < lX_[j]) {
145	0	x[j] = lX_[j];
146	0	v[j] = 0.0;
147	0	}
148	0	else if (x[j] > uX_[j]) {
149	0	x[j] = uX_[j];
150	0	v[j] = 0.0;
151	0	}
152	0	}
153		//Evaluate x
154	0	Real f = P.value(x);
155	0	if (f < pBF_[i]) {
156		//Update personal best
157	0	pBF_[i] = f;
158	0	pB = x;
159		//Check stationary condition
160	0	if (f < bestValue) {
161	0	bestValue = f;
162	0	bestPosition = i;
163	0	iterationStat = 0;
164	0	}
165	0	}
166	0	}
167	0	} while (true);
168	0	if (iteration > maxIteration)
169	0	ecType = EndCriteria::MaxIterations;
170	0	else
171	0	ecType = EndCriteria::StationaryPoint;
172
173		//Set result to best point
174	0	P.setCurrentValue(pBX_[bestPosition]);
175	0	P.setFunctionValue(bestValue);
176	0	return ecType;
177	0	}
178
179	0	void AdaptiveInertia::setValues() {
180	0	Real currBest = (*pBF_)[0];
181	0	for (Size i = 1; i < M_; i++) {
182	0	if (currBest >(pBF_)[i]) currBest = (pBF_)[i];
183	0	}
184	0	if (started_) { //First iteration leaves inertia unchanged
185	0	if (currBest < best_) {
186	0	best_ = currBest;
187	0	adaptiveCounter--;
188	0	}
189	0	else {
190	0	adaptiveCounter++;
191	0	}
192	0	if (adaptiveCounter > sh_) {
193	0	c0_ = std::max(minInertia_, std::min(maxInertia_, c0_*0.5));
194	0	}
195	0	else if (adaptiveCounter < sl_) {
196	0	c0_ = std::max(minInertia_, std::min(maxInertia_, c0_*2.0));
197	0	}
198	0	}
199	0	else {
200	0	best_ = currBest;
201	0	started_ = true;
202	0	}
203	0	for (Size i = 0; i < M_; i++) {
204	0	(V_)[i] = c0_;
205	0	}
206	0	}
207
208	0	void KNeighbors::findSocialBest() {
209	0	for (Size i = 0; i < M_; i++) {
210	0	Real bestF = (*pBF_)[i];
211	0	Size bestX = 0;
212		//Search K_ neightbors upwards
213	0	Size upper = std::min(i + K_, M_);
214		//Search K_ neighbors downwards
215	0	Size lower = std::max(i, K_ + 1) - K_ - 1;
216	0	for (Size j = lower; j < upper; j++) {
217	0	if ((*pBF_)[j] < bestF) {
218	0	bestF = (*pBF_)[j];
219	0	bestX = j;
220	0	}
221	0	}
222	0	if (i + K_ >= M_) { //loop around if i+K >= M_
223	0	for (Size j = 0; j < i + K_ - M_; j++) {
224	0	if ((*pBF_)[j] < bestF) {
225	0	bestF = (*pBF_)[j];
226	0	bestX = j;
227	0	}
228	0	}
229	0	}
230	0	else if (i < K_) {//loop around from above
231	0	for (Size j = M_ - (K_ - i) - 1; j < M_; j++) {
232	0	if ((*pBF_)[j] < bestF) {
233	0	bestF = (*pBF_)[j];
234	0	bestX = j;
235	0	}
236	0	}
237	0	}
238	0	(gBX_)[i] = (pBX_)[bestX];
239	0	(*gBF_)[i] = bestF;
240	0	}
241	0	}
242
243		ClubsTopology::ClubsTopology(Size defaultClubs,
244		Size totalClubs,
245		Size maxClubs,
246		Size minClubs,
247		Size resetIteration,
248		unsigned long seed)
249	0	: totalClubs_(totalClubs), maxClubs_(maxClubs), minClubs_(minClubs),
250	0	defaultClubs_(defaultClubs), resetIteration_(resetIteration), bestByClub_(totalClubs, 0),
251	0	worstByClub_(totalClubs, 0), generator_(seed), distribution_(1, totalClubs_) {
252	0	QL_REQUIRE(totalClubs_ >= defaultClubs_,
253	0	"Total number of clubs must be larger or equal than default clubs");
254	0	QL_REQUIRE(defaultClubs_ >= minClubs_,
255	0	"Number of default clubs must be larger or equal than minimum clubs");
256	0	QL_REQUIRE(maxClubs_ >= defaultClubs_,
257	0	"Number of maximum clubs must be larger or equal than default clubs");
258	0	QL_REQUIRE(totalClubs_ >= maxClubs_,
259	0	"Total number of clubs must be larger or equal than maximum clubs");
260	0	}
261
262	0	void ClubsTopology::setSize(Size M) {
263	0	M_ = M;
264
265	0	if (defaultClubs_ < totalClubs_) {
266	0	clubs4particles_ = std::vector<std::vector<bool> >(M_, std::vector<bool>(totalClubs_, false));
267	0	particles4clubs_ = std::vector<std::vector<bool> >(totalClubs_, std::vector<bool>(M_, false));
268		//Assign particles to clubs randomly
269	0	for (Size i = 0; i < M_; i++) {
270	0	std::vector<bool> &clubSet = clubs4particles_[i];
271	0	for (Size j = 0; j < defaultClubs_; j++) {
272	0	Size index = distribution_(generator_);
273	0	while (clubSet[index]) { index = distribution_(generator_); }
274	0	clubSet[index] = true;
275	0	particles4clubs_[index][i] = true;
276	0	}
277	0	}
278	0	}
279	0	else {
280		//Since totalClubs_ == defaultClubs_, then just initialize to true
281	0	clubs4particles_ = std::vector<std::vector<bool> >(M_, std::vector<bool>(totalClubs_, true));
282	0	particles4clubs_ = std::vector<std::vector<bool> >(totalClubs_, std::vector<bool>(M_, true));
283	0	}
284	0	}
285
286	0	void ClubsTopology::findSocialBest() {
287		//Update iteration
288	0	iteration_++;
289	0	bool reset = false;
290	0	if (iteration_ == resetIteration_) {
291	0	iteration_ = 0;
292	0	reset = true;
293	0	}
294
295		//Find best by current club
296	0	for (Size i = 0; i < totalClubs_; i++) {
297	0	Real bestByClub = QL_MAX_REAL;
298	0	Real worstByClub = -QL_MAX_REAL;
299	0	Size bestP = 0;
300	0	Size worstP = 0;
301	0	const std::vector<bool> &particlesSet = particles4clubs_[i];
302	0	for (Size j = 0; j < M_; j++) {
303	0	if (particlesSet[j]) {
304	0	if (bestByClub >(*pBF_)[j]) {
305	0	bestByClub = (*pBF_)[j];
306	0	bestP = j;
307	0	}
308	0	else if (worstByClub < (*pBF_)[j]) {
309	0	worstByClub = (*pBF_)[j];
310	0	worstP = j;
311	0	}
312	0	}
313	0	}
314	0	bestByClub_[i] = bestP;
315	0	worstByClub_[i] = worstP;
316	0	}
317
318		//Update clubs && global best
319	0	for (Size i = 0; i < M_; i++) {
320	0	std::vector<bool> &clubSet = clubs4particles_[i];
321	0	bool best = true;
322	0	bool worst = true;
323	0	Size currentClubs = 0;
324	0	for (Size j = 0; j < totalClubs_; j++) {
325	0	if (clubSet[j]) {
326		//If still thought of the best, check if best in club j
327	0	if (best && i != bestByClub_[j]) best = false;
328		//If still thought of the worst, check if worst in club j
329	0	if (worst && i != worstByClub_[j]) worst = false;
330		//Update currentClubs
331	0	currentClubs++;
332	0	}
333	0	}
334		//Update clubs
335	0	if (best) {
336		//Leave random club
337	0	leaveRandomClub(i, currentClubs);
338	0	}
339	0	else if (worst) {
340		//Join random club
341	0	joinRandomClub(i, currentClubs);
342	0	}
343	0	else if (reset && currentClubs != defaultClubs_) {
344		//If membership != defaultClubs_, then leave or join accordingly
345	0	if (currentClubs < defaultClubs_) {
346		//Join random club
347	0	joinRandomClub(i, currentClubs);
348	0	}
349	0	else {
350		//Leave random club
351	0	leaveRandomClub(i, currentClubs);
352	0	}
353	0	}
354
355		//Update global best
356	0	Real bestNeighborF = QL_MAX_REAL;
357	0	Size bestNeighborX = 0;
358	0	for (Size j = 0; j < totalClubs_; j++) {
359	0	if (clubSet[j] && bestNeighborF >(*pBF_)[bestByClub_[j]]) {
360	0	bestNeighborF = (*pBF_)[bestByClub_[j]];
361	0	bestNeighborX = j;
362	0	}
363	0	}
364	0	(gBX_)[i] = (pBX_)[bestNeighborX];
365	0	(*gBF_)[i] = bestNeighborF;
366	0	}
367	0	}
368
369	0	void ClubsTopology::leaveRandomClub(Size particle, Size currentClubs) {
370	0	Size randIndex = distribution_(generator_, param_type(1, currentClubs));
371	0	Size index = 1;
372	0	std::vector<bool> &clubSet = clubs4particles_[particle];
373	0	for (Size j = 0; j < totalClubs_; j++) {
374	0	if (clubSet[j]) {
375	0	if (index == randIndex) {
376	0	clubSet[j] = false;
377	0	particles4clubs_[j][particle] = false;
378	0	break;
379	0	}
380	0	index++;
381	0	}
382	0	}
383	0	}
384
385	0	void ClubsTopology::joinRandomClub(Size particle, Size currentClubs) {
386	0	Size randIndex = totalClubs_ == currentClubs ? 1 :
387	0	distribution_(generator_, param_type(1, totalClubs_ - currentClubs));
388	0	Size index = 1;
389	0	std::vector<bool> &clubSet = clubs4particles_[particle];
390	0	for (Size j = 0; j < totalClubs_; j++) {
391	0	if (!clubSet[j]) {
392	0	if (index == randIndex) {
393	0	clubSet[j] = true;
394	0	particles4clubs_[j][particle] = true;
395	0	break;
396	0	}
397	0	index++;
398	0	}
399	0	}
400	0	}
401		}
402