MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | Static Private Attributes | List of all members
PSO Class Reference

Class PSO, consist of a vector of object Class Particle and global best Takes in the size (number of particle) and number of iteration functions includes: finding global best, updating velocity, actual minimisation function. More...

#include <Fitters/PSO.h>

Inheritance diagram for PSO:
[legend]
Collaboration diagram for PSO:
[legend]

Public Member Functions

 PSO (manager *const fitMan)
 constructor
 
virtual ~PSO ()
 Destructor.
 
particleget_best_particle ()
 
void set_best_particle (particle *n)
 
std::vector< std::vector< double > > bisection (const std::vector< double > &position, const double minimum, const double range, const double precision)
 
std::vector< std::vector< double > > calc_uncertainty (const std::vector< double > &position, const double minimum)
 
void init ()
 
void uncertainty_check (const std::vector< double > &previous_pos)
 
void run ()
 
void WriteOutput ()
 
void RunMCMC () override
 Actual implementation of PSO Fit algorithm.
 
double CalcChi2 (const double *x)
 Chi2 calculation over all included samples and syst objects.
 
double rastriginFunc (const double *x)
 Evaluates the Rastrigin function for a given parameter values.
 
double swarmIterate ()
 
std::vector< double > vector_multiply (std::vector< double > velocity, const double mul)
 
std::vector< double > vector_add (const std::vector< double > &v1, const std::vector< double > &v2)
 
std::vector< double > vector_subtract (const std::vector< double > &v1, const std::vector< double > &v2)
 
std::vector< double > three_vector_addition (std::vector< double > vec1, const std::vector< double > &vec2, const std::vector< double > &vec3)
 
std::vector< double > four_vector_addition (std::vector< double > vec1, const std::vector< double > &vec2, const std::vector< double > &vec3, const std::vector< double > &vec4)
 
double CalcChi (std::vector< double > x)
 
std::string GetName () const
 Get name of class.
 
- Public Member Functions inherited from LikelihoodFit
 LikelihoodFit (manager *const fitMan)
 Constructor.
 
virtual ~LikelihoodFit ()
 Destructor.
 
virtual double CalcChi2 (const double *x)
 Chi2 calculation over all included samples and syst objects.
 
int GetNPars ()
 Get total number of params, this sums over all covariance objects.
 
virtual void RunMCMC ()=0
 Implementation of fitting algorithm.
 
virtual std::string GetName () const
 Get name of class.
 
- Public Member Functions inherited from FitterBase
 FitterBase (manager *const fitMan)
 Constructor.
 
virtual ~FitterBase ()
 Destructor for the FitterBase class.
 
void AddSampleHandler (SampleHandlerBase *sample)
 This function adds a sample PDF object to the analysis framework. The sample PDF object will be utilized in fitting procedures or likelihood scans.
 
void AddSystObj (ParameterHandlerBase *cov)
 This function adds a Covariance object to the analysis framework. The Covariance object will be utilized in fitting procedures or likelihood scans.
 
virtual void RunMCMC ()=0
 The specific fitting algorithm implemented in this function depends on the derived class. It could be Markov Chain Monte Carlo (MCMC), MinuitFit, or another algorithm.
 
void DragRace (const int NLaps=100)
 Calculates the required time for each sample or covariance object in a drag race simulation. Inspired by Dan's feature.
 
void RunLLHScan ()
 Perform a 1D likelihood scan.
 
void GetStepScaleBasedOnLLHScan ()
 LLH scan is good first estimate of step scale.
 
void Run2DLLHScan ()
 Perform a 2D likelihood scan.
 
void RunSigmaVar ()
 Perform a 2D and 1D sigma var for all samples.
 
virtual void StartFromPreviousFit (const std::string &FitName)
 Allow to start from previous fit/chain.
 
virtual std::string GetName () const
 Get name of class.
 

Private Attributes

particlebest_particle
 
double fBestValue
 
std::vector< double > prior
 
std::vector< bool > fixed
 
std::vector< double > ranges_max
 
std::vector< double > ranges_min
 
std::vector< particle * > system
 
double fInertia
 
double fOne
 
double fTwo
 
double fConvergence
 
int fIterations
 
double fConstriction
 
std::vector< std::vector< double > > uncertainties
 
int fParticles
 
std::vector< double * > paramlist
 
double vel [kMaxParticles]
 
double * par
 
int fDim
 

Static Private Attributes

static constexpr const int kMaxParticles = 10000
 

Additional Inherited Members

