20 fDim = Get<int>(
fitMan->
raw()[
"General"][
"PSO"][
"TestLikelihoodDim"], __FILE__, __LINE__);
34 for(
int i = 0; i <
fDim; ++i){
37 outTree->Branch(Form(
"Parameter_%d", i),
paramlist[i], Form(
"Parameter_%d[nParts]/D",i));
61 for (
int i = 0; i <
fDim; i++){
69 for (std::vector<ParameterHandlerBase*>::iterator it =
systematics.begin(); it !=
systematics.end(); ++it){
72 fDim += (*it)->GetNumParams();
73 for(
int i = 0; i < (*it)->GetNumParams(); ++i)
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);
81 if(high < curr + lim)
ranges_min.push_back(high);
83 prior.push_back(curr);
85 if((*it)->IsParameterFixed(i)){
95 fDim += (*it)->GetNParameters();
96 for(
int i = 0; i < (*it)->GetNParameters(); ++i)
100 prior.push_back((*it)->GetParInit(i));
101 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i)){
112 MACH3LOG_INFO(
"Printing Minimums and Maximums of Variables to be minimized");
113 for (
int i = 0; i <
fDim; i++){
119 std::vector<double> init_position;
120 std::vector<double> init_velocity;
123 for (
int j=0; j<
fDim; ++j){
125 init_position.push_back(
prior[j]);
126 init_velocity.push_back(0.0);
133 init_velocity.push_back((2.0*
random->Rndm()-1.0));
139 double new_value =
CalcChi(init_position);
142 system.push_back(new_particle);
151std::vector<std::vector<double> >
PSO::bisection(
const std::vector<double>& position,
const double minimum,
152 const double range,
const double precision) {
154 std::vector<std::vector<double>> uncertainties_list;
155 for (
unsigned int i = 0; i< position.size(); ++i) {
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;
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};
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]);
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]);
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;
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};
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]);
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]);
227 uncertainties_list.push_back({std::abs(position[i]-position_list[1][i]),std::abs(position[i]-position_list_p[1][i])});
229 MACH3LOG_INFO(
"LLR values for ± positive and negative uncertainties are {:<10.2f} and {:<10.2f}",
232 return uncertainties_list;
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];
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) {
254 double LLR =
CalcChi(position) - minimum - 1.0;
261 int closest_index = 0;
262 double closest_value = std::abs(y[0]);
263 for (
unsigned int ii = 1; ii < y.size(); ++ii) {
264 double abs_y = std::abs(y[ii]);
265 if (abs_y < closest_value) {
267 closest_value = abs_y;
270 neg_uncertainty[i] = x[closest_index];
274 StepPoint = (pos_stop-start) / (num - 1);
276 for (
int j = 0; j < num; ++j) {
278 double LLR =
CalcChi(position) - minimum - 1.0;
285 closest_value = std::abs(y[0]);
286 for (
unsigned int ii = 1; ii < y.size(); ++ii) {
287 double abs_y = std::abs(y[ii]);
288 if (abs_y < closest_value) {
290 closest_value = abs_y;
293 pos_uncertainty[i] = x[closest_index];
295 std::vector<std::vector<double>> res{neg_uncertainty,pos_uncertainty};
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;
315 for (
int j = 0;j < num; ++j) {
317 double LLR =
CalcChi(position);
322 position[i] = curr_ival;
325 for (
unsigned int k = 0;k<x.size(); ++k){
330 for (
unsigned int k = 0; k < x.size(); ++k){
340 std::vector<double> total_pos(
fDim,0.0);
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;});
350 for (
int j = 0; j <
fDim; ++j) {
359 if(
fixed[j]) new_pos[j] =
system[i]->get_position()[j];
364 for (
int j = 0; j <
fDim; ++j) {
366 velo += new_velocity[j]*new_velocity[j];
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);
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;});
389 double mean_dist_sq = 0;
390 for (
int i = 0; i<
fDim; i++){
391 mean_dist_sq += result[i]*result[i];
400 double mean_dist_sq = 0;
422 for (
int j = 0; j <
fDim; ++j){
437 for (
int i = 0; i<
fDim; ++i){
447 TVectorD* PSOParValue =
new TVectorD(
fDim);
448 TVectorD* PSOParError =
new TVectorD(
fDim);
450 for(
int i = 0; i <
fDim; ++i)
452 (*PSOParValue)(i) = 0;
453 (*PSOParError)(i) = 0;
461 for(
int i = 0; i <
fDim; ++i){
462 (*PSOParValue)(i) = minimum[i];
467 for (std::vector<ParameterHandlerBase*>::iterator it =
systematics.begin(); it !=
systematics.end(); ++it)
471 for(
int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
473 double ParVal = minimum[ParCounter];
475 (*PSOParValue)(ParCounter) = ParVal;
478 if((*it)->IsParameterFixed(i))
480 (*PSOParError)(ParCounter) = (*it)->GetDiagonalError(i);
487 TVectorD ParVals((*it)->GetNumParams());
488 TVectorD ParVals_PCA((*it)->GetNParameters());
490 TVectorD ErrorVals((*it)->GetNumParams());
491 TVectorD ErrorVals_PCA((*it)->GetNParameters());
495 const int StartVal = ParCounter;
496 for(
int i = 0; i < (*it)->GetNParameters(); ++i, ++ParCounter)
498 ParVals_PCA(i) = minimum[ParCounter];
501 ParVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ParVals_PCA;
502 ErrorVals = ((*it)->GetPCAHandler()->GetTransferMatrix())*ErrorVals_PCA;
504 ParCounter = StartVal;
506 for(
int i = 0; i < (*it)->GetNumParams(); ++i, ++ParCounter)
508 (*PSOParValue)(ParCounter) = ParVals(i);
509 (*PSOParError)(ParCounter) = std::fabs(ErrorVals(i));
512 if((*it)->GetPCAHandler()->IsParameterFixedPCA(i))
514 (*PSOParError)(ParCounter) = (*it)->GetDiagonalError(i);
521 PSOParValue->Write(
"PSOParValue");
522 PSOParError->Write(
"PSOParError");
545 constexpr double A = 10.0;
547 for (
int i = 0; i <
fDim; ++i) {
548 sum += x[i] * x[i] - A * cos(2.0 * 3.14 * x[i]);
550 double llh = A *
fDim + sum;
MaCh3Plotting::PlottingManager * man
std::unique_ptr< TRandom3 > random
Random number.
int accCount
counts accepted steps
void SaveOutput()
Save output and close files.
TFile * outputFile
Output.
manager * fitMan
The manager.
unsigned int step
current state
double accProb
current acceptance prob
double stepTime
Time of single step.
std::unique_ptr< TStopwatch > stepClock
tells how long single step/fit iteration took
double logLCurr
current likelihood
int auto_save
auto save every N steps
bool fTestLikelihood
Necessary for some fitting algorithms like PSO.
TTree * outTree
Output tree with posteriors.
void SanitiseInputs()
Remove obsolete memory and make other checks before fit starts.
std::vector< ParameterHandlerBase * > systematics
Systematic holder.
Implementation of base Likelihood Fit class, it is mostly responsible for likelihood calculation whil...
void PrepareFit()
prepare output and perform sanity checks
virtual double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
void RunMCMC() override
Actual implementation of PSO Fit algorithm.
PSO(manager *const fitMan)
constructor
std::vector< double > ranges_max
std::vector< double > three_vector_addition(std::vector< double > vec1, const std::vector< double > &vec2, const std::vector< double > &vec3)
double CalcChi2(const double *x)
Chi2 calculation over all included samples and syst objects.
std::vector< std::vector< double > > bisection(const std::vector< double > &position, const double minimum, const double range, const double precision)
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 > vector_multiply(std::vector< double > velocity, const double mul)
std::vector< bool > fixed
particle * get_best_particle()
std::vector< std::vector< double > > calc_uncertainty(const std::vector< double > &position, const double minimum)
void set_best_particle(particle *n)
std::vector< double * > paramlist
std::vector< particle * > system
double vel[kMaxParticles]
double CalcChi(std::vector< double > x)
std::vector< std::vector< double > > uncertainties
std::vector< double > ranges_min
double rastriginFunc(const double *x)
Evaluates the Rastrigin function for a given parameter values.
std::vector< double > prior
void uncertainty_check(const std::vector< double > &previous_pos)
The manager class is responsible for managing configurations and settings.
YAML::Node const & raw()
Return config.
Class particle - stores the position, velocity and personal best With functions which move particle a...
void set_personal_best_value(const double new_val)
void set_personal_best_position(const std::vector< double > &new_pos)
std::vector< double > get_personal_best_position()
void set_value(const double val)
static constexpr const double _LARGE_LOGL_
Large Likelihood is used it parameter go out of physical boundary, this indicates in MCMC that such s...