MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
PSO.cpp
Go to the documentation of this file.
1#include "PSO.h"
2
3#include <cmath>
4
5// ***************
7// ***************
8 fConstriction = Get<double>(fitMan->raw()["General"]["PSO"]["Constriction"], __FILE__, __LINE__);
9 fInertia = Get<double>(fitMan->raw()["General"]["PSO"]["Inertia"], __FILE__, __LINE__) * fConstriction;
10 fOne = Get<double>(fitMan->raw()["General"]["PSO"]["One"], __FILE__, __LINE__) * fConstriction;
11 fTwo = Get<double>(fitMan->raw()["General"]["PSO"]["Two"], __FILE__, __LINE__) * fConstriction;
12 fParticles = Get<int>(fitMan->raw()["General"]["PSO"]["Particles"], __FILE__, __LINE__);
13 fIterations = Get<int>(fitMan->raw()["General"]["PSO"]["Iterations"], __FILE__, __LINE__);
14 fConvergence = Get<double>(fitMan->raw()["General"]["PSO"]["Convergence"], __FILE__, __LINE__);
15
16 fDim = 0;
17
19 {
20 fDim = Get<int>(fitMan->raw()["General"]["PSO"]["TestLikelihoodDim"], __FILE__, __LINE__);
21 }
22}
23
24// ***************
26// ***************
27 PrepareFit();
28
29 // Remove obsolete memory and make other checks before fit starts
31
33 outTree->Branch("nParts", &fParticles, "nParts/I");
34 for(int i = 0; i < fDim; ++i){
35 par = new double[fParticles];
36 paramlist.push_back(par);
37 outTree->Branch(Form("Parameter_%d", i), paramlist[i], Form("Parameter_%d[nParts]/D",i));
38 }
39// vel = new double[fParticles];
40 outTree->Branch("vel", vel, "vel[nParts]/D");
41 }
42
43 init();
44 run();
46 return;
47}
48
49// *************************
50void PSO::init(){
51// *************************
53
54 //KS: For none PCA this will be eqaul to normal parameters
55 //const int NparsPSOFull = NPars;
56 //const int NparsPSO = NParsPCA;
57
58 MACH3LOG_INFO("Preparing PSO");
59 // Initialise bounds on parameters
61 for (int i = 0; i < fDim; i++){
62 // Test function ranges
63 ranges_min.push_back(-5);
64 ranges_max.push_back(5);
65 fixed.push_back(0);
66 }
67 }
68 else{
69 for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it){
70 if(!(*it)->IsPCA())
71 {
72 fDim += (*it)->GetNumParams();
73 for(int i = 0; i < (*it)->GetNumParams(); ++i)
74 {
75 double curr = (*it)->GetParInit(i);
76 double lim = 10.0*(*it)->GetDiagonalError(i);
77 double low = (*it)->GetLowerBound(i);
78 double high = (*it)->GetUpperBound(i);
79 if(low > curr - lim) ranges_min.push_back(low);
80 else ranges_min.push_back(curr - lim);
81 if(high < curr + lim) ranges_min.push_back(high);
82 else ranges_min.push_back(curr + lim);
83 prior.push_back(curr);
84
85 if((*it)->IsParameterFixed(i)){
86 fixed.push_back(1);
87 }
88 else{
89 fixed.push_back(0);
90 }
91 }
92 }
93 else
94 {
95 fDim += (*it)->GetNParameters();
96 for(int i = 0; i < (*it)->GetNParameters(); ++i)
97 {
98 ranges_min.push_back(-100.0);
99 ranges_max.push_back(100.0);
100 prior.push_back((*it)->GetParInit(i));
101 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i)){
102 fixed.push_back(1);
103 }
104 else{
105 fixed.push_back(0);
106 }
107 }
108 }
109 }
110 }
111
112 MACH3LOG_INFO("Printing Minimums and Maximums of Variables to be minimized");
113 for (int i = 0; i < fDim; i++){
114 MACH3LOG_INFO("Variable {} : {:.2f}, {:.2f}", i, ranges_min[i], ranges_max[i]);
115 }
116
117 // Initialise particle positions
118 for (int i = 0; i < fParticles; ++i){
119 std::vector<double> init_position;
120 std::vector<double> init_velocity;
121
122 //Initialising in +/- 5sigma of prior value from BANFF interface
123 for (int j=0; j<fDim; ++j){
124 if(fixed[j]){
125 init_position.push_back(prior[j]);
126 init_velocity.push_back(0.0);
127 }
128 else{
129 double dist = fabs(ranges_max[j]-ranges_min[j]);
130 //Initialise to random position uniform in space
131 init_position.push_back(ranges_min[j] + random->Rndm()*dist);
132 //Initialise velocity to random position uniform in space
133 init_velocity.push_back((2.0*random->Rndm()-1.0));//*dist);
134 }
135 }
136
137 particle* new_particle = new particle(init_position,init_velocity);
138 new_particle->set_personal_best_position(init_position);
139 double new_value = CalcChi(init_position);
140 new_particle->set_personal_best_value(new_value);
141 new_particle->set_value(new_value);
142 system.push_back(new_particle);
143 if(new_value < fBestValue){
144 fBestValue = new_value;
145 set_best_particle(new_particle);
146 }
147 }
148}
149
150// *************************
151std::vector<std::vector<double> > PSO::bisection(const std::vector<double>& position, const double minimum,
152 const double range, const double precision) {
153// *************************
154 std::vector<std::vector<double>> uncertainties_list;
155 for (unsigned int i = 0; i< position.size(); ++i) {
156 MACH3LOG_INFO("{}", i);
157 //std::vector<double> uncertainties;
158 std::vector<double> new_position = position; new_position[i] = position[i]-range;
159 double val_1 = CalcChi(new_position)-minimum-1.0;
160 while (val_1*-1.0 > 0.0){
161 new_position[i] -= range;
162 val_1 = CalcChi(new_position)-minimum-1.0;
163 }
164 std::vector<double> bisect_position = position; bisect_position[i] = bisect_position[i] - (position[i]-new_position[i])/2;
165 std::vector<std::vector<double>> position_list{new_position,bisect_position,position};
166 double val_2 = CalcChi(bisect_position)-minimum-1.0;
167 std::vector<double> value_list{val_1,val_2, -1.0};
168 double res = 1.0;
169 while (res > precision){
170 if (value_list[0] * value_list[1] < 0){
171 std::vector<double> new_bisect_position = position_list[0];
172 new_bisect_position[i] =new_bisect_position[i]+ (position_list[1][i]-position_list[0][i])/2;
173 double new_val = CalcChi(new_bisect_position)-minimum-1.0;
174 position_list[2] = position_list[1];
175 value_list[2] = value_list[1];
176 value_list[1] = new_val;
177 position_list[1] = new_bisect_position;
178 res = std::abs(position[2]-position[0]);
179 }
180 else{
181 std::vector<double> new_bisect_position = position_list[1];
182 new_bisect_position[i] += (position_list[2][i]-position_list[1][i])/2;
183 double new_val = CalcChi(new_bisect_position)-minimum-1.0;
184 position_list[0] = position_list[1];
185 value_list[0] = value_list[1];
186 value_list[1] = new_val;
187 position_list[1] = new_bisect_position;
188 res = std::abs(position_list[2][i]-position_list[1][i]);
189 }
190 }
191 //do the same thing for position uncertainty
192 std::vector<double> new_position_p = position; new_position_p[i] = position[i]+range;
193 double val_1_p = CalcChi(new_position_p)-minimum-1.0;
194 while (val_1_p * -1.0 > 0.0){
195 new_position_p[i] += range;
196 val_1_p = CalcChi(new_position_p)-minimum-1.0;
197 }
198 std::vector<double> bisect_position_p = position; bisect_position_p[i] = bisect_position_p[i] += (new_position_p[i]-position[i])/2;
199 std::vector<std::vector<double>> position_list_p{position,bisect_position_p,new_position_p};
200 double val_2_p = CalcChi(bisect_position_p)-minimum-1.0;
201 std::vector<double> value_list_p{-1.0,val_2_p, val_1_p};
202 double res_p = 1.0;
203 while (res_p > precision){
204 if (value_list_p[0] * value_list_p[1] < 0){
205 std::vector<double> new_bisect_position_p = position_list_p[0];new_bisect_position_p[i] += (position_list_p[1][i]-position_list_p[0][i])/2;
206 double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
207 position_list_p[2] = position_list_p[1];
208 value_list_p[2] = value_list_p[1];
209 value_list_p[1] = new_val_p;
210 position_list_p[1] = new_bisect_position_p;
211 res = std::abs(position[2]-position[0]);
212 res_p = std::abs(position_list_p[1][i]-position_list_p[0][i]);
213 MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]);
214
215 }
216 else{
217 std::vector<double> new_bisect_position_p = position_list_p[1];new_bisect_position_p[i] += (position_list_p[2][i]-position_list_p[1][i])/2;
218 double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
219 position_list_p[0] = position_list_p[1];
220 value_list_p[0] = value_list_p[1];
221 value_list_p[1] = new_val_p;
222 position_list_p[1] = new_bisect_position_p;
223 res_p = std::abs(position_list_p[2][i]-position_list_p[1][i]);
224 MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]);
225 }
226 }
227 uncertainties_list.push_back({std::abs(position[i]-position_list[1][i]),std::abs(position[i]-position_list_p[1][i])});
228 MACH3LOG_INFO("Uncertainty finished for d = {}", i);
229 MACH3LOG_INFO("LLR values for ± positive and negative uncertainties are {:<10.2f} and {:<10.2f}",
230 CalcChi(position_list[1]), CalcChi(position_list_p[1]));
231 }
232 return uncertainties_list;
233}
234
235// *************************
236std::vector<std::vector<double>> PSO::calc_uncertainty(const std::vector<double>& position, const double minimum) {
237// *************************
238 std::vector<double> pos_uncertainty(position.size());
239 std::vector<double> neg_uncertainty(position.size());
240 constexpr int num = 200;
241 std::vector<double> pos = position;
242 for (unsigned int i = 0; i < position.size(); ++i) {
243 double curr_ival = pos[i];
244
245 double neg_stop = position[i] - 5e-2;
246 double pos_stop = position[i] + 5e-2;
247 double start = position[i];
248 std::vector<double> x(num);
249 std::vector<double> y(num);
250 double StepPoint = (start-neg_stop) / (num - 1);
251 double value = start;
252 for (int j = 0; j < num; ++j) {
253 pos[i] = value;
254 double LLR = CalcChi(position) - minimum - 1.0;
255 x[j] = value;
256 y[j] = LLR;
257 value -= StepPoint;
258 }
259 pos[i] = curr_ival;
260
261 int closest_index = 0;
262 double closest_value = std::abs(y[0]); // Initialize with the first element
263 for (unsigned int ii = 1; ii < y.size(); ++ii) {
264 double abs_y = std::abs(y[ii]);
265 if (abs_y < closest_value) {
266 closest_index = ii;
267 closest_value = abs_y;
268 }
269 }
270 neg_uncertainty[i] = x[closest_index];
271 MACH3LOG_INFO("Neg");
272 x.assign(num, 0);
273 y.assign(num, 0);
274 StepPoint = (pos_stop-start) / (num - 1);
275 value = start;
276 for (int j = 0; j < num; ++j) {
277 pos[i] = value;
278 double LLR = CalcChi(position) - minimum - 1.0;
279 x[j] = value;
280 y[j] = LLR;
281 value += StepPoint;
282 }
283 pos[i] = curr_ival;
284 closest_index = 0;
285 closest_value = std::abs(y[0]); // Initialize with the first element
286 for (unsigned int ii = 1; ii < y.size(); ++ii) {
287 double abs_y = std::abs(y[ii]);
288 if (abs_y < closest_value) {
289 closest_index = ii;
290 closest_value = abs_y;
291 }
292 }
293 pos_uncertainty[i] = x[closest_index];
294 }
295 std::vector<std::vector<double>> res{neg_uncertainty,pos_uncertainty};
296 return res;
297}
298
299// *************************
300void PSO::uncertainty_check(const std::vector<double>& previous_pos){
301// *************************
302 std::vector<std::vector<double >> x_list;
303 std::vector<std::vector<double >> y_list;
304 std::vector<double> position = previous_pos;
305 constexpr int num = 5000;
306 for (unsigned int i = 0;i<previous_pos.size();++i){
307 double curr_ival = position[i];
308 double start = previous_pos[i] - 1e-1;
309 double stop = previous_pos[i] + 1e-1;
310 std::vector<double> x(num);
311 std::vector<double> y(num);
312 double StepPoint = (stop - start) / (num - 1);
313 double value = start;
314 MACH3LOG_TRACE("result for fDim: {}", 1);
315 for (int j = 0;j < num; ++j) {
316 position[i] = value;
317 double LLR = CalcChi(position);
318 x[j] = value;
319 y[j] = LLR;
320 value += StepPoint;
321 }
322 position[i] = curr_ival;
323 MACH3LOG_INFO("");
324 MACH3LOG_INFO("For fDim{} x list is", i+1);
325 for (unsigned int k = 0;k<x.size(); ++k){
326 MACH3LOG_INFO(" {}", x[k]);
327 }
328 MACH3LOG_INFO("");
329 MACH3LOG_INFO("For fDim{} y list is", i+1);
330 for (unsigned int k = 0; k < x.size(); ++k){
331 MACH3LOG_INFO(" {}", y[k]);
332 }
333 MACH3LOG_INFO("");
334 }
335}
336
337// *************************
339// *************************
340 std::vector<double> total_pos(fDim,0.0);
341
342 for (int i = 0; i < fParticles; ++i) {
343 std::vector<double> part1 = vector_multiply(system[i]->get_velocity(), fInertia);
344 std::vector<double> part2 = vector_multiply(vector_subtract(system[i]->get_personal_best_position(), system[i]->get_position()), (fOne * random->Rndm()));
345 std::vector<double> part3 = vector_multiply(vector_subtract(get_best_particle()->get_personal_best_position(), system[i]->get_position()),(fTwo * random->Rndm()));
346 std::vector<double> new_velocity = three_vector_addition(part1, part2, part3);
347 std::vector<double> new_pos = vector_add(system[i]->get_position(), new_velocity);
348 transform(total_pos.begin(), total_pos.end(), new_pos.begin(), total_pos.begin(),[](double x, double y) {return x+y;});
349
350 for (int j = 0; j < fDim; ++j) {
351 // Check if out of bounds and reflect if so
352 if(ranges_min[j] > new_pos[j]){
353 new_pos[j] = ranges_min[j];
354 }
355 else if(new_pos[j] > ranges_max[j]) {
356 new_pos[j] = ranges_max[j];
357 }
358 // If parameter fixed don't update it
359 if(fixed[j]) new_pos[j] = system[i]->get_position()[j];
360 }
361
362 if(fTestLikelihood){
363 double velo = 0.0;
364 for (int j = 0; j < fDim; ++j) {
365 paramlist[j][i] = new_pos[j];
366 velo += new_velocity[j]*new_velocity[j];
367 }
368 vel[i] = sqrt(velo);
369 }
370
371 system[i]->set_velocity(new_velocity);
372 system[i]->set_position(new_pos);
373 double new_value = CalcChi(new_pos);
374 if(new_value <= system[i]->get_personal_best_value()) {
375 system[i]->set_personal_best_value(new_value);
376 system[i]->set_personal_best_position(new_pos);
377 if(new_value < fBestValue){
378 fBestValue = new_value;
380 }
381 }
382 }
383
384 std::vector<double> best_pos = get_best_particle()->get_personal_best_position();
385 std::vector<double> result(best_pos.size(), 0.0);
386 transform(total_pos.begin(), total_pos.end(), total_pos.begin(), [=](double val){return val/fParticles;});
387 transform(total_pos.begin(),total_pos.end(),best_pos.begin(),result.begin(),[](double x, double y) {return x-y;});
388
389 double mean_dist_sq = 0;
390 for (int i = 0; i<fDim; i++){
391 mean_dist_sq += result[i]*result[i];
392 }
393
394 return mean_dist_sq;
395}
396
397// *************************
398void PSO::run() {
399// *************************
400 double mean_dist_sq = 0;
401
402 int iter = 0;
403 for(int i = 0; i < fIterations; ++i, ++iter){
404 mean_dist_sq = swarmIterate();
405 //double meanVel = std::accumulate(vel, vel + fParticles, 0) / fParticles;
406
407 // Weight inertia randomly but scaled by total distance of swarm from global minimum - proxy for total velocity
408 // fWeight = ((random->Rndm()+1.0)*0.5)*(10.0/meanVel);
409
411
412 outTree->Fill();
413 // Auto save the output
414 if (step % auto_save == 0) outTree->AutoSave();
415 step++;
416 accCount++;
417
418 if (i%100 == 0){
419 MACH3LOG_INFO("Mean Dist Sq = {:.2f}", mean_dist_sq);
420 MACH3LOG_INFO("Current LLR = {:.2f}", fBestValue);
421 MACH3LOG_INFO("Position:");
422 for (int j = 0; j < fDim; ++j){
423 MACH3LOG_INFO(" Dim {} = {:<10.2f}", j, get_best_particle()->get_personal_best_position()[j]);
424 }
425 }
426 if(fConvergence > 0.0){
427 if(mean_dist_sq < fConvergence){
428 break;
429 }
430 }
431 }
432 MACH3LOG_INFO("Finished after {} runs out of {}", iter, fIterations);
433 MACH3LOG_INFO("Mean Dist: {:.2f}", mean_dist_sq);
434 MACH3LOG_INFO("Best LLR: {:.2f}", get_best_particle()->get_personal_best_value());
435 uncertainties = bisection(get_best_particle()->get_personal_best_position(),get_best_particle()->get_personal_best_value(),0.5,0.005);
436 MACH3LOG_INFO("Position for Global Minimum = ");
437 for (int i = 0; i< fDim; ++i){
438 MACH3LOG_INFO(" Dim {} = {:<10.2f} +{:.2f}, -{:.2f}", i, get_best_particle()->get_personal_best_position()[i], uncertainties[i][1], uncertainties[i][0]);
439 }
440}
441
442// *************************
444// *************************
445 outputFile->cd();
446
447 TVectorD* PSOParValue = new TVectorD(fDim);
448 TVectorD* PSOParError = new TVectorD(fDim);
449
450 for(int i = 0; i < fDim; ++i)
451 {
452 (*PSOParValue)(i) = 0;
453 (*PSOParError)(i) = 0;
454 }
455
456 std::vector<double> minimum = get_best_particle()->get_personal_best_position();
457
458 int ParCounter = 0;
459
460 if(fTestLikelihood){
461 for(int i = 0; i < fDim; ++i){
462 (*PSOParValue)(i) = minimum[i];
463 (*PSOParError)(i) = (uncertainties[i][0]+uncertainties[i][1])/2.0;
464 }
465 }
466 else{
467 for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
468 {
469 if(!(*it)->IsPCA())
470 {
471 for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
472 {
473 double ParVal = minimum[ParCounter];
474 //KS: Basically apply mirroring for parameters out of bounds
475 (*PSOParValue)(ParCounter) = ParVal;
476 (*PSOParError)(ParCounter) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
477 //KS: For fixed params HESS will not calcuate error so we need to pass prior error
478 if((*it)->IsParameterFixed(i))
479 {
480 (*PSOParError)(ParCounter) = (*it)->GetDiagonalError(i);
481 }
482 }
483 }
484 else
485 {
486 //KS: We need to convert parameters from PCA to normal base
487 TVectorD ParVals((*it)->GetNumParams());
488 TVectorD ParVals_PCA((*it)->GetNParameters());
489
490 TVectorD ErrorVals((*it)->GetNumParams());
491 TVectorD ErrorVals_PCA((*it)->GetNParameters());
492
493 //First save them
494 //KS: This code is super convoluted as MaCh3 can store separate matrices while PSO has one matrix. In future this will be simplified, keep it like this for now.
495 const int StartVal = ParCounter;
496 for(int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
497 {
498 ParVals_PCA(i) = minimum[ParCounter];
499 ErrorVals_PCA(i) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
500 }
501 ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
502 ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
503
504 ParCounter = StartVal;
505 //KS: Now after going from PCA to normal let';s save it
506 for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
507 {
508 (*PSOParValue)(ParCounter) = ParVals(i);
509 (*PSOParError)(ParCounter) = std::fabs(ErrorVals(i));
510 //int ParCounterMatrix = StartVal;
511 //If fixed take prior
512 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
513 {
514 (*PSOParError)(ParCounter) = (*it)->GetDiagonalError(i);
515 }
516 }
517 }
518 }
519 }
520
521 PSOParValue->Write("PSOParValue");
522 PSOParError->Write("PSOParError");
523 delete PSOParValue;
524 delete PSOParError;
525 // Save all the output
526 SaveOutput();
527}
528
529// *******************
530double PSO::CalcChi2(const double* x) {
531// *******************
532 if(fTestLikelihood) {
533 return rastriginFunc(x);
534 } else {
535 return LikelihoodFit::CalcChi2(x);
536 }
537}
538
539// *************************
540double PSO::rastriginFunc(const double* x) {
541// *************************
542 stepClock->Start();
543
544 //Search range: [-5.12, 5.12]
545 constexpr double A = 10.0;
546 double sum = 0.0;
547 for (int i = 0; i < fDim; ++i) {
548 sum += x[i] * x[i] - A * cos(2.0 * 3.14 * x[i]);
549 }
550 double llh = A * fDim + sum;
551
552 accProb = 1;
553
554 stepClock->Stop();
555 stepTime = stepClock->RealTime();
556
557 return llh;
558}
MaCh3Plotting::PlottingManager * man
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_TRACE
Definition: MaCh3Logger.h:21
std::unique_ptr< TRandom3 > random
Random number.
Definition: FitterBase.h:128
int accCount
counts accepted steps
Definition: FitterBase.h:103
void SaveOutput()
Save output and close files.
Definition: FitterBase.cpp:221
TFile * outputFile
Output.
Definition: FitterBase.h:131
manager * fitMan
The manager.
Definition: FitterBase.h:92
unsigned int step
current state
Definition: FitterBase.h:95
double accProb
current acceptance prob
Definition: FitterBase.h:101
double stepTime
Time of single step.
Definition: FitterBase.h:125
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
Definition: FitterBase.h:123
double logLCurr
current likelihood
Definition: FitterBase.h:97
int auto_save
auto save every N steps
Definition: FitterBase.h:139
bool fTestLikelihood
Necessary for some fitting algorithms like PSO.
Definition: FitterBase.h:142
TTree * outTree
Output tree with posteriors.
Definition: FitterBase.h:137
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
Definition: FitterBase.cpp:213
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:118
Implementation of base Likelihood Fit class, it is mostly responsible for likelihood calculation whil...
Definition: LikelihoodFit.h:6
void PrepareFit()
prepare output and perform sanity checks
virtual double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
void RunMCMC() override
Actual implementation of PSO Fit algorithm.
Definition: PSO.cpp:25
PSO(manager *const fitMan)
constructor
Definition: PSO.cpp:6
std::vector< double > ranges_max
Definition: PSO.h:151
double fConstriction
Definition: PSO.h:159
void init()
Definition: PSO.cpp:50
std::vector< double > three_vector_addition(std::vector< double > vec1, const std::vector< double > &vec2, const std::vector< double > &vec3)
Definition: PSO.h:123
double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
Definition: PSO.cpp:530
std::vector< std::vector< double > > bisection(const std::vector< double > &position, const double minimum, const double range, const double precision)
Definition: PSO.cpp:151
double fTwo
Definition: PSO.h:156
std::vector< double > vector_add(const std::vector< double > &v1, const std::vector< double > &v2)
Definition: PSO.h:113
std::vector< double > vector_subtract(const std::vector< double > &v1, const std::vector< double > &v2)
Definition: PSO.h:118
std::vector< double > vector_multiply(std::vector< double > velocity, const double mul)
Definition: PSO.h:105
std::vector< bool > fixed
Definition: PSO.h:150
double fInertia
Definition: PSO.h:154
particle * get_best_particle()
Definition: PSO.h:84
int fParticles
Definition: PSO.h:162
double fOne
Definition: PSO.h:155
int fIterations
Definition: PSO.h:158
std::vector< std::vector< double > > calc_uncertainty(const std::vector< double > &position, const double minimum)
Definition: PSO.cpp:236
void set_best_particle(particle *n)
Definition: PSO.h:87
std::vector< double * > paramlist
Definition: PSO.h:163
std::vector< particle * > system
Definition: PSO.h:153
double vel[kMaxParticles]
Definition: PSO.h:165
double fConvergence
Definition: PSO.h:157
double CalcChi(std::vector< double > x)
Definition: PSO.h:139
double fBestValue
Definition: PSO.h:148
std::vector< std::vector< double > > uncertainties
Definition: PSO.h:160
std::vector< double > ranges_min
Definition: PSO.h:152
double * par
Definition: PSO.h:166
int fDim
Definition: PSO.h:168
void run()
Definition: PSO.cpp:398
void WriteOutput()
Definition: PSO.cpp:443
double swarmIterate()
Definition: PSO.cpp:338
double rastriginFunc(const double *x)
Evaluates the Rastrigin function for a given parameter values.
Definition: PSO.cpp:540
std::vector< double > prior
Definition: PSO.h:149
void uncertainty_check(const std::vector< double > &previous_pos)
Definition: PSO.cpp:300
The manager class is responsible for managing configurations and settings.
Definition: Manager.h:16
YAML::Node const & raw()
Return config.
Definition: Manager.h:41
Class particle - stores the position, velocity and personal best With functions which move particle a...
Definition: PSO.h:17
void set_personal_best_value(const double new_val)
Definition: PSO.h:40
void set_personal_best_position(const std::vector< double > &new_pos)
Definition: PSO.h:32
std::vector< double > get_personal_best_position()
Definition: PSO.h:36
void set_value(const double val)
Definition: PSO.h:59
static constexpr const double _LARGE_LOGL_
Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such s...
Definition: Core.h:70