00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "config.h"
00011
00012 #include "BAT/BCIntegrate.h"
00013 #include "BAT/BCLog.h"
00014 #include "BAT/BCMath.h"
00015
00016 #include <TH1D.h>
00017 #include <TH2D.h>
00018 #include <TMinuit.h>
00019 #include <TRandom3.h>
00020 #include <TTree.h>
00021
00022 #ifdef HAVE_CUBA_H
00023 #include "cuba.h"
00024 #endif
00025
00026
00027
00028 class BCIntegrate * global_this;
00029
00030
00031 BCIntegrate::BCIntegrate() : BCEngineMCMC()
00032 {
00033 fNvar=0;
00034 fNiterPerDimension = 100;
00035 fNSamplesPer2DBin = 100;
00036 fRandom = new TRandom3(0);
00037
00038 fNIterationsMax = 1000000;
00039 fRelativePrecision = 1e-3;
00040
00041 #ifdef HAVE_CUBA_H
00042 fIntegrationMethod = BCIntegrate::kIntCuba;
00043 #else
00044 fIntegrationMethod = BCIntegrate::kIntMonteCarlo;
00045 #endif
00046 fMarginalizationMethod = BCIntegrate::kMargMetropolis;
00047 fOptimizationMethod = BCIntegrate::kOptMinuit;
00048 fOptimizationMethodMode = BCIntegrate::kOptMinuit;
00049
00050 fNbins=100;
00051
00052 fDataPointLowerBoundaries = 0;
00053 fDataPointUpperBoundaries = 0;
00054
00055 fFillErrorBand = false;
00056
00057 fFitFunctionIndexX = -1;
00058 fFitFunctionIndexY = -1;
00059
00060 fErrorBandNbinsX = 100;
00061 fErrorBandNbinsY = 500;
00062
00063 fErrorBandContinuous = true;
00064
00065 fMinuit = 0;
00066 fMinuitArglist[0] = 20000;
00067 fMinuitArglist[1] = 0.01;
00068 fFlagIgnorePrevOptimization = false;
00069
00070 fFlagWriteMarkovChain = false;
00071 fMarkovChainTree = 0;
00072
00073 fMarkovChainStepSize = 0.1;
00074
00075 fMarkovChainAutoN = true;
00076
00077 fSAT0 = 100.0;
00078 fSATmin = 0.1;
00079 fSASchedule = BCIntegrate::kSACauchy;
00080 fFlagWriteSAToFile = false;
00081 }
00082
00083
00084 BCIntegrate::BCIntegrate(BCParameterSet * par) : BCEngineMCMC(1)
00085 {
00086 fNvar=0;
00087 fNiterPerDimension = 100;
00088 fNSamplesPer2DBin = 100;
00089 fRandom = new TRandom3(0);
00090
00091 fNbins=100;
00092
00093 this->SetParameters(par);
00094
00095 fDataPointLowerBoundaries = 0;
00096 fDataPointUpperBoundaries = 0;
00097
00098 fFillErrorBand = false;
00099
00100 fFitFunctionIndexX = -1;
00101 fFitFunctionIndexY = -1;
00102
00103 fErrorBandNbinsX = 100;
00104 fErrorBandNbinsY = 200;
00105
00106 fErrorBandContinuous = true;
00107
00108 fMinuit = 0;
00109 fMinuitArglist[0] = 20000;
00110 fMinuitArglist[1] = 0.01;
00111 fFlagIgnorePrevOptimization = false;
00112
00113 fFlagWriteMarkovChain = false;
00114 fMarkovChainTree = 0;
00115
00116 fMarkovChainStepSize = 0.1;
00117
00118 fMarkovChainAutoN = true;
00119
00120 fSAT0 = 100.0;
00121 fSATmin = 0.1;
00122 fTreeSA = false;
00123 }
00124
00125
00126 BCIntegrate::~BCIntegrate()
00127 {
00128 DeleteVarList();
00129
00130 fx=0;
00131
00132 delete fRandom;
00133 fRandom=0;
00134
00135 if (fMinuit)
00136 delete fMinuit;
00137
00138 int n1 = fHProb1D.size();
00139 if(n1>0)
00140 {
00141 for (int i=0;i<n1;i++)
00142 delete fHProb1D.at(i);
00143 }
00144
00145 int n2 = fHProb2D.size();
00146 if(n2>0)
00147 {
00148 for (int i=0;i<n2;i++)
00149 delete fHProb2D.at(i);
00150 }
00151 }
00152
00153
00154 void BCIntegrate::SetParameters(BCParameterSet * par)
00155 {
00156 DeleteVarList();
00157
00158 fx = par;
00159 fNvar = fx->size();
00160 fMin = new double[fNvar];
00161 fMax = new double[fNvar];
00162 fVarlist = new int[fNvar];
00163
00164 this->ResetVarlist(1);
00165
00166 for(int i=0;i<fNvar;i++)
00167 {
00168 fMin[i]=(fx->at(i))->GetLowerLimit();
00169 fMax[i]=(fx->at(i))->GetUpperLimit();
00170 }
00171
00172 fXmetro0.clear();
00173 for(int i=0;i<fNvar;i++)
00174 fXmetro0.push_back((fMin[i]+fMax[i])/2.0);
00175
00176 fXmetro1.clear();
00177 fXmetro1.assign(fNvar,0.);
00178
00179 fMCMCBoundaryMin.clear();
00180 fMCMCBoundaryMax.clear();
00181
00182 for(int i=0;i<fNvar;i++)
00183 {
00184 fMCMCBoundaryMin.push_back(fMin[i]);
00185 fMCMCBoundaryMax.push_back(fMax[i]);
00186 fMCMCFlagsFillHistograms.push_back(true);
00187 }
00188
00189 for (int i = int(fMCMCH1NBins.size()); i<fNvar; ++i)
00190 fMCMCH1NBins.push_back(100);
00191
00192 fMCMCNParameters = fNvar;
00193 }
00194
00195
00196 void BCIntegrate::SetMarkovChainInitialPosition(std::vector<double> position)
00197 {
00198 int n = position.size();
00199
00200 fXmetro0.clear();
00201
00202 for (int i = 0; i < n; ++i)
00203 fXmetro0.push_back(position.at(i));
00204 }
00205
00206
00207 void BCIntegrate::DeleteVarList()
00208 {
00209 if(fNvar)
00210 {
00211 delete[] fVarlist;
00212 fVarlist=0;
00213
00214 delete[] fMin;
00215 fMin=0;
00216
00217 delete[] fMax;
00218 fMax=0;
00219
00220 fx=0;
00221 fNvar=0;
00222 }
00223 }
00224
00225
00226 void BCIntegrate::SetNbins(int nbins, int index)
00227 {
00228 if (fNvar == 0)
00229 return;
00230
00231
00232 if (index >= fNvar)
00233 {
00234 BCLog::Out(BCLog::warning, BCLog::warning,"BCIntegrate::SetNbins : Index out of range.");
00235 return;
00236 }
00237
00238 else if (index < 0)
00239 {
00240 for (int i = 0; i < fNvar; ++i)
00241 this -> SetNbins(nbins, i);
00242 return;
00243 }
00244
00245
00246 if (nbins <= 0)
00247 nbins = 10;
00248
00249 fMCMCH1NBins[index] = nbins;
00250
00251 return;
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265 }
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288
00289
00290
00291
00292 void BCIntegrate::SetVarList(int * varlist)
00293 {
00294 for(int i=0;i<fNvar;i++)
00295 fVarlist[i]=varlist[i];
00296 }
00297
00298
00299 void BCIntegrate::ResetVarlist(int v)
00300 {
00301 for(int i=0;i<fNvar;i++)
00302 fVarlist[i]=v;
00303 }
00304
00305
00306 double BCIntegrate::Eval(std::vector <double> x)
00307 {
00308 BCLog::Out(BCLog::warning, BCLog::warning, "BCIntegrate::Eval. No function. Function needs to be overloaded.");
00309 return 0;
00310 }
00311
00312
00313 double BCIntegrate::LogEval(std::vector <double> x)
00314 {
00315
00316 return log(this->Eval(x));
00317 }
00318
00319
00320 double BCIntegrate::EvalSampling(std::vector <double> x)
00321 {
00322 BCLog::Out(BCLog::warning, BCLog::warning, "BCIntegrate::EvalSampling. No function. Function needs to be overloaded.");
00323 return 0;
00324 }
00325
00326
00327 double BCIntegrate::LogEvalSampling(std::vector <double> x)
00328 {
00329 return log(this->EvalSampling(x));
00330 }
00331
00332
00333 double BCIntegrate::EvalPrint(std::vector <double> x)
00334 {
00335 double val=this->Eval(x);
00336 BCLog::Out(BCLog::detail, BCLog::detail, Form("BCIntegrate::EvalPrint. Value: %d.", val));
00337
00338 return val;
00339 }
00340
00341
00342 void BCIntegrate::SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method)
00343 {
00344 #ifdef HAVE_CUBA_H
00345 fIntegrationMethod = method;
00346 #else
00347 BCLog::Out(BCLog::warning,BCLog::warning,"!!! This version of BAT is compiled without Cuba.");
00348 BCLog::Out(BCLog::warning,BCLog::warning," Monte Carlo Sampled Mean integration method will be used.");
00349 BCLog::Out(BCLog::warning,BCLog::warning," To be able to use Cuba integration, install Cuba and recompile BAT.");
00350 #endif
00351 }
00352
00353
00354 double BCIntegrate::Integrate()
00355 {
00356 std::vector <double> parameter;
00357 parameter.assign(fNvar, 0.0);
00358
00359 switch(fIntegrationMethod)
00360 {
00361 case BCIntegrate::kIntMonteCarlo:
00362 return IntegralMC(parameter);
00363
00364 case BCIntegrate::kIntMetropolis:
00365 return this -> IntegralMetro(parameter);
00366
00367 case BCIntegrate::kIntImportance:
00368 return this -> IntegralImportance(parameter);
00369
00370 case BCIntegrate::kIntCuba:
00371 #ifdef HAVE_CUBA_H
00372 return this -> CubaIntegrate();
00373 #else
00374 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
00375 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT.");
00376 #endif
00377 }
00378
00379 BCLog::Out(BCLog::error, BCLog::error,
00380 Form("BCIntegrate::Integrate : Invalid integration method: %d", fIntegrationMethod));
00381
00382 return 0;
00383 }
00384
00385
00386 double BCIntegrate::IntegralMC(std::vector <double> x, int * varlist)
00387 {
00388 this->SetVarList(varlist);
00389 return IntegralMC(x);
00390 }
00391
00392
00393 double BCIntegrate::IntegralMC(std::vector <double> x)
00394 {
00395
00396 int NvarNow=0;
00397
00398 for(int i = 0; i < fNvar; i++)
00399 if(fVarlist[i])
00400 NvarNow++;
00401
00402
00403 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Running MC integation over %i dimensions.", NvarNow));
00404 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Maximum number of iterations: %i", this->GetNIterationsMax()));
00405 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Aimed relative precision: %e", this->GetRelativePrecision()));
00406
00407
00408 double D = 1.0;
00409 for(int i = 0; i < fNvar; i++)
00410 if (fVarlist[i])
00411 D = D * (fMax[i] - fMin[i]);
00412
00413
00414 double pmax = 0.0;
00415 double sumW = 0.0;
00416 double sumW2 = 0.0;
00417 double integral = 0.0;
00418 double variance = 0.0;
00419 double relprecision = 1.0;
00420
00421 std::vector <double> randx;
00422 randx.assign(fNvar, 0.0);
00423
00424
00425 fNIterations = 0;
00426
00427
00428 while ((fRelativePrecision < relprecision && fNIterationsMax > fNIterations) || fNIterations < 10)
00429 {
00430
00431 fNIterations++;
00432
00433
00434 this -> GetRandomVector(randx);
00435
00436
00437 for(int i = 0; i < fNvar; i++)
00438 {
00439 if(fVarlist[i])
00440 randx[i]=fMin[i]+randx[i]*(fMax[i]-fMin[i]);
00441 else
00442 randx[i]=x[i];
00443 }
00444
00445
00446 double value = this->Eval(randx);
00447
00448
00449 sumW += value;
00450 sumW2 += value * value;
00451
00452
00453 if (value > pmax)
00454 {
00455
00456 pmax = value;
00457
00458
00459 fBestFitParameters.clear();
00460
00461
00462 for(int i = 0; i < fNvar; i++)
00463 fBestFitParameters.push_back(randx.at(i));
00464 }
00465
00466
00467 integral = D * sumW / fNIterations;
00468
00469 if (fNIterations%10000 == 0)
00470 {
00471 variance = (1.0 / double(fNIterations)) * (D * D * sumW2 / double(fNIterations) - integral * integral);
00472 double error = sqrt(variance);
00473 relprecision = error / integral;
00474
00475 BCLog::Out(BCLog::debug, BCLog::debug,
00476 Form("BCIntegrate::IntegralMC. Iteration %i, integral: %e +- %e.", fNIterations, integral, error));
00477 }
00478 }
00479
00480 fError = variance / fNIterations;
00481
00482
00483 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Result of integration: %e +- %e in %i iterations.", integral, sqrt(variance), fNIterations));
00484 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Obtained relative precision: %e. ", sqrt(variance) / integral));
00485
00486 return integral;
00487 }
00488
00489
00490
00491 double BCIntegrate::IntegralMetro(std::vector <double> x)
00492 {
00493
00494 BCLog::Out(BCLog::debug, BCLog::debug, Form("BCIntegrate::IntegralMetro. Integate over %i dimensions.", fNvar));
00495
00496
00497 double Niter = pow(fNiterPerDimension, fNvar);
00498
00499
00500 if(Niter>1e5)
00501 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Total number of iterations in Metropolis: %d.", Niter));
00502
00503
00504 double sumI = 0;
00505
00506
00507 std::vector <double> randx;
00508 randx.assign(fNvar,0.);
00509 InitMetro();
00510
00511
00512 double pmax = 0.0;
00513
00514
00515 for(int i = 0; i <= Niter; i++)
00516 {
00517
00518 this -> GetRandomPointSamplingMetro(randx);
00519
00520
00521 double val_f = this -> Eval(randx);
00522
00523
00524 double val_g = this -> EvalSampling(randx);
00525
00526
00527 if (val_g > 0)
00528 sumI += val_f / val_g;
00529
00530
00531 if (val_f > pmax)
00532 {
00533
00534 pmax = val_f;
00535
00536
00537 fBestFitParameters.clear();
00538
00539
00540 for(int i = 0; i < fNvar; i++)
00541 fBestFitParameters.push_back(randx.at(i));
00542 }
00543
00544
00545 if((i+1)%100000 == 0)
00546 BCLog::Out(BCLog::debug, BCLog::debug,
00547 Form("BCIntegrate::IntegralMetro. Iteration %d, integral: %d.", i+1, sumI/(i+1)));
00548 }
00549
00550
00551 double result=sumI/Niter;
00552
00553
00554 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Integral is %f in %i iterations.", result, Niter));
00555
00556 return result;
00557 }
00558
00559
00560 double BCIntegrate::IntegralImportance(std::vector <double> x)
00561 {
00562
00563 BCLog::Out(BCLog::debug, BCLog::debug, Form("BCIntegrate::IntegralImportance. Integate over %i dimensions.", fNvar));
00564
00565
00566 double Niter = pow(fNiterPerDimension, fNvar);
00567
00568
00569 if(Niter>1e5)
00570 BCLog::Out(BCLog::detail, BCLog::detail, Form("BCIntegrate::IntegralImportance. Total number of iterations: %d.", Niter));
00571
00572
00573 double sumI = 0;
00574
00575 std::vector <double> randx;
00576 randx.assign(fNvar,0.);
00577
00578
00579 double pmax = 0.0;
00580
00581
00582 for(int i = 0; i <= Niter; i++)
00583 {
00584
00585 this -> GetRandomPointImportance(randx);
00586
00587
00588 double val_f = this -> Eval(randx);
00589
00590
00591 double val_g = this -> EvalSampling(randx);
00592
00593
00594 if (val_g > 0)
00595 sumI += val_f / val_g;
00596
00597
00598 if (val_f > pmax)
00599 {
00600
00601 pmax = val_f;
00602
00603
00604 fBestFitParameters.clear();
00605
00606
00607 for(int i = 0; i < fNvar; i++)
00608 fBestFitParameters.push_back(randx.at(i));
00609 }
00610
00611
00612 if((i+1)%100000 == 0)
00613 BCLog::Out(BCLog::debug, BCLog::debug,
00614 Form("BCIntegrate::IntegralImportance. Iteration %d, integral: %d.", i+1, sumI/(i+1)));
00615 }
00616
00617
00618 double result=sumI/Niter;
00619
00620
00621 BCLog::Out(BCLog::debug, BCLog::debug, Form("BCIntegrate::IntegralImportance. Integral %f in %i iterations.", result, Niter));
00622
00623 return result;
00624 }
00625
00626
00627 TH1D* BCIntegrate::Marginalize(BCParameter * parameter)
00628 {
00629 BCLog::Out(BCLog::detail, BCLog::detail,
00630 Form(" --> Marginalizing model wrt. parameter %s using method %d.", parameter->GetName().data(), fMarginalizationMethod));
00631
00632 switch(fMarginalizationMethod)
00633 {
00634 case BCIntegrate::kMargMonteCarlo:
00635 return MarginalizeByIntegrate(parameter);
00636
00637 case BCIntegrate::kMargMetropolis:
00638 return MarginalizeByMetro(parameter);
00639 }
00640
00641 BCLog::Out(BCLog::warning, BCLog::warning,
00642 Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod));
00643
00644 return 0;
00645 }
00646
00647
00648 TH2D * BCIntegrate::Marginalize(BCParameter * parameter1, BCParameter * parameter2)
00649 {
00650 switch(fMarginalizationMethod)
00651 {
00652 case BCIntegrate::kMargMonteCarlo:
00653 return MarginalizeByIntegrate(parameter1,parameter2);
00654
00655 case BCIntegrate::kMargMetropolis:
00656 return MarginalizeByMetro(parameter1,parameter2);
00657 }
00658
00659 BCLog::Out(BCLog::warning, BCLog::warning,
00660 Form("BCIntegrate::Marginalize. Invalid marginalization method: %d. Return 0.", fMarginalizationMethod));
00661
00662 return 0;
00663 }
00664
00665
00666 TH1D* BCIntegrate::MarginalizeByIntegrate(BCParameter * parameter)
00667 {
00668
00669 this->ResetVarlist(1);
00670 int index = parameter->GetIndex();
00671 this->UnsetVar(index);
00672
00673
00674 double hmin = parameter -> GetLowerLimit();
00675 double hmax = parameter -> GetUpperLimit();
00676 double hdx = (hmax - hmin) / double(fNbins);
00677 TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax);
00678
00679
00680 std::vector <double> randx;
00681 randx.assign(fNvar, 0.0);
00682
00683 for(int i=0;i<=fNbins;i++)
00684 {
00685 randx[index] = hmin + (double)i * hdx;
00686
00687 double val = IntegralMC(randx);
00688 hist->Fill(randx[index], val);
00689 }
00690
00691
00692 hist -> Scale( 1./hist->Integral("width") );
00693
00694 return hist;
00695 }
00696
00697
00698 TH2D * BCIntegrate::MarginalizeByIntegrate(BCParameter * parameter1, BCParameter * parameter2)
00699 {
00700
00701 this->ResetVarlist(1);
00702 int index1 = parameter1->GetIndex();
00703 this->UnsetVar(index1);
00704 int index2 = parameter2->GetIndex();
00705 this->UnsetVar(index2);
00706
00707
00708 double hmin1 = parameter1 -> GetLowerLimit();
00709 double hmax1 = parameter1 -> GetUpperLimit();
00710 double hdx1 = (hmax1 - hmin1) / double(fNbins);
00711
00712 double hmin2 = parameter2 -> GetLowerLimit();
00713 double hmax2 = parameter2 -> GetUpperLimit();
00714 double hdx2 = (hmax2 - hmin2) / double(fNbins);
00715
00716 TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1 -> GetName().data(), parameter2 -> GetName().data()),"",
00717 fNbins, hmin1, hmax1,
00718 fNbins, hmin2, hmax2);
00719
00720
00721 std::vector <double> randx;
00722 randx.assign(fNvar, 0.0);
00723
00724 for(int i=0;i<=fNbins;i++)
00725 {
00726 randx[index1] = hmin1 + (double)i * hdx1;
00727 for(int j=0;j<=fNbins;j++)
00728 {
00729 randx[index2] = hmin2 + (double)j * hdx2;
00730
00731 double val = IntegralMC(randx);
00732 hist->Fill(randx[index1],randx[index2], val);
00733 }
00734 }
00735
00736
00737 hist -> Scale(1.0/hist->Integral("width"));
00738
00739 return hist;
00740 }
00741
00742
00743 TH1D * BCIntegrate::MarginalizeByMetro(BCParameter * parameter)
00744 {
00745 int niter = fMarkovChainNIterations;
00746
00747 if (fMarkovChainAutoN == true)
00748 niter = fNbins*fNbins*fNSamplesPer2DBin*fNvar;
00749
00750 BCLog::Out(BCLog::detail, BCLog::detail,
00751 Form(" --> Number of samples in Metropolis marginalization: %d.", niter));
00752
00753
00754 int index = parameter->GetIndex();
00755
00756
00757 double hmin = parameter -> GetLowerLimit();
00758 double hmax = parameter -> GetUpperLimit();
00759 TH1D * hist = new TH1D("hist","", fNbins, hmin, hmax);
00760
00761
00762 std::vector <double> randx;
00763 randx.assign(fNvar, 0.0);
00764 InitMetro();
00765
00766 for(int i=0;i<=niter;i++)
00767 {
00768 GetRandomPointMetro(randx);
00769 hist->Fill(randx[index]);
00770 }
00771
00772
00773 hist -> Scale(1.0/hist->Integral("width"));
00774
00775 return hist;
00776 }
00777
00778
00779 TH2D * BCIntegrate::MarginalizeByMetro(BCParameter * parameter1, BCParameter * parameter2)
00780 {
00781 int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar;
00782
00783
00784 int index1 = parameter1->GetIndex();
00785 int index2 = parameter2->GetIndex();
00786
00787
00788 double hmin1 = parameter1 -> GetLowerLimit();
00789 double hmax1 = parameter1 -> GetUpperLimit();
00790
00791 double hmin2 = parameter2 -> GetLowerLimit();
00792 double hmax2 = parameter2 -> GetUpperLimit();
00793
00794 TH2D * hist = new TH2D(Form("hist_%s_%s", parameter1 -> GetName().data(), parameter2 -> GetName().data()),"",
00795 fNbins, hmin1, hmax1,
00796 fNbins, hmin2, hmax2);
00797
00798
00799 std::vector <double> randx;
00800 randx.assign(fNvar, 0.0);
00801 InitMetro();
00802
00803 for(int i=0;i<=niter;i++)
00804 {
00805 GetRandomPointMetro(randx);
00806 hist->Fill(randx[index1],randx[index2]);
00807 }
00808
00809
00810 hist -> Scale(1.0/hist->Integral("width"));
00811
00812 return hist;
00813 }
00814
00815
00816 int BCIntegrate::MarginalizeAllByMetro(const char * name="")
00817 {
00818 int niter=fNbins*fNbins*fNSamplesPer2DBin*fNvar;
00819
00820 BCLog::Out(BCLog::detail, BCLog::detail,
00821 Form(" --> Number of samples in Metropolis marginalization: %d.", niter));
00822
00823
00824 for(int i=0;i<fNvar;i++)
00825 {
00826 double hmin1 = fx->at(i) -> GetLowerLimit();
00827 double hmax1 = fx->at(i) -> GetUpperLimit();
00828
00829 TH1D * h1 = new TH1D(
00830 TString::Format("h%s_%s", name, fx->at(i) -> GetName().data()),"",
00831 fNbins, hmin1, hmax1);
00832
00833 fHProb1D.push_back(h1);
00834 }
00835
00836
00837 for(int i=0;i<fNvar-1;i++)
00838 for(int j=i+1;j<fNvar;j++)
00839 {
00840 double hmin1 = fx->at(i) -> GetLowerLimit();
00841 double hmax1 = fx->at(i) -> GetUpperLimit();
00842
00843 double hmin2 = fx->at(j) -> GetLowerLimit();
00844 double hmax2 = fx->at(j) -> GetUpperLimit();
00845
00846 TH2D * h2 = new TH2D(
00847 TString::Format("h%s_%s_%s", name, fx->at(i) -> GetName().data(), fx->at(j) -> GetName().data()),"",
00848 fNbins, hmin1, hmax1,
00849 fNbins, hmin2, hmax2);
00850
00851 fHProb2D.push_back(h2);
00852 }
00853
00854
00855 int nh2d = fHProb2D.size();
00856
00857 BCLog::Out(BCLog::detail, BCLog::detail,
00858 Form(" --> Marginalizing %d 1D distributions and %d 2D distributions.", fNvar, nh2d));
00859
00860
00861 double dx = 0.0;
00862 double dy = 0.0;
00863
00864 if (fFitFunctionIndexX >= 0)
00865 {
00866 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) -
00867 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX))
00868 / double(fErrorBandNbinsX);
00869
00870 dx = (fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) -
00871 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY))
00872 / double(fErrorBandNbinsY);
00873
00874 fErrorBandXY = new TH2D(
00875 TString::Format("errorbandxy_%d",BCLog::GetHIndex()), "",
00876 fErrorBandNbinsX,
00877 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexX) - 0.5 * dx,
00878 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexX) + 0.5 * dx,
00879 fErrorBandNbinsY,
00880 fDataPointLowerBoundaries -> GetValue(fFitFunctionIndexY) - 0.5 * dy,
00881 fDataPointUpperBoundaries -> GetValue(fFitFunctionIndexY) + 0.5 * dy);
00882 fErrorBandXY -> SetStats(kFALSE);
00883
00884 for (int ix = 1; ix <= fErrorBandNbinsX; ++ix)
00885 for (int iy = 1; iy <= fErrorBandNbinsX; ++iy)
00886 fErrorBandXY -> SetBinContent(ix, iy, 0.0);
00887 }
00888
00889
00890 std::vector <double> randx;
00891
00892 randx.assign(fNvar, 0.0);
00893 InitMetro();
00894
00895
00896 fNacceptedMCMC = 0;
00897
00898
00899 for(int i=0;i<=niter;i++)
00900 {
00901 GetRandomPointMetro(randx);
00902
00903
00904 if (fFlagWriteMarkovChain)
00905 {
00906 fMCMCIteration = i;
00907 fMarkovChainTree -> Fill();
00908 }
00909
00910 for(int j=0;j<fNvar;j++)
00911 fHProb1D[j] -> Fill( randx[j] );
00912
00913 int ih=0;
00914 for(int j=0;j<fNvar-1;j++)
00915 for(int k=j+1;k<fNvar;k++)
00916 {
00917 fHProb2D[ih] -> Fill(randx[j],randx[k]);
00918 ih++;
00919 }
00920
00921 if((i+1)%100000==0)
00922 BCLog::Out(BCLog::debug, BCLog::debug,
00923 Form("BCIntegrate::MarginalizeAllByMetro. %d samples done.",i+1));
00924
00925
00926
00927 if (fFitFunctionIndexX >= 0)
00928 {
00929
00930 if (fErrorBandContinuous)
00931 {
00932 double x = 0;
00933
00934 for (int ix = 0; ix < fErrorBandNbinsX; ix++)
00935 {
00936
00937 x = fErrorBandXY -> GetXaxis() -> GetBinCenter(ix + 1);
00938
00939
00940 std::vector <double> xvec;
00941 xvec.push_back(x);
00942 double y = this -> FitFunction(xvec, randx);
00943 xvec.clear();
00944
00945
00946 fErrorBandXY -> Fill(x, y);
00947 }
00948 }
00949
00950
00951 else
00952 {
00953 int ndatapoints = int(fErrorBandX.size());
00954 double x = 0;
00955
00956 for (int ix = 0; ix < ndatapoints; ++ix)
00957 {
00958
00959 x = fErrorBandX.at(ix);
00960
00961
00962 std::vector <double> xvec;
00963 xvec.push_back(x);
00964 double y = this -> FitFunction(xvec, randx);
00965 xvec.clear();
00966
00967
00968 fErrorBandXY -> Fill(x, y);
00969 }
00970 }
00971 }
00972 }
00973
00974
00975 for(int i=0;i<fNvar;i++)
00976 fHProb1D[i] -> Scale( 1./fHProb1D[i]->Integral("width") );
00977
00978 for (int i=0;i<nh2d;i++)
00979 fHProb2D[i] -> Scale( 1./fHProb2D[i]->Integral("width") );
00980
00981 if (fFitFunctionIndexX >= 0)
00982 fErrorBandXY -> Scale(1.0/fErrorBandXY -> Integral() * fErrorBandXY -> GetNbinsX());
00983
00984 if (fFitFunctionIndexX >= 0)
00985 {
00986 for (int ix = 1; ix <= fErrorBandNbinsX; ix++)
00987 {
00988 double sum = 0;
00989
00990 for (int iy = 1; iy <= fErrorBandNbinsY; iy++)
00991 sum += fErrorBandXY -> GetBinContent(ix, iy);
00992
00993 for (int iy = 1; iy <= fErrorBandNbinsY; iy++)
00994 {
00995 double newvalue = fErrorBandXY -> GetBinContent(ix, iy) / sum;
00996 fErrorBandXY -> SetBinContent(ix, iy, newvalue);
00997 }
00998 }
00999 }
01000
01001 BCLog::Out(BCLog::detail, BCLog::detail,
01002 Form("BCIntegrate::MarginalizeAllByMetro done with %i trials and %i accepted trials. Efficiency is %f",fNmetro, fNacceptedMCMC, double(fNacceptedMCMC)/double(fNmetro)));
01003
01004 return fNvar+nh2d;
01005 }
01006
01007
01008 TH1D * BCIntegrate::GetH1D(int parIndex)
01009 {
01010 if(fHProb1D.size()==0)
01011 {
01012 BCLog::Out(BCLog::warning, BCLog::warning,
01013 "BCModel::GetH1D. MarginalizeAll() has to be run prior to this to fill the distributions.");
01014 return 0;
01015 }
01016
01017 if(parIndex<0 || parIndex>fNvar)
01018 {
01019 BCLog::Out(BCLog::warning, BCLog::warning,
01020 Form("BCIntegrate::GetH1D. Parameter index %d is invalid.",parIndex));
01021 return 0;
01022 }
01023
01024 return fHProb1D[parIndex];
01025 }
01026
01027
01028 int BCIntegrate::GetH2DIndex(int parIndex1, int parIndex2)
01029 {
01030 if(parIndex1>fNvar-1 || parIndex1<0)
01031 {
01032 BCLog::Out(BCLog::warning, BCLog::warning,
01033 Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex1));
01034 return -1;
01035 }
01036
01037 if(parIndex2>fNvar-1 || parIndex2<0)
01038 {
01039 BCLog::Out(BCLog::warning, BCLog::warning,
01040 Form("BCIntegrate::GetH2DIndex. Parameter index %d is invalid", parIndex2));
01041 return -1;
01042 }
01043
01044 if(parIndex1==parIndex2)
01045 {
01046 BCLog::Out(BCLog::warning, BCLog::warning,
01047 Form("BCIntegrate::GetH2DIndex. Parameters have equal indeces: %d , %d", parIndex1, parIndex2));
01048 return -1;
01049 }
01050
01051 if(parIndex1>parIndex2)
01052 {
01053 BCLog::Out(BCLog::warning, BCLog::warning,
01054 "BCIntegrate::GetH2DIndex. First parameters must be smaller than second (sorry).");
01055 return -1;
01056 }
01057
01058 int index=0;
01059 for(int i=0;i<fNvar-1;i++)
01060 for(int j=i+1;j<fNvar;j++)
01061 {
01062 if(i==parIndex1 && j==parIndex2)
01063 return index;
01064 index++;
01065 }
01066
01067 BCLog::Out(BCLog::warning, BCLog::warning,
01068 "BCIntegrate::GetH2DIndex. Invalid index combination.");
01069
01070 return -1;
01071 }
01072
01073
01074 TH2D * BCIntegrate::GetH2D(int parIndex1, int parIndex2)
01075 {
01076 if(fHProb2D.size()==0)
01077 {
01078 BCLog::Out(BCLog::warning, BCLog::warning,
01079 "BCModel::GetH2D. MarginalizeAll() has to be run prior to this to fill the distributions.");
01080 return 0;
01081 }
01082
01083 int hindex = this -> GetH2DIndex(parIndex1, parIndex2);
01084 if(hindex==-1)
01085 return 0;
01086
01087 if(hindex>(int)fHProb2D.size()-1)
01088 {
01089 BCLog::Out(BCLog::warning, BCLog::warning,
01090 "BCIntegrate::GetH2D. Got invalid index from GetH2DIndex(). Something went wrong.");
01091 return 0;
01092 }
01093
01094 return fHProb2D[hindex];
01095 }
01096
01097
01098 double BCIntegrate::GetRandomPoint(std::vector <double> &x)
01099 {
01100 this->GetRandomVector(x);
01101
01102 for(int i=0;i<fNvar;i++)
01103 x[i]=fMin[i]+x[i]*(fMax[i]-fMin[i]);
01104
01105 return this->Eval(x);
01106 }
01107
01108
01109 double BCIntegrate::GetRandomPointImportance(std::vector <double> &x)
01110 {
01111 double p = 1.1;
01112 double g = 1.0;
01113
01114 while (p>g)
01115 {
01116 this->GetRandomVector(x);
01117
01118 for(int i=0;i<fNvar;i++)
01119 x[i] = fMin[i] + x[i] * (fMax[i]-fMin[i]);
01120
01121 p = fRandom->Rndm();
01122
01123 g = this -> EvalSampling(x);
01124 }
01125
01126 return this->Eval(x);
01127 }
01128
01129
01130 void BCIntegrate::InitMetro()
01131 {
01132 fNmetro=0;
01133
01134 if (fXmetro0.size() <= 0)
01135 {
01136
01137 for(int i=0;i<fNvar;i++)
01138 fXmetro0.push_back((fMin[i]+fMax[i])/2.0);
01139 }
01140
01141 fMarkovChainValue = this -> LogEval(fXmetro0);
01142
01143
01144 std::vector <double> x;
01145 x.assign(fNvar,0.);
01146
01147 for(int i=0;i<1000;i++)
01148 GetRandomPointMetro(x);
01149
01150 fNmetro = 0;
01151 }
01152
01153
01154 void BCIntegrate::GetRandomPointMetro(std::vector <double> &x)
01155 {
01156
01157 this->GetRandomVectorMetro(fXmetro1);
01158
01159
01160 int in=1;
01161 for(int i=0;i<fNvar;i++)
01162 {
01163 fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]);
01164
01165
01166 if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] )
01167 in=0;
01168 }
01169
01170
01171
01172 double p0 = fMarkovChainValue;
01173 double p1 = 0;
01174 int accept=0;
01175
01176 if(in)
01177 {
01178 p1 = this -> LogEval(fXmetro1);
01179
01180 if(p1>=p0)
01181 accept=1;
01182 else
01183 {
01184 double r=log(fRandom->Rndm());
01185 if(r<p1-p0)
01186 accept=1;
01187 }
01188 }
01189
01190
01191 if(accept)
01192 {
01193
01194 fNacceptedMCMC++;
01195
01196 for(int i=0;i<fNvar;i++)
01197 {
01198 fXmetro0[i]=fXmetro1[i];
01199 x[i]=fXmetro1[i];
01200 fMarkovChainValue = p1;
01201 }
01202 }
01203 else
01204 for(int i=0;i<fNvar;i++)
01205 {
01206 x[i]=fXmetro0[i];
01207 fXmetro1[i] = x[i];
01208 fMarkovChainValue = p0;
01209 }
01210
01211 fNmetro++;
01212 }
01213
01214
01215 void BCIntegrate::GetRandomPointSamplingMetro(std::vector <double> &x)
01216 {
01217
01218 this->GetRandomVectorMetro(fXmetro1);
01219
01220
01221 int in=1;
01222 for(int i=0;i<fNvar;i++)
01223 {
01224 fXmetro1[i] = fXmetro0[i] + fXmetro1[i] * (fMax[i]-fMin[i]);
01225
01226
01227 if( fXmetro1[i]<fMin[i] || fXmetro1[i]>fMax[i] )
01228 in=0;
01229 }
01230
01231
01232 double p0 = this -> LogEvalSampling(fXmetro0);
01233 double p1=0;
01234 if(in)
01235 p1= this -> LogEvalSampling(fXmetro1);
01236
01237
01238 int accept=0;
01239 if(in)
01240 {
01241 if(p1>=p0)
01242 accept=1;
01243 else
01244 {
01245 double r=log(fRandom->Rndm());
01246 if(r<p1-p0)
01247 accept=1;
01248 }
01249 }
01250
01251
01252 if(accept)
01253 for(int i=0;i<fNvar;i++)
01254 {
01255 fXmetro0[i]=fXmetro1[i];
01256 x[i]=fXmetro1[i];
01257 }
01258 else
01259 for(int i=0;i<fNvar;i++)
01260 x[i]=fXmetro0[i];
01261
01262 fNmetro++;
01263 }
01264
01265
01266 void BCIntegrate::GetRandomVector(std::vector <double> &x)
01267 {
01268 double * randx = new double[fNvar];
01269
01270 fRandom -> RndmArray(fNvar, randx);
01271
01272 for(int i=0;i<fNvar;i++)
01273 x[i] = randx[i];
01274
01275 delete[] randx;
01276 randx = 0;
01277 }
01278
01279
01280 void BCIntegrate::GetRandomVectorMetro(std::vector <double> &x)
01281 {
01282 double * randx = new double[fNvar];
01283
01284 fRandom -> RndmArray(fNvar, randx);
01285
01286 for(int i=0;i<fNvar;i++)
01287 x[i] = 2.0 * (randx[i] - 0.5) * fMarkovChainStepSize;
01288
01289 delete[] randx;
01290 randx = 0;
01291 }
01292
01293
01294 TMinuit * BCIntegrate::GetMinuit()
01295 {
01296 if (!fMinuit)
01297 fMinuit = new TMinuit();
01298
01299 return fMinuit;
01300 }
01301
01302
01303
01304 void BCIntegrate::FindModeMinuit(std::vector<double> start, int printlevel)
01305 {
01306 bool have_start = true;
01307
01308
01309 if (int(start.size()) != fNvar)
01310 have_start = false;
01311
01312
01313 global_this = this;
01314
01315
01316 if (fMinuit)
01317 delete fMinuit;
01318 fMinuit = new TMinuit(fNvar);
01319
01320
01321 fMinuit -> SetPrintLevel(printlevel);
01322
01323
01324 fMinuit -> SetFCN(&BCIntegrate::FCNLikelihood);
01325
01326
01327 fMinuit -> SetErrorDef(0.5);
01328
01329
01330 int flag;
01331 for (int i = 0; i < fNvar; i++)
01332 {
01333 double starting_point = (fMin[i]+fMax[i])/2.;
01334 if(have_start)
01335 starting_point = start[i];
01336 fMinuit -> mnparm(i,
01337 fx -> at(i) -> GetName().data(),
01338 starting_point,
01339 (fMax[i]-fMin[i])/100.0,
01340 fMin[i],
01341 fMax[i],
01342 flag);
01343 }
01344
01345
01346
01347
01348
01349 fMinuit -> mnexcm("MIGRAD", fMinuitArglist, 2, flag);
01350
01351
01352
01353
01354
01355 fMinuitErrorFlag = flag;
01356
01357
01358 double probmax = 0;
01359 bool valid = false;
01360
01361 if ( int(fBestFitParameters.size()) == fNvar)
01362 {
01363 valid = true;
01364 for (int i = 0; i < fNvar; ++i)
01365 if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i])
01366 valid= false;
01367
01368 if (valid)
01369 probmax = this -> Eval( fBestFitParameters );
01370 }
01371
01372 std::vector<double> tempvec;
01373 for (int i = 0; i < fNvar; i++)
01374 {
01375 double par;
01376 double parerr;
01377 fMinuit -> GetParameter(i, par, parerr);
01378 tempvec.push_back(par);
01379 }
01380 double probmaxminuit = this -> Eval( tempvec );
01381
01382
01383 if ( (probmaxminuit > probmax && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization)
01384 {
01385 fBestFitParameters.clear();
01386 fBestFitParameterErrors.clear();
01387
01388 for (int i = 0; i < fNvar; i++)
01389 {
01390 double par;
01391 double parerr;
01392 fMinuit -> GetParameter(i, par, parerr);
01393 fBestFitParameters.push_back(par);
01394 fBestFitParameterErrors.push_back(parerr);
01395 }
01396
01397
01398 fOptimizationMethodMode = BCIntegrate::kOptMinuit;
01399 }
01400
01401
01402 delete fMinuit;
01403 fMinuit = 0;
01404
01405 return;
01406 }
01407
01408
01409
01410 void BCIntegrate::FindModeSA(std::vector<double> start)
01411 {
01412
01413
01414
01415 bool have_start = true;
01416 std::vector<double> x, y, best_fit;
01417 double fval_x, fval_y, fval_best_fit;
01418 int t = 1;
01419
01420
01421 if (int(start.size()) != fNvar)
01422 have_start = false;
01423
01424
01425 if ( !have_start )
01426 {
01427 start.clear();
01428 for (int i = 0; i < fNvar; i++)
01429 start.push_back((fMin[i]+fMax[i])/2.);
01430 }
01431
01432
01433 x.clear();
01434 best_fit.clear();
01435 for (int i = 0; i < fNvar; i++)
01436 {
01437 x.push_back(start[i]);
01438 best_fit.push_back(start[i]);
01439 }
01440
01441 fval_x = fval_best_fit = this->LogEval(x);
01442
01443
01444 while ( this->SATemperature(t) > fSATmin )
01445 {
01446
01447 y = this->GetProposalPointSA(x, t);
01448
01449
01450
01451 bool is_in_ranges = true;
01452
01453 for (int i = 0; i < fNvar; i++)
01454 if (y[i] > fMax[i] || y[i] < fMin[i])
01455 is_in_ranges = false;
01456
01457 if ( !is_in_ranges )
01458 ;
01459 else
01460 {
01461
01462 fval_y = this->LogEval(y);
01463
01464
01465
01466 if (fval_y >= fval_x)
01467 {
01468 x.clear();
01469 for (int i = 0; i < fNvar; i++)
01470 x.push_back(y[i]);
01471
01472 fval_x = fval_y;
01473
01474 if (fval_y > fval_best_fit)
01475 {
01476 best_fit.clear();
01477 for (int i = 0; i < fNvar; i++)
01478 best_fit.push_back(y[i]);
01479
01480 fval_best_fit = fval_y;
01481 }
01482 }
01483
01484 else
01485 {
01486 if (fRandom->Rndm() <= exp( (fval_y - fval_x)
01487 / this->SATemperature(t) ))
01488 {
01489 x.clear();
01490 for (int i = 0; i < fNvar; i++)
01491 x.push_back(y[i]);
01492
01493 fval_x = fval_y;
01494 }
01495 }
01496 }
01497
01498
01499 fSANIterations = t;
01500 fSATemperature = this -> SATemperature(t);
01501 fSALogProb = fval_x;
01502 fSAx.clear();
01503 for (int i = 0; i < fNvar; ++i)
01504 fSAx.push_back(x[i]);
01505
01506
01507 if (fFlagWriteSAToFile)
01508 fTreeSA -> Fill();
01509
01510
01511 t++;
01512 }
01513
01514
01515 double probmax = 0;
01516 bool valid = false;
01517
01518 if ( int(fBestFitParameters.size()) == fNvar)
01519 {
01520 valid = true;
01521 for (int i = 0; i < fNvar; ++i)
01522 if (fBestFitParameters.at(i) < fMin[i] || fBestFitParameters.at(i) > fMax[i])
01523 valid= false;
01524
01525 if (valid)
01526 probmax = this -> Eval( fBestFitParameters );
01527 }
01528
01529 double probmaxsa = this -> Eval( best_fit);
01530
01531 if ( (probmaxsa > probmax && !fFlagIgnorePrevOptimization) || !valid || fFlagIgnorePrevOptimization)
01532 {
01533
01534 fBestFitParameters.clear();
01535 fBestFitParameterErrors.clear();
01536
01537 for (int i = 0; i < fNvar; i++)
01538 {
01539 fBestFitParameters.push_back(best_fit[i]);
01540 fBestFitParameterErrors.push_back(0.0);
01541 }
01542
01543
01544 fOptimizationMethodMode = BCIntegrate::kOptSA;
01545 }
01546
01547 return;
01548 }
01549
01550
01551
01552 double BCIntegrate::SATemperature(double t)
01553 {
01554
01555 if (this->fSASchedule == BCIntegrate::kSABoltzmann)
01556 return this->SATemperatureBoltzmann(t);
01557 else if (this->fSASchedule == BCIntegrate::kSACauchy)
01558 return this->SATemperatureCauchy(t);
01559 else
01560 return this->SATemperatureCustom(t);
01561 }
01562
01563
01564
01565 double BCIntegrate::SATemperatureBoltzmann(double t)
01566 {
01567 return fSAT0 / log((double)(t + 1));
01568 }
01569
01570
01571
01572 double BCIntegrate::SATemperatureCauchy(double t)
01573 {
01574 return fSAT0 / (double)t;
01575 }
01576
01577
01578
01579 double BCIntegrate::SATemperatureCustom(double t)
01580 {
01581 BCLog::OutError("BCIntegrate::SATemperatureCustom : No custom temperature schedule defined");
01582 return 0.;
01583 }
01584
01585
01586
01587 std::vector<double> BCIntegrate::GetProposalPointSA(std::vector<double> x, int t)
01588 {
01589
01590 if (this->fSASchedule == BCIntegrate::kSABoltzmann)
01591 return this->GetProposalPointSABoltzmann(x, t);
01592 else if (this->fSASchedule == BCIntegrate::kSACauchy)
01593 return this->GetProposalPointSACauchy(x, t);
01594 else
01595 return this->GetProposalPointSACustom(x, t);
01596 }
01597
01598
01599
01600 std::vector<double> BCIntegrate::GetProposalPointSABoltzmann(std::vector<double> x, int t)
01601 {
01602 std::vector<double> y;
01603 y.clear();
01604 double new_val, norm;
01605
01606 for (int i = 0; i < fNvar; i++)
01607 {
01608 norm = (fMax[i] - fMin[i]) * this->SATemperature(t) / 2.;
01609 new_val = x[i] + norm * fRandom->Gaus();
01610 y.push_back(new_val);
01611 }
01612
01613 return y;
01614 }
01615
01616
01617
01618 std::vector<double> BCIntegrate::GetProposalPointSACauchy(std::vector<double> x, int t)
01619 {
01620 std::vector<double> y;
01621 y.clear();
01622
01623 if (fNvar == 1)
01624 {
01625 double cauchy, new_val, norm;
01626
01627 norm = (fMax[0] - fMin[0]) * this->SATemperature(t) / 2.;
01628 cauchy = tan(3.14159 * (fRandom->Rndm() - 0.5));
01629 new_val = x[0] + norm * cauchy;
01630 y.push_back(new_val);
01631 }
01632 else
01633 {
01634
01635
01636
01637
01638 y = this->SAHelperGetRandomPointOnHypersphere();
01639
01640
01641
01642 double radial = this->SATemperature(t) * this->SAHelperGetRadialCauchy();
01643
01644
01645
01646 for (int i = 0; i < fNvar; i++)
01647 y[i] = (fMax[i] - fMin[i]) * y[i] * radial / 2. + x[i];
01648 }
01649
01650 return y;
01651 }
01652
01653
01654
01655 std::vector<double> BCIntegrate::GetProposalPointSACustom(std::vector<double> x, int t)
01656 {
01657 BCLog::OutError("BCIntegrate::GetProposalPointSACustom : No custom proposal function defined");
01658 return std::vector<double>(fNvar);
01659 }
01660
01661
01662
01663 std::vector<double> BCIntegrate::SAHelperGetRandomPointOnHypersphere()
01664 {
01665 std::vector<double> rand_point(fNvar);
01666
01667
01668
01669
01670
01671
01672
01673 if (fNvar == 2)
01674 {
01675 double x1, x2, s;
01676 do {
01677 x1 = fRandom->Rndm() * 2. - 1.;
01678 x2 = fRandom->Rndm() * 2. - 1.;
01679 s = x1*x1 + x2*x2;
01680 } while (s >= 1);
01681
01682 rand_point[0] = (x1*x1 - x2*x2) / s;
01683 rand_point[1] = (2.*x1*x2) / s;
01684 }
01685 else if (fNvar == 3)
01686 {
01687 fRandom->Sphere(rand_point[0], rand_point[1], rand_point[2], 1.0);
01688 }
01689 else
01690 {
01691 double s = 0.,
01692 gauss_num;
01693
01694 for (int i = 0; i < fNvar; i++)
01695 {
01696 gauss_num = fRandom->Gaus();
01697 rand_point[i] = gauss_num;
01698 s += gauss_num * gauss_num;
01699 }
01700 s = sqrt(s);
01701
01702 for (int i = 0; i < fNvar; i++)
01703 rand_point[i] = rand_point[i] / s;
01704 }
01705
01706 return rand_point;
01707 }
01708
01709
01710
01711 double BCIntegrate::SAHelperGetRadialCauchy()
01712 {
01713
01714
01715
01716
01717 double theta;
01718
01719
01720 static double map_u[10001];
01721 static double map_theta[10001];
01722 static bool initialized = false;
01723 static int map_dimension = 0;
01724
01725
01726 if (!initialized || map_dimension != fNvar)
01727 {
01728 double init_theta;
01729 double beta = this->SAHelperSinusToNIntegral(fNvar - 1, 1.57079632679);
01730
01731 for (int i = 0; i <= 10000; i++)
01732 {
01733 init_theta = 3.14159265 * (double)i / 5000.;
01734 map_theta[i] = init_theta;
01735
01736 map_u[i] = this->SAHelperSinusToNIntegral(fNvar - 1, init_theta) / beta;
01737 }
01738
01739 map_dimension = fNvar;
01740 initialized = true;
01741 }
01742
01743
01744 double u = fRandom->Rndm();
01745
01746
01747
01748 int lo = 0;
01749 int up = 10000;
01750 int mid;
01751
01752 while (up != lo)
01753 {
01754 mid = ((up - lo + 1) / 2) + lo;
01755
01756 if (u >= map_u[mid])
01757 lo = mid;
01758 else
01759 up = mid - 1;
01760 }
01761 up++;
01762
01763
01764 theta = map_theta[lo] + (u - map_u[lo]) / (map_u[up] - map_u[lo])
01765 * (map_theta[up] - map_theta[lo]);
01766
01767 return tan(theta);
01768 }
01769
01770
01771
01772 double BCIntegrate::SAHelperSinusToNIntegral(int dim, double theta)
01773 {
01774 if (dim < 1)
01775 return theta;
01776 else if (dim == 1)
01777 return (1. - cos(theta));
01778 else if (dim == 2)
01779 return 0.5 * (theta - sin(theta) * cos(theta));
01780 else if (dim == 3)
01781 return (2. - sin(theta) * sin(theta) * cos(theta) - 2. * cos(theta)) / 3.;
01782 else
01783 return - pow(sin(theta), (double)(dim - 1)) * cos(theta) / (double)dim
01784 + (double)(dim - 1) / (double)dim
01785 * this->SAHelperSinusToNIntegral(dim - 2, theta);
01786 }
01787
01788
01789
01790 void BCIntegrate::SetMode(std::vector <double> mode)
01791 {
01792 if((int)mode.size() == fNvar)
01793 {
01794 fBestFitParameters.clear();
01795 for (int i = 0; i < fNvar; i++)
01796 fBestFitParameters.push_back(mode[i]);
01797 }
01798 }
01799
01800
01801
01802 void BCIntegrate::FCNLikelihood(int &npar, double * grad, double &fval, double * par, int flag)
01803 {
01804
01805 std::vector <double> parameters;
01806
01807 int n = global_this -> GetNvar();
01808
01809 for (int i = 0; i < n; i++)
01810 parameters.push_back(par[i]);
01811
01812 fval = - global_this -> LogEval(parameters);
01813
01814
01815 parameters.clear();
01816 }
01817
01818
01819
01820
01821 void BCIntegrate::FindModeMCMC()
01822 {
01823
01824
01825
01826
01827
01828
01829
01830
01831 double probmax = 0;
01832
01833 if ( int(fBestFitParameters.size()) == fNvar)
01834 {
01835 probmax = this -> Eval( fBestFitParameters );
01836
01837 }
01838
01839
01840 for (int i = 0; i < fMCMCNChains; ++i)
01841 {
01842 double prob = exp( (this -> MCMCGetMaximumLogProb()).at(i));
01843
01844
01845 if ( (prob >= probmax && !fFlagIgnorePrevOptimization) || fFlagIgnorePrevOptimization)
01846 {
01847 probmax = prob;
01848
01849 fBestFitParameters.clear();
01850 fBestFitParameterErrors.clear();
01851 fBestFitParameters = this -> MCMCGetMaximumPoint(i);
01852
01853 for (int j = 0; j < fNvar; ++j)
01854 fBestFitParameterErrors.push_back(0.0);
01855
01856 fOptimizationMethodMode = BCIntegrate::kOptMetropolis;
01857 }
01858 }
01859
01860 }
01861
01862
01863 void BCIntegrate::CubaIntegrand(const int *ndim, const double xx[],
01864 const int *ncomp, double ff[])
01865 {
01866 #ifdef HAVE_CUBA_H
01867
01868 double jacobian = 1.0;
01869
01870 std::vector<double> scaled_parameters;
01871
01872 for (int i = 0; i < *ndim; i++)
01873 {
01874 double range = global_this -> fx -> at(i) -> GetUpperLimit() - global_this -> fx -> at(i) -> GetLowerLimit();
01875
01876
01877 jacobian *= range;
01878
01879
01880 scaled_parameters.push_back(global_this -> fx -> at(i) -> GetLowerLimit() + xx[i] * range);
01881 }
01882
01883
01884 ff[0] = global_this -> Eval(scaled_parameters);
01885
01886
01887 ff[0] *= jacobian;
01888
01889
01890 ff[0] *= 1e99;
01891
01892
01893 scaled_parameters.clear();
01894 #else
01895 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
01896 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT.");
01897 #endif
01898 }
01899
01900
01901 double BCIntegrate::CubaIntegrate()
01902 {
01903 #ifdef HAVE_CUBA_H
01904 double EPSREL = 1e-3;
01905 double EPSABS = 1e-12;
01906 double VERBOSE = 0;
01907 double MINEVAL = 0;
01908 double MAXEVAL = 2000000;
01909 double NSTART = 25000;
01910 double NINCREASE = 25000;
01911
01912 std::vector<double> parameters_double;
01913 std::vector<double> parameters_int;
01914
01915 parameters_double.push_back(EPSREL);
01916 parameters_double.push_back(EPSABS);
01917
01918 parameters_int.push_back(VERBOSE);
01919 parameters_int.push_back(MINEVAL);
01920 parameters_int.push_back(MAXEVAL);
01921 parameters_int.push_back(NSTART);
01922 parameters_int.push_back(NINCREASE);
01923
01924
01925 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Running Cuba/Vegas integation over %i dimensions.", fNvar));
01926 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Maximum number of iterations: %i", (int)MAXEVAL));
01927 BCLog::Out(BCLog::detail, BCLog::detail, Form(" --> Aimed relative precision: %e", EPSREL));
01928
01929 return this -> CubaIntegrate(0, parameters_double, parameters_int);
01930 #else
01931 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
01932 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT.");
01933 return 0.;
01934 #endif
01935 }
01936
01937
01938 double BCIntegrate::CubaIntegrate(int method, std::vector<double> parameters_double, std::vector<double> parameters_int)
01939 {
01940 #ifdef HAVE_CUBA_H
01941 const int NDIM = int(fx ->size());
01942 const int NCOMP = 1;
01943
01944 const double EPSREL = parameters_double[0];
01945 const double EPSABS = parameters_double[1];
01946 const int VERBOSE = int(parameters_int[0]);
01947 const int MINEVAL = int(parameters_int[1]);
01948 const int MAXEVAL = int(parameters_int[2]);
01949
01950 int neval;
01951 int fail;
01952 int nregions;
01953 double integral[NCOMP];
01954 double error[NCOMP];
01955 double prob[NCOMP];
01956
01957 global_this = this;
01958
01959 integrand_t an_integrand = &BCIntegrate::CubaIntegrand;
01960
01961 if (method == 0)
01962 {
01963
01964 const int NSTART = int(parameters_int[3]);
01965 const int NINCREASE = int(parameters_int[4]);
01966
01967
01968 Vegas(NDIM, NCOMP, an_integrand,
01969 EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
01970 NSTART, NINCREASE,
01971 &neval, &fail, integral, error, prob);
01972 }
01973 else if (method == 1)
01974 {
01975
01976 const int LAST = int(parameters_int[3]);
01977 const int NNEW = int(parameters_int[4]);
01978 const int FLATNESS = int(parameters_int[5]);
01979
01980
01981 Suave(NDIM, NCOMP, an_integrand,
01982 EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL,
01983 NNEW, FLATNESS,
01984 &nregions, &neval, &fail, integral, error, prob);
01985 }
01986 else if (method == 2)
01987 {
01988
01989 const int KEY1 = int(parameters_int[3]);
01990 const int KEY2 = int(parameters_int[4]);
01991 const int KEY3 = int(parameters_int[5]);
01992 const int MAXPASS = int(parameters_int[6]);
01993 const int BORDER = int(parameters_int[7]);
01994 const int MAXCHISQ = int(parameters_int[8]);
01995 const int MINDEVIATION = int(parameters_int[9]);
01996 const int NGIVEN = int(parameters_int[10]);
01997 const int LDXGIVEN = int(parameters_int[11]);
01998 const int NEXTRA = int(parameters_int[12]);
01999
02000
02001 Divonne(NDIM, NCOMP, an_integrand,
02002 EPSREL, EPSABS, VERBOSE, MINEVAL, MAXEVAL,
02003 KEY1, KEY2, KEY3, MAXPASS, BORDER, MAXCHISQ, MINDEVIATION,
02004 NGIVEN, LDXGIVEN, NULL, NEXTRA, NULL,
02005 &nregions, &neval, &fail, integral, error, prob);
02006 }
02007 else if (method == 3)
02008 {
02009
02010 const int LAST = int(parameters_int[3]);
02011 const int KEY = int(parameters_int[4]);
02012
02013
02014 Cuhre(NDIM, NCOMP, an_integrand,
02015 EPSREL, EPSABS, VERBOSE | LAST, MINEVAL, MAXEVAL, KEY,
02016 &nregions, &neval, &fail, integral, error, prob);
02017 }
02018 else
02019 {
02020 std::cout << " Integration method not available. " << std::endl;
02021 integral[0] = -1e99;
02022 }
02023
02024 if (fail != 0)
02025 {
02026 std::cout << " Warning, integral did not converge with the given set of parameters. "<< std::endl;
02027 std::cout << " neval = " << neval << std::endl;
02028 std::cout << " fail = " << fail << std::endl;
02029 std::cout << " integral = " << integral[0] << std::endl;
02030 std::cout << " error = " << error[0] << std::endl;
02031 std::cout << " prob = " << prob[0] << std::endl;
02032 }
02033
02034 return integral[0] / 1e99;
02035 #else
02036 BCLog::Out(BCLog::error,BCLog::error,"!!! This version of BAT is compiled without Cuba.");
02037 BCLog::Out(BCLog::error,BCLog::error," Use other integration methods or install Cuba and recompile BAT.");
02038 return 0.;
02039 #endif
02040 }
02041
02042
02043 void BCIntegrate::MCMCIterationInterface()
02044 {
02045
02046
02047
02048
02049 this -> MCMCFillErrorBand();
02050
02051
02052 this -> MCMCUserIterationInterface();
02053 }
02054
02055
02056 void BCIntegrate::MCMCFillErrorBand()
02057 {
02058 if (!fFillErrorBand)
02059 return;
02060
02061
02062 if (fFitFunctionIndexX < 0)
02063 return;
02064
02065
02066 if (fErrorBandContinuous)
02067 {
02068 double x = 0;
02069 for (int ix = 0; ix < fErrorBandNbinsX; ix++)
02070 {
02071
02072 x = fErrorBandXY -> GetXaxis() -> GetBinCenter(ix + 1);
02073
02074
02075 std::vector <double> xvec;
02076 xvec.push_back(x);
02077
02078
02079 for (int ichain = 0; ichain < this -> MCMCGetNChains(); ++ichain)
02080 {
02081
02082 double y = this -> FitFunction(xvec, this -> MCMCGetx(ichain));
02083
02084
02085 fErrorBandXY -> Fill(x, y);
02086 }
02087
02088 xvec.clear();
02089 }
02090 }
02091
02092 else
02093 {
02094 int ndatapoints = int(fErrorBandX.size());
02095 double x = 0;
02096
02097 for (int ix = 0; ix < ndatapoints; ++ix)
02098 {
02099
02100 x = fErrorBandX.at(ix);
02101
02102
02103 std::vector <double> xvec;
02104 xvec.push_back(x);
02105
02106
02107 for (int ichain = 0; ichain < this -> MCMCGetNChains(); ++ichain)
02108 {
02109
02110 double y = this -> FitFunction(xvec, this -> MCMCGetx(ichain));
02111
02112
02113 fErrorBandXY -> Fill(x, y);
02114 }
02115
02116 xvec.clear();
02117 }
02118 }
02119 }
02120
02121
02122 void BCIntegrate::SAInitialize()
02123 {
02124 fSAx.clear();
02125 fSAx.assign(fNvar, 0.0);
02126 }
02127
02128
02129