- Protected Member Functions inherited from LikelihoodFit
void PrepareFit ()
 prepare output and perform sanity checks
 
- Protected Member Functions inherited from FitterBase
void ProcessMCMC ()
 Process MCMC output.
 
void PrepareOutput ()
 Prepare the output file.
 
void SaveOutput ()
 Save output and close files.
 
void SanitiseInputs ()
 Remove obsolete memory and make other checks before fit starts.
 
void SaveSettings ()
 Save the settings that the MCMC was run with.
 
bool GetScaneRange (std::map< std::string, std::vector< double > > &scanRanges)
 YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow.
 
bool CheckSkipParameter (const std::vector< std::string > &SkipVector, const std::string &ParamName) const
 KS: Check whether we want to skip parameter using skip vector.
 
- Protected Attributes inherited from LikelihoodFit
int NPars
 Number of all parameters from all covariances.
 
int NParsPCA
 Number of all parameters from all covariances in PCA base.
 
bool fMirroring
 Flag telling if mirroring is used or not.
 
- Protected Attributes inherited from FitterBase
managerfitMan
 The manager.
 
unsigned int step
 current state
 
double logLCurr
 current likelihood
 
double logLProp
 proposed likelihood
 
double accProb
 current acceptance prob
 
int accCount
 counts accepted steps
 
int stepStart
 step start if restarting
 
std::vector< double > sample_llh
 store the llh breakdowns
 
std::vector< double > syst_llh
 systematic llh breakdowns
 
std::vector< SampleHandlerBase * > samples
 Sample holder.
 
unsigned int TotalNSamples
 Total number of samples used.
 
std::vector< ParameterHandlerBase * > systematics
 Systematic holder.
 
std::unique_ptr< TStopwatch > clock
 tells global time how long fit took
 
std::unique_ptr< TStopwatch > stepClock
 tells how long single step/fit iteration took
 
double stepTime
 Time of single step.
 
std::unique_ptr< TRandom3 > random
 Random number.
 
TFile * outputFile
 Output.
 
TDirectory * CovFolder
 Output cov folder.
 
TDirectory * SampleFolder
 Output sample folder.
 
TTree * outTree
 Output tree with posteriors.
 
int auto_save
 auto save every N steps
 
bool fTestLikelihood
 Necessary for some fitting algorithms like PSO.
 
bool FileSaved
 Checks if file saved not repeat some operations.
 
bool SettingsSaved
 Checks if setting saved not repeat some operations.
 
bool OutputPrepared
 Checks if output prepared not repeat some operations.
 

Detailed Description

Class PSO, consist of a vector of object Class Particle and global best Takes in the size (number of particle) and number of iteration functions includes: finding global best, updating velocity, actual minimisation function.

Author
Emily Ip
Mark Scott
Date
24/2/2023

Definition at line 77 of file PSO.h.

Constructor & Destructor Documentation

◆ PSO()

PSO::PSO ( manager *const  fitMan)

constructor

Definition at line 6 of file PSO.cpp.

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}
MaCh3Plotting::PlottingManager * man
manager * fitMan
The manager.
Definition: FitterBase.h:92
bool fTestLikelihood
Necessary for some fitting algorithms like PSO.
Definition: FitterBase.h:142
Implementation of base Likelihood Fit class, it is mostly responsible for likelihood calculation whil...
Definition: LikelihoodFit.h:6
double fConstriction
Definition: PSO.h:159
double fTwo
Definition: PSO.h:156
double fInertia
Definition: PSO.h:154
int fParticles
Definition: PSO.h:162
double fOne
Definition: PSO.h:155
int fIterations
Definition: PSO.h:158
double fConvergence
Definition: PSO.h:157
int fDim
Definition: PSO.h:168
YAML::Node const & raw()
Return config.
Definition: Manager.h:41

◆ ~PSO()

virtual PSO::~PSO ( )
inlinevirtual

Destructor.

Definition at line 82 of file PSO.h.

82{};

Member Function Documentation

◆ bisection()

std::vector< std::vector< double > > PSO::bisection ( const std::vector< double > &  position,
const double  minimum,
const double  range,
const double  precision 
)

Definition at line 151 of file PSO.cpp.

152 {
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}
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_TRACE
Definition: MaCh3Logger.h:21
double CalcChi(std::vector< double > x)
Definition: PSO.h:139

◆ calc_uncertainty()

std::vector< std::vector< double > > PSO::calc_uncertainty ( const std::vector< double > &  position,
const double  minimum 
)

