00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "BAT/BCModelManager.h"
00011
00012 #include "BAT/BCModel.h"
00013 #include "BAT/BCDataSet.h"
00014 #include "BAT/BCDataPoint.h"
00015 #include "BAT/BCLog.h"
00016 #include "BAT/BCErrorCodes.h"
00017
00018 #include <fstream>
00019 #include <iostream>
00020
00021 #include <TString.h>
00022
00023
00024
00025 BCModelManager::BCModelManager()
00026 {
00027 fModelContainer = new BCModelContainer();
00028 fDataSet = 0;
00029 }
00030
00031
00032
00033 BCModelManager::~BCModelManager()
00034 {
00035 delete fModelContainer;
00036
00037 if (fDataSet)
00038 delete fDataSet;
00039 }
00040
00041
00042
00043 BCModelManager::BCModelManager(const BCModelManager & modelmanager)
00044 {
00045 modelmanager.Copy(*this);
00046 }
00047
00048
00049
00050 BCModelManager & BCModelManager::operator = (const BCModelManager & modelmanager)
00051 {
00052 if (this != &modelmanager)
00053 modelmanager.Copy(* this);
00054
00055 return * this;
00056 }
00057
00058
00059
00060 void BCModelManager::SetDataSet(BCDataSet * dataset)
00061 {
00062
00063 fDataSet = dataset;
00064
00065
00066 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00067 this -> GetModel(i) -> SetDataSet(fDataSet);
00068 }
00069
00070
00071
00072 void BCModelManager::SetSingleDataPoint(BCDataPoint * datapoint)
00073 {
00074
00075 BCDataSet * dataset = new BCDataSet();
00076
00077
00078 dataset -> AddDataPoint(datapoint);
00079
00080
00081 this -> SetDataSet(dataset);
00082 }
00083
00084
00085
00086 void BCModelManager::SetSingleDataPoint(BCDataSet * dataset, unsigned int index)
00087 {
00088 if (index < 0 || index > dataset -> GetNDataPoints())
00089 return;
00090
00091 this -> SetSingleDataPoint(dataset -> GetDataPoint(index));
00092 }
00093
00094
00095
00096 void BCModelManager::AddModel(BCModel * model, double probability)
00097 {
00098
00099 unsigned int index = fModelContainer -> size();
00100
00101
00102 model -> SetIndex(index);
00103
00104
00105 model -> SetModelAPrioriProbability(probability);
00106
00107
00108 model -> SetDataSet(fDataSet);
00109
00110
00111 fModelContainer -> push_back(model);
00112 }
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 void BCModelManager::SetIntegrationMethod(BCIntegrate::BCIntegrationMethod method)
00129 {
00130
00131
00132 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00133 this -> GetModel(i) -> SetIntegrationMethod(method);
00134 }
00135
00136
00137
00138 void BCModelManager::SetMarginalizationMethod(BCIntegrate::BCMarginalizationMethod method)
00139 {
00140
00141 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00142 this -> GetModel(i) -> SetMarginalizationMethod(method);
00143 };
00144
00145
00146
00147 void BCModelManager::SetOptimizationMethod(BCIntegrate::BCOptimizationMethod method)
00148 {
00149
00150 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00151 this -> GetModel(i) -> SetOptimizationMethod(method);
00152 }
00153
00154
00155
00156 void BCModelManager::SetNiterationsPerDimension(unsigned int niterations)
00157 {
00158
00159 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00160 this -> GetModel(i) -> SetNiterationsPerDimension(niterations);
00161 }
00162
00163
00164
00165 void BCModelManager::SetNSamplesPer2DBin(unsigned int n)
00166 {
00167
00168 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00169 this -> GetModel(i) -> SetNSamplesPer2DBin(n);
00170 }
00171
00172
00173
00174 void BCModelManager::SetRelativePrecision(double relprecision)
00175 {
00176
00177 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00178 this -> GetModel(i) -> SetRelativePrecision(relprecision);
00179 }
00180
00181
00182
00183 void BCModelManager::SetNbins(unsigned int n)
00184 {
00185
00186 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00187 this -> GetModel(i) -> BCIntegrate::SetNbins(n);
00188 }
00189
00190
00191
00192 void BCModelManager::SetFillErrorBand(bool flag)
00193 {
00194 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00195 this -> GetModel(i) -> SetFillErrorBand(flag);
00196 }
00197
00198
00199
00200 void BCModelManager::SetFitFunctionIndexX(int index)
00201 {
00202
00203 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00204 this -> GetModel(i) -> SetFitFunctionIndexX(index);
00205 }
00206
00207
00208
00209 void BCModelManager::SetFitFunctionIndexY(int index)
00210 {
00211
00212 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00213 this -> GetModel(i) -> SetFitFunctionIndexY(index);
00214 }
00215
00216
00217
00218 void BCModelManager::SetFitFunctionIndices(int indexx, int indexy)
00219 {
00220
00221 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00222 this -> GetModel(i) -> SetFitFunctionIndices(indexx, indexy);
00223 }
00224
00225
00226
00227 void BCModelManager::SetDataPointLowerBoundaries(BCDataPoint * datasetlowerboundaries)
00228 {
00229
00230 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00231 this -> GetModel(i) -> SetDataPointLowerBoundaries(datasetlowerboundaries);
00232 }
00233
00234
00235
00236 void BCModelManager::SetDataPointUpperBoundaries(BCDataPoint * datasetupperboundaries)
00237 {
00238
00239 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00240 this -> GetModel(i) -> SetDataPointUpperBoundaries(datasetupperboundaries);
00241 }
00242
00243
00244
00245 void BCModelManager::SetDataPointLowerBoundary(int index, double lowerboundary)
00246 {
00247
00248 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00249 this -> GetModel(i) -> SetDataPointLowerBoundary(index, lowerboundary);
00250 }
00251
00252
00253
00254 void BCModelManager::SetDataPointUpperBoundary(int index, double upperboundary)
00255 {
00256
00257 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00258 this -> GetModel(i) -> SetDataPointUpperBoundary(index, upperboundary);
00259 }
00260
00261
00262
00263 void BCModelManager::SetDataBoundaries(int index, double lowerboundary, double upperboundary)
00264 {
00265
00266 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00267 this -> GetModel(i) -> SetDataBoundaries(index, lowerboundary, upperboundary);
00268 }
00269
00270
00271
00272 void BCModelManager::FixDataAxis(int index, bool fixed)
00273 {
00274
00275 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00276 this -> GetModel(i) -> FixDataAxis(index, fixed);
00277 }
00278
00279
00280
00281 void BCModelManager::SetNChains(unsigned int n)
00282 {
00283
00284 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00285 this -> GetModel(i) -> MCMCSetNChains(n);
00286 }
00287
00288
00289
00290 int BCModelManager::ReadDataFromFileTree(const char * filename, const char * treename, const char * branchnames)
00291 {
00292 if (fModelContainer -> size() < 0)
00293 {
00294 BCLog::Out(BCLog::warning, BCLog::warning, "BCModelManager::ReadDataFromFileTree : No model defined.");
00295 return ERROR_NOMODELS;
00296 }
00297
00298
00299 if (!fDataSet)
00300 fDataSet = new BCDataSet();
00301 else
00302 fDataSet -> Reset();
00303
00304
00305 int read_file = fDataSet -> ReadDataFromFileTree(filename, treename, branchnames);
00306
00307 if (read_file >=0)
00308 {
00309 this -> SetDataSet(fDataSet);
00310
00311 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00312 fModelContainer -> at(i) -> SetDataSet(fDataSet);
00313 }
00314 else if (read_file == ERROR_FILENOTFOUND)
00315 {
00316 delete fDataSet;
00317 return ERROR_FILENOTFOUND;
00318 }
00319
00320 return 0;
00321 }
00322
00323
00324
00325 int BCModelManager::ReadDataFromFileTxt(const char * filename, int nbranches)
00326 {
00327 if (fModelContainer -> size() < 0)
00328 {
00329 BCLog::Out(BCLog::warning, BCLog::warning, "BCModelManager::ReadDataFromFileTree. No model defined.");
00330 return ERROR_NOMODELS;
00331 }
00332
00333
00334 if (!fDataSet)
00335 fDataSet = new BCDataSet();
00336 else
00337 fDataSet -> Reset();
00338
00339
00340 int read_file = fDataSet -> ReadDataFromFileTxt(filename, nbranches);
00341
00342 if (read_file >=0)
00343 {
00344 this -> SetDataSet(fDataSet);
00345
00346 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00347 fModelContainer -> at(i) -> SetDataSet(fDataSet);
00348 }
00349 else
00350 {
00351 delete fDataSet;
00352 return ERROR_FILENOTFOUND;
00353 }
00354
00355 return 0;
00356 }
00357
00358
00359
00360 void BCModelManager::Normalize()
00361 {
00362
00363 double normalization = 0.0;
00364
00365
00366 BCLog::OutSummary("Running normalization of all models.");
00367
00368 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00369 {
00370 fModelContainer -> at(i) -> Normalize();
00371
00372
00373 normalization += (fModelContainer -> at(i) -> GetNormalization() *
00374 fModelContainer -> at(i) -> GetModelAPrioriProbability());
00375 }
00376
00377
00378 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00379 fModelContainer -> at(i) -> SetModelAPosterioriProbability(
00380 (fModelContainer -> at(i) -> GetNormalization() *
00381 fModelContainer -> at(i) -> GetModelAPrioriProbability()) /
00382 normalization);
00383 }
00384
00385
00386
00387 double BCModelManager::BayesFactor(const unsigned int imodel1, const unsigned int imodel2)
00388 {
00389
00390
00391
00392
00393
00394 const double norm1 = fModelContainer -> at(imodel1) -> GetNormalization();
00395 const double norm2 = fModelContainer -> at(imodel2) -> GetNormalization();
00396
00397
00398 if(norm1<0.)
00399 {
00400 BCLog::OutError(
00401 Form("Model %s (index %d) not normalized. Cannot calculate Bayes factor.",
00402 fModelContainer->at(imodel1)->GetName().data(),imodel1));
00403 return -1.;
00404 }
00405
00406
00407 if(norm2<0.)
00408 {
00409 BCLog::OutError(
00410 Form("Model %s (index %d) not normalized. Cannot calculate Bayes factor.",
00411 fModelContainer->at(imodel2)->GetName().data(),imodel2));
00412 return -1.;
00413 }
00414
00415
00416 if(norm2==0. && norm1!=0.)
00417 {
00418 BCLog::OutError(
00419 Form("Model %s (index %d) has ZERO probability. Bayes factor is infinite.",
00420 fModelContainer->at(imodel2)->GetName().data(),imodel2));
00421 return -1.;
00422 }
00423
00424
00425 if(norm2==0. && norm1==0.)
00426 {
00427 BCLog::OutWarning(
00428 Form("Models %s and %s have ZERO probability. Bayes factor is unknown. Returning 1.",
00429 fModelContainer->at(imodel2)->GetName().data(),fModelContainer->at(imodel1)->GetName().data()));
00430 return 1.;
00431 }
00432
00433
00434 return norm1/norm2;
00435 }
00436
00437
00438
00439 void BCModelManager::FindMode()
00440 {
00441
00442 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00443 this -> GetModel(i) -> FindMode();
00444 }
00445
00446
00447
00448 void BCModelManager::MarginalizeAll()
00449 {
00450
00451 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00452 this -> GetModel(i) -> MarginalizeAll();
00453 }
00454
00455
00456
00457 void BCModelManager::WriteMarkovChain(bool flag)
00458 {
00459
00460 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00461 this -> GetModel(i) -> WriteMarkovChain(flag);
00462 }
00463
00464
00465
00466 void BCModelManager::CalculatePValue(bool flag_histogram)
00467 {
00468
00469 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00470 this -> GetModel(i) -> CalculatePValue(this -> GetModel(i) -> GetBestFitParameters(), flag_histogram);
00471 }
00472
00473
00474
00475 void BCModelManager::PrintSummary(const char * file)
00476 {
00477 ofstream out;
00478 std::streambuf * old_buffer = 0;
00479
00480 if(file)
00481 {
00482 out.open(file);
00483 if (!out.is_open())
00484 {
00485 std::cerr<<"Couldn't open file "<<file<<std::endl;
00486 return;
00487 }
00488 old_buffer = std::cout.rdbuf(out.rdbuf());
00489 }
00490
00491
00492 int nmodels = fModelContainer -> size();
00493 std::cout
00494 <<std::endl
00495 <<"======================================"<<std::endl
00496 <<" Summary"<<std::endl
00497 <<"======================================"<<std::endl
00498 <<std::endl
00499 <<" Number of models : "<<nmodels<<std::endl
00500 <<std::endl
00501 <<" - Models:"<<std::endl;
00502
00503 for (int i = 0; i < nmodels; i++)
00504 fModelContainer -> at(i) -> PrintSummary();
00505
00506
00507 std::cout
00508 <<" - Data:"<<std::endl
00509 <<std::endl
00510 <<" Number of entries: "<<fDataSet -> GetNDataPoints()<<std::endl
00511 <<std::endl;
00512
00513 std::cout
00514 <<"======================================"<<std::endl
00515 <<" Model comparison"<<std::endl
00516 <<std::endl;
00517
00518
00519 std::cout
00520 <<" - A priori probabilities:"<<std::endl
00521 <<std::endl;
00522
00523 for (int i=0; i<nmodels; i++)
00524 std::cout
00525 <<" p("<< fModelContainer -> at(i) -> GetName()
00526 <<") = "<< fModelContainer -> at(i) -> GetModelAPrioriProbability()
00527 <<std::endl;
00528 std::cout<<std::endl;
00529
00530 std::cout
00531 <<" - A posteriori probabilities:"<<std::endl
00532 <<std::endl;
00533
00534 for (int i = 0; i < nmodels; i++)
00535 std::cout
00536 <<" p("<< fModelContainer -> at(i) -> GetName()
00537 <<" | data) = "<< fModelContainer -> at(i) -> GetModelAPosterioriProbability()
00538 <<std::endl;
00539 std::cout<<std::endl;
00540
00541 std::cout
00542 <<"======================================"<<std::endl
00543 <<std::endl;
00544
00545 if (file)
00546 std::cout.rdbuf(old_buffer);
00547
00548 }
00549
00550
00551
00552 void BCModelManager::PrintModelComparisonSummary(const char * file)
00553 {
00554 ofstream out;
00555 std::streambuf * old_buffer = 0;
00556
00557 if(file)
00558 {
00559 out.open(file);
00560 if (!out.is_open())
00561 {
00562 std::cerr<<"Couldn't open file "<<file<<std::endl;
00563 return;
00564 }
00565 old_buffer = std::cout.rdbuf(out.rdbuf());
00566 }
00567
00568
00569 int nmodels = fModelContainer -> size();
00570 std::cout
00571 <<std::endl
00572 <<"==========================================="<<std::endl
00573 <<" Model Comparison Summary"<<std::endl
00574 <<"==========================================="<<std::endl
00575 <<std::endl
00576 <<" Number of models : "<<nmodels<<std::endl
00577 <<std::endl;
00578
00579
00580 std::cout
00581 <<" - A priori probabilities:"<<std::endl
00582 <<std::endl;
00583
00584 for (int i=0; i<nmodels; i++)
00585 std::cout
00586 <<" p("<< fModelContainer -> at(i) -> GetName()
00587 <<") = "<< fModelContainer -> at(i) -> GetModelAPrioriProbability()
00588 <<std::endl;
00589 std::cout<<std::endl;
00590
00591 std::cout
00592 <<" - A posteriori probabilities:"<<std::endl
00593 <<std::endl;
00594
00595 for (int i = 0; i < nmodels; i++)
00596 std::cout
00597 <<" p("<< fModelContainer -> at(i) -> GetName()
00598 <<" | data) = "<< fModelContainer -> at(i) -> GetModelAPosterioriProbability()
00599 <<std::endl;
00600 std::cout<<std::endl;
00601
00602
00603 std::cout
00604 <<" - Bayes factors:"<<std::endl
00605 <<std::endl;
00606 for (int i = 0; i < nmodels-1; i++)
00607 for (int j = i+1; j < nmodels; j++)
00608 std::cout
00609 <<" K = p(data | "<<fModelContainer -> at(i) -> GetName()<<") / "
00610 <<"p(data | "<<fModelContainer -> at(j) -> GetName()<<") = "
00611 <<BayesFactor(i,j)<<std::endl;
00612 std::cout<<std::endl;
00613
00614
00615 std::cout
00616 <<" - p-values:"<<std::endl
00617 <<std::endl;
00618
00619 for (int i = 0; i < nmodels; i++)
00620 {
00621 double p = fModelContainer -> at(i) -> GetPValue();
00622 std::cout
00623 <<" "<< fModelContainer -> at(i) -> GetName();
00624 if(p>=0.)
00625 std::cout<<": p-value = "<< p;
00626 else
00627 std::cout<<": p-value not calculated";
00628 std::cout<<std::endl;
00629 }
00630 std::cout<<std::endl;
00631
00632 std::cout
00633 <<"==========================================="<<std::endl
00634 <<std::endl;
00635
00636 if (file)
00637 std::cout.rdbuf(old_buffer);
00638
00639 }
00640
00641
00642
00643 void BCModelManager::PrintResults()
00644 {
00645
00646 for (unsigned int i = 0; i < this -> GetNModels(); i++)
00647 this -> GetModel(i) -> PrintResults(Form("%s.txt", this -> GetModel(i) -> GetName().data()));
00648 }
00649
00650
00651
00652 void BCModelManager::Copy(BCModelManager & modelmanager) const
00653 {
00654
00655 modelmanager.fModelContainer = this -> fModelContainer;
00656 modelmanager.fDataSet = this -> fDataSet;
00657 }
00658
00659