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 More...

virtual	~PSO ()
	Destructor. More...

particle *	get_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 () final
	Actual implementation of PSO Fit algorithm. More...

double	CalcChi2 (const double *x) final
	Chi2 calculation over all included samples and syst objects. More...

double	rastriginFunc (const double *x)
	Evaluates the Rastrigin function for a given parameter values. More...

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)

Public Member Functions inherited from LikelihoodFit
	LikelihoodFit (Manager *const fitMan)
	Constructor. More...

virtual	~LikelihoodFit ()
	Destructor. More...

int	GetNPars () const
	Get total number of params, this sums over all covariance objects. More...

Public Member Functions inherited from FitterBase
	FitterBase (Manager *const fitMan)
	Constructor. More...

virtual	~FitterBase ()
	Destructor for the FitterBase class. More...

void	AddSampleHandler (SampleHandlerInterface *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. More...

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. More...

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. More...

void	RunLLHScan ()
	Perform a 1D likelihood scan. More...

void	RunLLHMap ()
	Perform a general multi-dimensional likelihood scan. More...

void	GetStepScaleBasedOnLLHScan (const std::string &filename="")
	LLH scan is good first estimate of step scale. More...

void	Run2DLLHScan ()
	Perform a 2D likelihood scan. More...

void	RunSigmaVar ()
	Perform a 1D/2D sigma var for all samples. More...

virtual void	StartFromPreviousFit (const std::string &FitName)
	Allow to start from previous fit/chain. More...

std::string	GetName () const
	Get name of class. More...

Private Attributes
particle *	best_particle

double	fBestValue

std::vector< double >	prior

std::vector< bool >	fixed

std::vector< double >	ranges_max

std::vector< double >	ranges_min

std::vector< std::unique_ptr< particle > >	system

double	fInertia

double	fOne

double	fTwo

double	fConvergence

int	fIterations

double	fConstriction

std::vector< std::vector< double > >	uncertainties

int	fParticles

std::vector< std::vector< double > >	paramlist

double	vel [kMaxParticles]

int	fDim

Static Private Attributes
constexpr static const int	kMaxParticles = 10000

Additional Inherited Members
Protected Member Functions inherited from LikelihoodFit
void	PrepareFit ()
	prepare output and perform sanity checks More...

Protected Member Functions inherited from FitterBase
void	ProcessMCMC ()
	Process MCMC output. More...

void	PrepareOutput ()
	Prepare the output file. More...

void	SaveOutput ()
	Save output and close files. More...

void	SanitiseInputs ()
	Remove obsolete memory and make other checks before fit starts. More...

void	SaveSettings ()
	Save the settings that the MCMC was run with. More...

bool	GetScanRange (std::map< std::string, std::vector< double >> &scanRanges) const
	YSP: Set up a mapping to store parameters with user-specified ranges, suggested by D. Barrow. More...

void	GetParameterScanRange (const ParameterHandlerBase *cov, const int i, double &CentralValue, double &lower, double &upper, const int n_points, const std::string &suffix="") const
	Helper function to get parameter scan range, central value. More...

bool	CheckSkipParameter (const std::vector< std::string > &SkipVector, const std::string &ParamName) const
	KS: Check whether we want to skip parameter using skip vector. More...

void	CustomRange (const std::string &ParName, const double sigma, double &ParamShiftValue) const
	For comparison with other fitting frameworks (like P-Theta) we usually have to apply different parameter values then usual 1, 3 sigma. More...

Protected Attributes inherited from LikelihoodFit
int	NPars
	Number of all parameters from all covariances. More...

int	NParsPCA
	Number of all parameters from all covariances in PCA base. More...

bool	fMirroring
	Flag telling if mirroring is used or not. More...

Protected Attributes inherited from FitterBase
Manager *	fitMan
	The manager for configuration handling. More...

unsigned int	step
	current state More...

double	logLCurr
	current likelihood More...

double	logLProp
	proposed likelihood More...

double	accProb
	current acceptance prob More...

int	accCount
	counts accepted steps More...

unsigned int	stepStart
	step start, by default 0 if we start from previous chain then it will be different More...

std::vector< double >	sample_llh
	store the llh breakdowns More...

std::vector< double >	syst_llh
	systematic llh breakdowns More...

std::vector< SampleHandlerInterface * >	samples
	Sample holder. More...

unsigned int	TotalNSamples
	Total number of samples used, single SampleHandler can store more than one analysis sample! More...

std::vector< ParameterHandlerBase * >	systematics
	Systematic holder. More...

std::unique_ptr< TStopwatch >	clock
	tells global time how long fit took More...

std::unique_ptr< TStopwatch >	stepClock
	tells how long single step/fit iteration took More...

double	stepTime
	Time of single step. More...

std::unique_ptr< TRandom3 >	random
	Random number. More...

TFile *	outputFile
	Output. More...

TDirectory *	CovFolder
	Output cov folder. More...

TDirectory *	SampleFolder
	Output sample folder. More...

TTree *	outTree
	Output tree with posteriors. More...

int	auto_save
	auto save every N steps More...

bool	fTestLikelihood
	Necessary for some fitting algorithms like PSO. More...

bool	FileSaved
	Checks if file saved not repeat some operations. More...

bool	SettingsSaved
	Checks if setting saved not repeat some operations. More...

bool	OutputPrepared
	Checks if output prepared not repeat some operations. More...

std::string	AlgorithmName
	Name of fitting algorithm that is being used. More...

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 79 of file PSO.h.

Constructor & Destructor Documentation

◆ PSO()

PSO::PSO ( Manager *const fitMan )

constructor

Definition at line 6 of file PSO.cpp.

                      : LikelihoodFit(man) {
 // ***************
   AlgorithmName = "PSO";
   fConstriction = Get<double>(fitMan->raw()["General"]["PSO"]["Constriction"], __FILE__, __LINE__);
   fInertia = Get<double>(fitMan->raw()["General"]["PSO"]["Inertia"], __FILE__, __LINE__) * fConstriction;
   fOne = Get<double>(fitMan->raw()["General"]["PSO"]["One"], __FILE__, __LINE__) * fConstriction;
   fTwo = Get<double>(fitMan->raw()["General"]["PSO"]["Two"], __FILE__, __LINE__) * fConstriction;
   fParticles = Get<int>(fitMan->raw()["General"]["PSO"]["Particles"], __FILE__, __LINE__);
   fIterations = Get<int>(fitMan->raw()["General"]["PSO"]["Iterations"], __FILE__, __LINE__);
   fConvergence = Get<double>(fitMan->raw()["General"]["PSO"]["Convergence"], __FILE__, __LINE__);
  
   fDim = 0;
  
   if(fTestLikelihood)
   {
     fDim = Get<int>(fitMan->raw()["General"]["PSO"]["TestLikelihoodDim"], __FILE__, __LINE__);
   }
 }

◆ ~PSO()

virtual PSO::~PSO ( )

inlinevirtual

Destructor.

Definition at line 84 of file PSO.h.

84 {};

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 158 of file PSO.cpp.

                                                                                              {
 // *************************
   std::vector<std::vector<double>> uncertainties_list;
   for (unsigned int i = 0; i< position.size(); ++i) {
     MACH3LOG_INFO("{}", i);
     //std::vector<double> uncertainties;
     std::vector<double> new_position = position; new_position[i] = position[i]-range;
     double val_1 = CalcChi(new_position)-minimum-1.0;
     while (val_1*-1.0 > 0.0){
       new_position[i] -= range;
       val_1 = CalcChi(new_position)-minimum-1.0;
     }
     std::vector<double> bisect_position = position; bisect_position[i] = bisect_position[i] - (position[i]-new_position[i])/2;
     std::vector<std::vector<double>> position_list{new_position,bisect_position,position};
     double val_2 = CalcChi(bisect_position)-minimum-1.0;
     std::vector<double> value_list{val_1,val_2, -1.0};
     double res = 1.0;
     while (res > precision){
       if (value_list[0] * value_list[1] < 0){
         std::vector<double> new_bisect_position = position_list[0];
         new_bisect_position[i] =new_bisect_position[i]+ (position_list[1][i]-position_list[0][i])/2;
         double new_val = CalcChi(new_bisect_position)-minimum-1.0;
         position_list[2] = position_list[1];
         value_list[2] = value_list[1];
         value_list[1] = new_val;
         position_list[1] = new_bisect_position;
         res = std::abs(position[2]-position[0]);
       }
       else{
         std::vector<double> new_bisect_position = position_list[1];
         new_bisect_position[i] += (position_list[2][i]-position_list[1][i])/2;
         double new_val = CalcChi(new_bisect_position)-minimum-1.0;
         position_list[0] = position_list[1];
         value_list[0] = value_list[1];
         value_list[1] = new_val;
         position_list[1] = new_bisect_position;
         res = std::abs(position_list[2][i]-position_list[1][i]);
       }
     }
     //do the same thing for position uncertainty
     std::vector<double> new_position_p = position; new_position_p[i] = position[i]+range;
     double val_1_p = CalcChi(new_position_p)-minimum-1.0;
     while (val_1_p * -1.0 > 0.0){
       new_position_p[i] += range;
       val_1_p = CalcChi(new_position_p)-minimum-1.0;
     }
     std::vector<double> bisect_position_p = position; bisect_position_p[i] = bisect_position_p[i] += (new_position_p[i]-position[i])/2;
     std::vector<std::vector<double>> position_list_p{position,bisect_position_p,new_position_p};
     double val_2_p = CalcChi(bisect_position_p)-minimum-1.0;
     std::vector<double> value_list_p{-1.0,val_2_p, val_1_p};
     double res_p = 1.0;
     while (res_p > precision){
       if (value_list_p[0] * value_list_p[1] < 0){
         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;
         double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
         position_list_p[2] = position_list_p[1];
         value_list_p[2] = value_list_p[1];
         value_list_p[1] = new_val_p;
         position_list_p[1] = new_bisect_position_p;
         res = std::abs(position[2]-position[0]);
         res_p = std::abs(position_list_p[1][i]-position_list_p[0][i]);
         MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]);
       }
       else{
         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;
         double new_val_p = CalcChi(new_bisect_position_p)-minimum-1.0;
         position_list_p[0] = position_list_p[1];
         value_list_p[0] = value_list_p[1];
         value_list_p[1] = new_val_p;
         position_list_p[1] = new_bisect_position_p;
         res_p = std::abs(position_list_p[2][i]-position_list_p[1][i]);
         MACH3LOG_TRACE("Pos midpoint is {:.2f}", position_list_p[1][i]);
       }
     }
     uncertainties_list.push_back({std::abs(position[i]-position_list[1][i]),std::abs(position[i]-position_list_p[1][i])});
     MACH3LOG_INFO("Uncertainty finished for d = {}", i);
     MACH3LOG_INFO("LLR values for ± positive and negative uncertainties are {:<10.2f} and {:<10.2f}",
                   CalcChi(position_list[1]), CalcChi(position_list_p[1]));
   }
   return uncertainties_list;
 }