Definition at line 236 of file PSO.cpp.

236 {
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}

◆ CalcChi()

double PSO::CalcChi ( std::vector< double >  x)
inline

Definition at line 139 of file PSO.h.

139 {
140 double* a = &x[0];
141 return CalcChi2(a);
142 };
double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
Definition: PSO.cpp:530

◆ CalcChi2()

double PSO::CalcChi2 ( const double *  x)
virtual

Chi2 calculation over all included samples and syst objects.

Reimplemented from LikelihoodFit.

Definition at line 530 of file PSO.cpp.

530 {
531// *******************
532 if(fTestLikelihood) {
533 return rastriginFunc(x);
534 } else {
535 return LikelihoodFit::CalcChi2(x);
536 }
537}
virtual double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
double rastriginFunc(const double *x)
Evaluates the Rastrigin function for a given parameter values.
Definition: PSO.cpp:540

◆ four_vector_addition()

std::vector< double > PSO::four_vector_addition ( std::vector< double >  vec1,
const std::vector< double > &  vec2,
const std::vector< double > &  vec3,
const std::vector< double > &  vec4 
)
inline

Definition at line 131 of file PSO.h.

132 {
133 for (size_t i = 0; i < vec1.size(); ++i) {
134 vec1[i] += vec2[i] + vec3[i] + vec4[i];
135 }
136 return vec1;
137 };

◆ get_best_particle()

particle * PSO::get_best_particle ( )
inline

Definition at line 84 of file PSO.h.

84 {
85 return best_particle;
86 }
particle * best_particle
Definition: PSO.h:147

◆ GetName()

std::string PSO::GetName ( ) const
inlinevirtual

Get name of class.

Reimplemented from LikelihoodFit.

Definition at line 144 of file PSO.h.

144{return "PSO";};

◆ init()

void PSO::init ( )

Definition at line 50 of file PSO.cpp.

50 {
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}
std::unique_ptr< TRandom3 > random
Random number.
Definition: FitterBase.h:128
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Definition: FitterBase.h:118
std::vector< double > ranges_max
Definition: PSO.h:151
std::vector< bool > fixed
Definition: PSO.h:150
void set_best_particle(particle *n)
Definition: PSO.h:87
std::vector< particle * > system
Definition: PSO.h:153
double fBestValue
Definition: PSO.h:148
std::vector< double > ranges_min
Definition: PSO.h:152
std::vector< double > prior
Definition: PSO.h:149
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
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

◆ rastriginFunc()

double PSO::rastriginFunc ( const double *  x)

Evaluates the Rastrigin function for a given parameter values.

Definition at line 540 of file PSO.cpp.

540 {
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}
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

◆ run()

void PSO::run ( )

Definition at line 398 of file PSO.cpp.

398 {
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}
int accCount
counts accepted steps
Definition: FitterBase.h:103
unsigned int step
current state
Definition: FitterBase.h:95
double logLCurr
current likelihood
Definition: FitterBase.h:97
int auto_save
auto save every N steps
Definition: FitterBase.h:139
TTree * outTree
Output tree with posteriors.
Definition: FitterBase.h:137
std::vector< std::vector< double > > bisection(const std::vector< double > &position, const double minimum, const double range, const double precision)
Definition: PSO.cpp:151
particle * get_best_particle()
Definition: PSO.h:84
std::vector< std::vector< double > > uncertainties
Definition: PSO.h:160
double swarmIterate()
Definition: PSO.cpp:338

◆ RunMCMC()

void PSO::RunMCMC ( )
overridevirtual

Actual implementation of PSO Fit algorithm.

Implements LikelihoodFit.

Definition at line 25 of file PSO.cpp.

25 {
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}
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
Definition: FitterBase.cpp:213
void PrepareFit()
prepare output and perform sanity checks
void init()
Definition: PSO.cpp:50
std::vector< double * > paramlist
Definition: PSO.h:163
double vel[kMaxParticles]
Definition: PSO.h:165
double * par
Definition: PSO.h:166
void run()
Definition: PSO.cpp:398
void WriteOutput()
Definition: PSO.cpp:443

◆ set_best_particle()

void PSO::set_best_particle ( particle n)
inline

Definition at line 87 of file PSO.h.

87 {
88 best_particle = n;
89 }

◆ swarmIterate()

double PSO::swarmIterate ( )

Definition at line 338 of file PSO.cpp.

