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.

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 () 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
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< 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
manager *	fitMan
	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.

                     : LikelihoodFit(man) {
// ***************
  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 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.

                                                                                             {
// *************************
  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 236 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 139 of file PSO.h.

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

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

                                    {
// *******************
  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 131 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 84 of file PSO.h.

                                 {
        return best_particle;
    }

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

              {
// *************************
  fBestValue = M3::_LARGE_LOGL_;
 
  //KS: For none PCA this will be eqaul 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)->GetParInit(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)->GetParInit(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);
      }
    }
 
    particle* new_particle = new particle(init_position,init_velocity);
    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);
    system.push_back(new_particle);
    if(new_value < fBestValue){
      fBestValue = new_value;
      set_best_particle(new_particle);
    }
  }
}

◆ rastriginFunc()

double PSO::rastriginFunc ( const double * x )

Evaluates the Rastrigin function for a given parameter values.

Definition at line 540 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 398 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 ( )

overridevirtual

Actual implementation of PSO Fit algorithm.

Implements LikelihoodFit.

Definition at line 25 of file PSO.cpp.

                 {
// ***************
  PrepareFit();
 
  // Remove obsolete memory and make other checks before fit starts
  SanitiseInputs();
 
  if(fTestLikelihood){
    outTree->Branch("nParts", &fParticles, "nParts/I");
    for(int i = 0; i < fDim; ++i){
      par = new double[fParticles];
      paramlist.push_back(par);
      outTree->Branch(Form("Parameter_%d", i), paramlist[i], 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 87 of file PSO.h.

                                       {
      best_particle = n;
    }

◆ swarmIterate()

double PSO::swarmIterate ( )

Definition at line 338 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]);
      }
    }
  }
 
  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(), [=](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 123 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 300 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 113 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 105 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 118 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 443 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 (std::vector<ParameterHandlerBase*>::iterator it = systematics.begin(); it != systematics.end(); ++it)
    {
      if(!(*it)->IsPCA())
      {
        for(int i = 0; i < (*it)->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((*it)->IsParameterFixed(i))
          {
            (*PSOParError)(ParCounter) = (*it)->GetDiagonalError(i);
          }
        }
      }
      else
      {
        //KS: We need to convert parameters from PCA to normal base
        TVectorD ParVals((*it)->GetNumParams());
        TVectorD ParVals_PCA((*it)->GetNParameters());
 
        TVectorD ErrorVals((*it)->GetNumParams());
        TVectorD ErrorVals_PCA((*it)->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 < (*it)->GetNParameters(); ++i, ++ParCounter)
        {
          ParVals_PCA(i) = minimum[ParCounter];
          ErrorVals_PCA(i) = (uncertainties[ParCounter][0]+uncertainties[ParCounter][1])/2.0;
        }
        ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
        ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
 
        ParCounter = StartVal;
        //KS: Now after going from PCA to normal let';s save it
        for(int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
        {
          (*PSOParValue)(ParCounter) = ParVals(i);
          (*PSOParError)(ParCounter) = std::fabs(ErrorVals(i));
          //int ParCounterMatrix = StartVal;
          //If fixed take prior
          if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
          {
            (*PSOParError)(ParCounter) = (*it)->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 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:

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()

◆ GetName()

◆ 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

◆ par

◆ paramlist

◆ prior

◆ ranges_max

◆ ranges_min

◆ system

◆ uncertainties

◆ vel