◆ calc_uncertainty()

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

Definition at line 242 of file PSO.cpp.

                                                                                                           {
 // *************************
   std::vector<double> pos_uncertainty(position.size());
   std::vector<double> neg_uncertainty(position.size());
   constexpr int num = 200;
   std::vector<double> pos = position;
   for (unsigned int i = 0; i < position.size(); ++i) {
     double curr_ival = pos[i];
  
     double neg_stop = position[i] - 5e-2;
     double pos_stop = position[i] + 5e-2;
     double start = position[i];
     std::vector<double> x(num);
     std::vector<double> y(num);
     double StepPoint = (start-neg_stop) / (num - 1);
     double value = start;
     for (int j = 0; j < num; ++j) {
       pos[i] = value;
       double LLR = CalcChi(position) - minimum - 1.0;
       x[j] = value;
       y[j] = LLR;
       value -= StepPoint;
     }
     pos[i] = curr_ival;
  
     int closest_index = 0;
     double closest_value = std::abs(y[0]); // Initialize with the first element
     for (unsigned int ii = 1; ii < y.size(); ++ii) {
       double abs_y = std::abs(y[ii]);
       if (abs_y < closest_value) {
         closest_index = ii;
         closest_value = abs_y;
       }
     }
     neg_uncertainty[i] = x[closest_index];
     MACH3LOG_INFO("Neg");
     x.assign(num, 0);
     y.assign(num, 0);
     StepPoint = (pos_stop-start) / (num - 1);
     value = start;
     for (int j = 0; j < num; ++j) {
       pos[i] = value;
       double LLR = CalcChi(position) - minimum - 1.0;
       x[j] = value;
       y[j] = LLR;
       value += StepPoint;
     }
     pos[i] = curr_ival;
     closest_index = 0;
     closest_value = std::abs(y[0]); // Initialize with the first element
     for (unsigned int ii = 1; ii < y.size(); ++ii) {
       double abs_y = std::abs(y[ii]);
       if (abs_y < closest_value) {
         closest_index = ii;
         closest_value = abs_y;
       }
     }
     pos_uncertainty[i] = x[closest_index];
   }
   std::vector<std::vector<double>> res{neg_uncertainty,pos_uncertainty};
   return res;
 }

◆ CalcChi()

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

inline

Definition at line 141 of file PSO.h.

                                        {
       double* a = &x[0];
       return CalcChi2(a);
     };

◆ CalcChi2()

double PSO::CalcChi2 ( const double * x )

finalvirtual

Chi2 calculation over all included samples and syst objects.

Reimplemented from LikelihoodFit.

Definition at line 536 of file PSO.cpp.

                                     {
 // *******************
   if(fTestLikelihood) {
     return rastriginFunc(x);
   } else {
     return LikelihoodFit::CalcChi2(x);
   }
 }

◆ 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 133 of file PSO.h.

                                                                                                           {
       for (size_t i = 0; i < vec1.size(); ++i) {
         vec1[i] += vec2[i] + vec3[i] + vec4[i];
       }
       return vec1;
     };

◆ get_best_particle()

particle* PSO::get_best_particle ( )

inline

Definition at line 86 of file PSO.h.

                                  {
         return best_particle;
     }

◆ init()

void PSO::init ( )

Definition at line 58 of file PSO.cpp.

               {
 // *************************
   fBestValue = M3::_LARGE_LOGL_;
  
   //KS: For none PCA this will be equal to normal parameters
   //const int NparsPSOFull = NPars;
   //const int NparsPSO = NParsPCA;
  
   MACH3LOG_INFO("Preparing PSO");
   // Initialise bounds on parameters
   if(fTestLikelihood){
     for (int i = 0; i < fDim; i++){
       // Test function ranges
       ranges_min.push_back(-5);
       ranges_max.push_back(5);
       fixed.push_back(0);
     }
   }
   else{
     for (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it){
       if(!(*it)->IsPCA())
       {
         fDim += (*it)->GetNumParams();
         for(int i = 0; i < (*it)->GetNumParams(); ++i)
         {
           double curr = (*it)->GetParPreFit(i);
           double lim = 10.0*(*it)->GetDiagonalError(i);
           double low = (*it)->GetLowerBound(i);
           double high = (*it)->GetUpperBound(i);
           if(low > curr - lim) ranges_min.push_back(low);
           else ranges_min.push_back(curr - lim);
           if(high < curr + lim) ranges_min.push_back(high);
           else ranges_min.push_back(curr + lim);
           prior.push_back(curr);
  
           if((*it)->IsParameterFixed(i)){
             fixed.push_back(1);
           }
           else{
             fixed.push_back(0);
           }
         }
       }
       else
       {
         fDim += (*it)->GetNParameters();
         for(int i = 0; i < (*it)->GetNParameters(); ++i)
         {
           ranges_min.push_back(-100.0);
           ranges_max.push_back(100.0);
           prior.push_back((*it)->GetParPreFit(i));
           if((*it)->GetPCAHandler()->IsParameterFixedPCA(i)){
             fixed.push_back(1);
           }
           else{
             fixed.push_back(0);
           }
         }
       }
     }
   }
  
   MACH3LOG_INFO("Printing Minimums and Maximums of Variables to be minimized");
   for (int i = 0; i < fDim; i++){
     MACH3LOG_INFO("Variable {} : {:.2f}, {:.2f}", i, ranges_min[i], ranges_max[i]);
   }
  
   // Initialise particle positions
   for (int i = 0; i < fParticles; ++i){
     std::vector<double> init_position;
     std::vector<double> init_velocity;
  
     //Initialising in +/- 5sigma of prior value from BANFF interface
     for (int j=0; j<fDim; ++j){
       if(fixed[j]){
         init_position.push_back(prior[j]);
         init_velocity.push_back(0.0);
       }
       else{
         double dist = fabs(ranges_max[j]-ranges_min[j]);
         //Initialise to random position uniform in space
         init_position.push_back(ranges_min[j] + random->Rndm()*dist);
         //Initialise velocity to random position uniform in space
         init_velocity.push_back((2.0*random->Rndm()-1.0));//*dist);
       }
     }
     system.emplace_back(std::make_unique<particle>(init_position, init_velocity));
     auto& new_particle = system.back();
     new_particle->set_personal_best_position(init_position);
     double new_value = CalcChi(init_position);
     new_particle->set_personal_best_value(new_value);
     new_particle->set_value(new_value);
     if(new_value < fBestValue){
       fBestValue = new_value;
       set_best_particle(system.back().get());
     }
   }
 }

◆ rastriginFunc()

double PSO::rastriginFunc ( const double * x )

Evaluates the Rastrigin function for a given parameter values.

Definition at line 546 of file PSO.cpp.

                                          {
 // *************************
   stepClock->Start();
  
   //Search range: [-5.12, 5.12]
   constexpr double A = 10.0;
   double sum = 0.0;
   for (int i = 0; i < fDim; ++i) {
       sum += x[i] * x[i] - A * cos(2.0 * 3.14 * x[i]);
   }
   double llh = A * fDim + sum;
  
   accProb = 1;
  
   stepClock->Stop();
   stepTime = stepClock->RealTime();
  
   return llh;
 }

◆ run()

void PSO::run ( )

Definition at line 404 of file PSO.cpp.

               {
 // *************************
   double mean_dist_sq = 0;
  
   int iter = 0;
   for(int i = 0; i < fIterations; ++i, ++iter){
     mean_dist_sq = swarmIterate();
     //double meanVel = std::accumulate(vel, vel + fParticles, 0) / fParticles;
  
     // Weight inertia randomly but scaled by total distance of swarm from global minimum - proxy for total velocity
     // fWeight = ((random->Rndm()+1.0)*0.5)*(10.0/meanVel);
  
     logLCurr = fBestValue;
  
     outTree->Fill();
     // Auto save the output
     if (step % auto_save == 0) outTree->AutoSave();
     step++;
     accCount++;
  
     if (i%100 == 0){
       MACH3LOG_INFO("Mean Dist Sq = {:.2f}", mean_dist_sq);
       MACH3LOG_INFO("Current LLR = {:.2f}", fBestValue);
       MACH3LOG_INFO("Position:");
       for (int j = 0; j < fDim; ++j){
         MACH3LOG_INFO("    Dim {} = {:<10.2f}", j, get_best_particle()->get_personal_best_position()[j]);
       }
     }
     if(fConvergence > 0.0){
       if(mean_dist_sq < fConvergence){
         break;
       }
     }
   }
   MACH3LOG_INFO("Finished after {} runs out of {}", iter, fIterations);
   MACH3LOG_INFO("Mean Dist: {:.2f}", mean_dist_sq);
   MACH3LOG_INFO("Best LLR: {:.2f}", get_best_particle()->get_personal_best_value());
   uncertainties = bisection(get_best_particle()->get_personal_best_position(),get_best_particle()->get_personal_best_value(),0.5,0.005);
   MACH3LOG_INFO("Position for Global Minimum = ");
   for (int i = 0; i< fDim; ++i){
     MACH3LOG_INFO("    Dim {} = {:<10.2f} +{:.2f}, -{:.2f}", i, get_best_particle()->get_personal_best_position()[i], uncertainties[i][1], uncertainties[i][0]);
   }
 }

◆ RunMCMC()

void PSO::RunMCMC ( )

finalvirtual

Actual implementation of PSO Fit algorithm.

Implements FitterBase.

Definition at line 26 of file PSO.cpp.

                  {
 // ***************
   PrepareFit();
  
   // Remove obsolete memory and make other checks before fit starts
   SanitiseInputs();
  
   // Sanitise the adaptive MCMC
   for (const auto &syst : systematics) {
     if (syst->GetDoAdaption()) {
       MACH3LOG_ERROR("Param Handler {} has enabled Adaption, this is not needed for {} so please turn it off", syst->GetName(), GetName());
       throw MaCh3Exception(__FILE__ , __LINE__ );
     }
   }
  
   if(fTestLikelihood){
     outTree->Branch("nParts", &fParticles, "nParts/I");
     for(int i = 0; i < fDim; ++i){
       paramlist.emplace_back(fParticles);
       outTree->Branch(Form("Parameter_%d", i), paramlist[i].data(), Form("Parameter_%d[nParts]/D",i));
     }
 //  vel = new double[fParticles];
     outTree->Branch("vel", vel, "vel[nParts]/D");
   }
  
   init();
   run();
   WriteOutput();
   return;
 }

◆ set_best_particle()

void PSO::set_best_particle ( particle * n )

inline

Definition at line 89 of file PSO.h.

                                        {
       best_particle = n;
     }

◆ swarmIterate()

double PSO::swarmIterate ( )

Definition at line 344 of file PSO.cpp.

                         {
 // *************************
   std::vector<double> total_pos(fDim,0.0);
  
   for (int i = 0; i < fParticles; ++i) {
     std::vector<double> part1 = vector_multiply(system[i]->get_velocity(), fInertia);
     std::vector<double> part2 = vector_multiply(vector_subtract(system[i]->get_personal_best_position(), system[i]->get_position()), (fOne * random->Rndm()));
     std::vector<double> part3 = vector_multiply(vector_subtract(get_best_particle()->get_personal_best_position(), system[i]->get_position()),(fTwo * random->Rndm()));
     std::vector<double> new_velocity = three_vector_addition(part1, part2, part3);
     std::vector<double> new_pos = vector_add(system[i]->get_position(), new_velocity);
     transform(total_pos.begin(), total_pos.end(), new_pos.begin(), total_pos.begin(),[](double x, double y) {return x+y;});
  
     for (int j = 0; j < fDim; ++j) {
       // Check if out of bounds and reflect if so
       if(ranges_min[j] > new_pos[j]){
           new_pos[j] = ranges_min[j];
       }
       else if(new_pos[j] > ranges_max[j]) {
           new_pos[j] = ranges_max[j];
       }
       // If parameter fixed don't update it
       if(fixed[j]) new_pos[j] = system[i]->get_position()[j];
     }
  
     if(fTestLikelihood){
       double velo = 0.0;
       for (int j = 0; j < fDim; ++j) {
         paramlist[j][i] = new_pos[j];
         velo += new_velocity[j]*new_velocity[j];
       }
       vel[i] = sqrt(velo);
     }
  
     system[i]->set_velocity(new_velocity);
     system[i]->set_position(new_pos);
     double new_value = CalcChi(new_pos);
     if(new_value <= system[i]->get_personal_best_value()) {
       system[i]->set_personal_best_value(new_value);
       system[i]->set_personal_best_position(new_pos);
       if(new_value < fBestValue){
         fBestValue = new_value;
         set_best_particle(system[i].get());
       }
     }
   }
  
   std::vector<double> best_pos = get_best_particle()->get_personal_best_position();
   std::vector<double> result(best_pos.size(), 0.0);
   transform(total_pos.begin(), total_pos.end(), total_pos.begin(), [this](double val){return val/fParticles;});
   transform(total_pos.begin(),total_pos.end(),best_pos.begin(),result.begin(),[](double x, double y) {return x-y;});
  
   double mean_dist_sq = 0;
   for (int i = 0; i<fDim; i++){
     mean_dist_sq += result[i]*result[i];
   }
  
   return mean_dist_sq;
 }

◆ 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 125 of file PSO.h.

                                                                              {
       for (size_t i = 0; i < vec1.size(); ++i) {
         vec1[i] += vec2[i] + vec3[i];
       }
       return vec1;
     };

◆ uncertainty_check()

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

Definition at line 306 of file PSO.cpp.

                                                                 {
 // *************************
   std::vector<std::vector<double >> x_list;
   std::vector<std::vector<double >> y_list;
   std::vector<double> position = previous_pos;
   constexpr int num = 5000;
   for (unsigned int i = 0;i<previous_pos.size();++i){
     double curr_ival = position[i];
     double start = previous_pos[i] - 1e-1;
     double stop = previous_pos[i] + 1e-1;
     std::vector<double> x(num);
     std::vector<double> y(num);
     double StepPoint = (stop - start) / (num - 1);
     double value = start;
     MACH3LOG_TRACE("result for fDim: {}", 1);
     for (int j = 0;j < num; ++j) {
       position[i] = value;
       double LLR = CalcChi(position);
       x[j] = value;
       y[j] = LLR;
       value += StepPoint;
     }
     position[i] = curr_ival;
     MACH3LOG_INFO("");
     MACH3LOG_INFO("For fDim{} x list is", i+1);
     for (unsigned int k = 0; k < x.size(); ++k){
       MACH3LOG_INFO("  {}", x[k]);
     }
     MACH3LOG_INFO("");
     MACH3LOG_INFO("For fDim{} y list is", i+1);
     for (unsigned int k = 0; k < x.size(); ++k){
       MACH3LOG_INFO("  {}", y[k]);
     }
     MACH3LOG_INFO("");
   }
 }

◆ vector_add()

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

inline

Definition at line 115 of file PSO.h.

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

◆ vector_multiply()

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

inline

Definition at line 107 of file PSO.h.

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

◆ vector_subtract()

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

inline

Definition at line 120 of file PSO.h.

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

◆ WriteOutput()

void PSO::WriteOutput ( )

Definition at line 449 of file PSO.cpp.

                      {
 // *************************
   outputFile->cd();
  
   TVectorD* PSOParValue = new TVectorD(fDim);
   TVectorD* PSOParError = new TVectorD(fDim);
  
   for(int i = 0; i < fDim; ++i)
   {
     (*PSOParValue)(i) = 0;
     (*PSOParError)(i) = 0;
   }
  
   std::vector<double> minimum = get_best_particle()->get_personal_best_position();
  
   int ParCounter = 0;
  
   if(fTestLikelihood){
     for(int i = 0; i < fDim; ++i){
       (*PSOParValue)(i) = minimum[i];
       (*PSOParError)(i) = (uncertainties[i][0]+uncertainties[i][1])/2.0;
     }
   }
   else{
     for (auto& ParHandler : systematics)
     {
       if(!ParHandler->IsPCA())
       {
         for(int i = 0; i < ParHandler->GetNumParams(); ++i, ++ParCounter)
         {
           double ParVal = minimum[ParCounter];
           //KS: Basically apply mirroring for parameters out of bounds
           (*PSOParValue)(ParCounter) = ParVal;
           (*PSOParError)(ParCounter) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
           //KS: For fixed params HESS will not calcuate error so we need to pass prior error
           if(ParHandler->IsParameterFixed(i))
           {
             (*PSOParError)(ParCounter) = ParHandler->GetDiagonalError(i);
           }
         }
       }
       else
       {
         //KS: We need to convert parameters from PCA to normal base
         TVectorD ParVals(ParHandler->GetNumParams());
         TVectorD ParVals_PCA(ParHandler->GetNParameters());
  
         TVectorD ErrorVals(ParHandler->GetNumParams());
         TVectorD ErrorVals_PCA(ParHandler->GetNParameters());
  
         //First save them
         //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.
         const int StartVal = ParCounter;
         for(int i = 0; i < ParHandler->GetNParameters(); ++i, ++ParCounter)
         {
           ParVals_PCA(i) = minimum[ParCounter];
           ErrorVals_PCA(i) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
         }
         ParVals = (ParHandler->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
         ErrorVals = (ParHandler->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
  
         ParCounter = StartVal;
         //KS: Now after going from PCA to normal let';s save it
         for(int i = 0; i < ParHandler->GetNumParams(); ++i, ++ParCounter)
         {
           (*PSOParValue)(ParCounter) = ParVals(i);
           (*PSOParError)(ParCounter) = std::fabs(ErrorVals(i));
           //int ParCounterMatrix = StartVal;
           //If fixed take prior
           if(ParHandler->GetPCAHandler()->IsParameterFixedPCA(i))
           {
             (*PSOParError)(ParCounter) = ParHandler->GetDiagonalError(i);
           }
         }
       }
     }
   }
  
   PSOParValue->Write("PSOParValue");
   PSOParError->Write("PSOParError");
   delete PSOParValue;
   delete PSOParError;
   // Save all the output
   SaveOutput();
 }

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 166 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 static const int PSO::kMaxParticles = 10000

staticconstexprprivate

Definition at line 164 of file PSO.h.

◆ paramlist

std::vector<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<std::unique_ptr<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:

Fitters/PSO.h
Fitters/PSO.cpp

Public Member Functions

Private Attributes

Static Private Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

◆ PSO()

◆ ~PSO()

Member Function Documentation

◆ bisection()

◆ calc_uncertainty()

◆ CalcChi()

◆ CalcChi2()

◆ four_vector_addition()

◆ get_best_particle()

◆ init()

◆ rastriginFunc()

◆ run()

◆ RunMCMC()

◆ set_best_particle()

◆ swarmIterate()

◆ three_vector_addition()

◆ uncertainty_check()

◆ vector_add()

◆ vector_multiply()

◆ vector_subtract()

◆ WriteOutput()

Member Data Documentation

◆ best_particle

◆ fBestValue

◆ fConstriction

◆ fConvergence

◆ fDim

◆ fInertia

◆ fIterations

◆ fixed

◆ fOne

◆ fParticles

◆ fTwo

◆ kMaxParticles

◆ paramlist

◆ prior

◆ ranges_max

◆ ranges_min

◆ system

◆ uncertainties

◆ vel