338 {
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}
std::vector< double > three_vector_addition(std::vector< double > vec1, const std::vector< double > &vec2, const std::vector< double > &vec3)
Definition: PSO.h:123
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< double > get_personal_best_position()
Definition: PSO.h:36

◆ three_vector_addition()

std::vector< double > PSO::three_vector_addition ( std::vector< double >  vec1,
const std::vector< double > &  vec2,
const std::vector< double > &  vec3 
)
inline

Definition at line 123 of file PSO.h.

125 {
126 for (size_t i = 0; i < vec1.size(); ++i) {
127 vec1[i] += vec2[i] + vec3[i];
128 }
129 return vec1;
130 };

◆ uncertainty_check()

void PSO::uncertainty_check ( const std::vector< double > &  previous_pos)

Definition at line 300 of file PSO.cpp.

300 {
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}

◆ vector_add()

std::vector< double > PSO::vector_add ( const std::vector< double > &  v1,
const std::vector< double > &  v2 
)
inline

Definition at line 113 of file PSO.h.

113 {
114 std::vector<double> v3;
115 transform(v1.begin(), v1.end(), v2.begin(), back_inserter(v3), std::plus<double>());
116 return v3;
117 };

◆ vector_multiply()

std::vector< double > PSO::vector_multiply ( std::vector< double >  velocity,
const double  mul 
)
inline

Definition at line 105 of file PSO.h.

105 {
106 // std::bind1st deprecated since C++11, removed in c++17
107 // transform(velocity.begin(),velocity.end(),velocity.begin(),std::bind1st(std::multiplies<double>(),mul));
108 std::transform(velocity.begin(), velocity.end(), velocity.begin(),
109 std::bind(std::multiplies<double>(), mul, std::placeholders::_1));
110 return velocity;
111 };

◆ vector_subtract()

std::vector< double > PSO::vector_subtract ( const std::vector< double > &  v1,
const std::vector< double > &  v2 
)
inline

Definition at line 118 of file PSO.h.

118 {
119 std::vector<double> v3 ;
120 transform(v1.begin(), v1.end(), v2.begin(), back_inserter(v3), std::minus<double>());
121 return v3;
122 };

◆ WriteOutput()

void PSO::WriteOutput ( )

Definition at line 443 of file PSO.cpp.

443 {
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}
void SaveOutput()
Save output and close files.
Definition: FitterBase.cpp:221
TFile * outputFile
Output.
Definition: FitterBase.h:131

Member Data Documentation

◆ best_particle

particle* PSO::best_particle
private

Definition at line 147 of file PSO.h.

◆ fBestValue

double PSO::fBestValue
private

Definition at line 148 of file PSO.h.

◆ fConstriction

double PSO::fConstriction
private

Definition at line 159 of file PSO.h.

◆ fConvergence

double PSO::fConvergence
private

Definition at line 157 of file PSO.h.

◆ fDim

int PSO::fDim
private

Definition at line 168 of file PSO.h.

◆ fInertia

double PSO::fInertia
private

Definition at line 154 of file PSO.h.

◆ fIterations

int PSO::fIterations
private

Definition at line 158 of file PSO.h.

◆ fixed

std::vector<bool> PSO::fixed
private

Definition at line 150 of file PSO.h.

◆ fOne

double PSO::fOne
private

Definition at line 155 of file PSO.h.

◆ fParticles

int PSO::fParticles
private

Definition at line 162 of file PSO.h.

◆ fTwo

double PSO::fTwo
private

Definition at line 156 of file PSO.h.

◆ kMaxParticles

constexpr const int PSO::kMaxParticles = 10000
staticconstexprprivate

Definition at line 164 of file PSO.h.

◆ par

double* PSO::par
private

Definition at line 166 of file PSO.h.

◆ paramlist

std::vector<double*> PSO::paramlist
private

Definition at line 163 of file PSO.h.

◆ prior

std::vector<double> PSO::prior
private

Definition at line 149 of file PSO.h.

◆ ranges_max

std::vector<double> PSO::ranges_max
private

Definition at line 151 of file PSO.h.

◆ ranges_min

std::vector<double> PSO::ranges_min
private

Definition at line 152 of file PSO.h.

◆ system

std::vector<particle*> PSO::system
private

Definition at line 153 of file PSO.h.

◆ uncertainties

std::vector<std::vector<double> > PSO::uncertainties
private

Definition at line 160 of file PSO.h.

◆ vel

double PSO::vel[kMaxParticles]
private

Definition at line 165 of file PSO.h.


The documentation for this class was generated from the following files: