BAT  0.9.4
The Bayesian analysis toolkit
 All Classes Namespaces Functions Variables Enumerations
BCSummaryTool.cxx
1 /*
2  * Copyright (C) 2007-2014, the BAT core developer team
3  * All rights reserved.
4  *
5  * For the licensing terms see doc/COPYING.
6  * For documentation see http://mpp.mpg.de/bat
7  */
8 
9 // ---------------------------------------------------------
10 
11 #include "BCSummaryTool.h"
12 #include <string>
13 
14 #include "BCH1D.h"
15 #include "BCH2D.h"
16 #include "BCLog.h"
17 #include "BCMath.h"
18 #include "BCModel.h"
19 #include "BCParameter.h"
20 #include "BCSummaryPriorModel.h"
21 
22 #include <TArrow.h>
23 #include <TCanvas.h>
24 #include <TF1.h>
25 #include <TGaxis.h>
26 #include <TGraphAsymmErrors.h>
27 #include <TGraphErrors.h>
28 #include <TH2D.h>
29 #include <TLatex.h>
30 #include <TLegend.h>
31 #include <TLine.h>
32 #include <TMarker.h>
33 #include <TPostScript.h>
34 #include <TStyle.h>
35 
36 #include <fstream>
37 #include <iostream>
38 
39 unsigned int BCSummaryTool::fHCounter=0;
40 
41 // ---------------------------------------------------------
42 BCSummaryTool::BCSummaryTool()
43  : fModel(0)
44  , fPriorModel(0)
45  , fFlagInfoMarg(false)
46  , fFlagInfoOpt(false)
47 {
48  // define sum of probabilities for quantiles
49  fSumProb.push_back(0.05);
50  fSumProb.push_back(0.10);
51  fSumProb.push_back(0.1587);
52  fSumProb.push_back(0.50);
53  fSumProb.push_back(0.8413);
54  fSumProb.push_back(0.90);
55  fSumProb.push_back(0.95);
56 
57  // set text style
58  gStyle->SetPaintTextFormat(".2g");
59 }
60 
61 // ---------------------------------------------------------
62 BCSummaryTool::BCSummaryTool(BCModel * model)
63  : fModel(model)
64  , fPriorModel(0)
65  , fFlagInfoMarg(false)
66  , fFlagInfoOpt(false)
67 {
68  // define sum of probabilities for quantiles
69  fSumProb.push_back(0.05);
70  fSumProb.push_back(0.10);
71  fSumProb.push_back(0.1587);
72  fSumProb.push_back(0.50);
73  fSumProb.push_back(0.8413);
74  fSumProb.push_back(0.90);
75  fSumProb.push_back(0.95);
76 
77  // set text style
78  gStyle->SetPaintTextFormat(".2g");
79 }
80 
81 // ---------------------------------------------------------
82 BCSummaryTool::~BCSummaryTool()
83 {
84  delete fPriorModel;
85 }
86 
87 // ---------------------------------------------------------
88 int BCSummaryTool::CopySummaryData()
89 {
90  // check if model exists
91  if (!fModel)
92  return 0;
93 
94  // clear information
95  fParName.clear();
96  fParMin.clear();
97  fParMax.clear();
98  fMean.clear();
99  fMargMode.clear();
100  fGlobalMode.clear();
101  fQuantiles.clear();
102  fSmallInt.clear();
103  fRMS.clear();
104  fCorrCoeff.clear();
105 
106  // copy information from marginalized distributions
107  for (unsigned i = 0; i < fModel->GetNParameters(); ++i) {
108 
109  // copy parameter information
110  fParName.push_back( (fModel->GetParameter(i)->GetLatexName()) );
111  fParMin.push_back( fModel->GetParameter(i)->GetLowerLimit() );
112  fParMax.push_back( fModel->GetParameter(i)->GetUpperLimit() );
113 
114  // copy 1D marginalized information
115  BCH1D * bch1d_temp = fModel->GetMarginalized(i);
116  if (bch1d_temp) {
117  fFlagInfoMarg = true;
118  fMean.push_back( bch1d_temp->GetMean() );
119  fRMS.push_back( bch1d_temp->GetRMS() );
120  fMargMode.push_back( bch1d_temp->GetMode() );
121  for (unsigned j = 0; j < fSumProb.size(); ++j)
122  fQuantiles.push_back( bch1d_temp->GetQuantile( fSumProb.at(j) ) );
123  std::vector<double> intervals = bch1d_temp->GetSmallestIntervals();
124  int nintervals = int(intervals.size() / 5);
125  fSmallInt.push_back(nintervals);
126  fSmallInt.insert( fSmallInt.end(), intervals.begin(), intervals.end() );
127  }
128  else {
129  double tmpval = fModel->GetParameter(i)->GetUpperLimit() - fModel->GetParameter(i)->GetLowerLimit();
130  tmpval = fModel->GetParameter(i)->GetLowerLimit() - 2. * tmpval;
131  fMean.push_back( tmpval );
132  fRMS.push_back( tmpval );
133  fMargMode.push_back( tmpval );
134  for (unsigned j = 0; j < fSumProb.size(); ++j)
135  fQuantiles.push_back( tmpval );
136  fSmallInt.push_back( 0 );
137  }
138 
139  // copy 2D marginal information
140  for (unsigned j = 0; j < fModel->GetNParameters(); ++j) {
141  if (i == j)
142  fCorrCoeff.push_back(1.);
143  else {
144  BCH2D * bch2d_temp = fModel->GetMarginalized(i, j);
145  if ( bch2d_temp ) {
146  fFlagInfoMarg = true;
147  fCorrCoeff.push_back( bch2d_temp->GetHistogram()->GetCorrelationFactor() );
148  }
149  else
150  fCorrCoeff.push_back( 0. );
151  }
152  }
153 
154  // copy optimization information
155  if ((fModel->GetBestFitParameters()).size() > 0) {
156  fFlagInfoOpt = true;
157  fGlobalMode.push_back ( (fModel->GetBestFitParameters()).at(i) );
158  }
159  }
160 
161  // no error
162  return 1;
163 }
164 
165 // ---------------------------------------------------------
166 int BCSummaryTool::PrintParameterPlot(const char * filename)
167 {
168  // copy summary data
169  if (!CopySummaryData())
170  return 0;
171 
172  // get number of parameters and quantiles
173  int npar = fModel->GetNParameters();
174  int nquantiles = int( fSumProb.size() );
175 
176  // create histogram
177  TH1D * hist_axes = new TH1D(
178  TString::Format("hist_axes_par_%d",getNextIndex()),
179  ";;Scaled parameter range [a.u.]",npar, -0.5, npar-0.5);
180  hist_axes->SetStats(kFALSE);
181  for (int i = 0; i < npar; ++i)
182  hist_axes->GetXaxis()->SetBinLabel( i+1, fParName.at(i).c_str() );
183  hist_axes->GetXaxis()->SetLabelOffset(0.015);
184  hist_axes->GetXaxis()->SetLabelSize(0.06);
185  hist_axes->GetXaxis()->SetTickLength(0.0);
186  hist_axes->GetYaxis()->SetRangeUser(-0.1, 1.1);
187  hist_axes->GetYaxis()->SetTickLength(0.0);
188 
189  // create graphs
190  TGraphErrors * graph_quantiles = new TGraphErrors(npar*nquantiles);
191  graph_quantiles->SetMarkerSize(0);
192  graph_quantiles->SetLineColor(38);
193  graph_quantiles->SetLineStyle(2);
194 
195  TGraphErrors * graph_mean = new TGraphErrors(npar);
196  graph_mean->SetMarkerColor(kBlack);
197  graph_mean->SetMarkerStyle(20);
198 
199  TGraphErrors * graph_mode = new TGraphErrors(npar);
200  graph_mode->SetMarkerColor(kRed);
201  graph_mode->SetMarkerStyle(20);
202 
203  TGraphAsymmErrors * graph_intervals = new TGraphAsymmErrors(0);
204  graph_intervals->SetFillColor(kYellow);
205  graph_intervals->SetLineStyle(2);
206  graph_intervals->SetLineColor(kRed);
207  graph_intervals->SetMarkerSize(0);
208 
209  // fill graphs
210  int indexintervals = 0;
211 
212  // fill graph quantiles
213  if (fFlagInfoMarg) {
214  for (int i = 0; i < npar; ++i) {
215  for (int j = 0; j < nquantiles; ++j) {
216  graph_quantiles->SetPoint(i*nquantiles+j,double(i),
217  (fQuantiles.at(i*nquantiles+j) - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
218  graph_quantiles->SetPointError(i*nquantiles+j, 0.5, 0.0);
219  }
220  }
221  }
222 
223  // fill graph mean and rms
224  if (fFlagInfoMarg) {
225  for (int i = 0; i < npar; ++i) {
226  // fill graph mean
227  graph_mean->SetPoint(i, double(i), (fMean.at(i) - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
228  graph_mean->SetPointError(i, 0.0, fRMS.at(i)/(fParMax.at(i)-fParMin.at(i)));
229  }
230  }
231 
232  // fill graph mode
233  if (fFlagInfoOpt) {
234  for (int i = 0; i < npar; ++i)
235  graph_mode->SetPoint(i, double(i), (fGlobalMode.at(i) - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
236  }
237 
238  // fill graph smallest intervals
239  if (fFlagInfoMarg) {
240  for (int i = 0; i < npar; ++i) {
241  int nintervals = int(fSmallInt.at(indexintervals++));
242  for (int j = 0; j < nintervals; ++j) {
243  double xmin = fSmallInt.at(indexintervals++);
244  double xmax = fSmallInt.at(indexintervals++);
245  indexintervals++;
246  double xlocalmaxpos = fSmallInt.at(indexintervals++);
247  indexintervals++;
248  int npoints = graph_intervals->GetN();
249  graph_intervals->SetPoint(npoints,double(i),
250  (xlocalmaxpos - fParMin.at(i))/(fParMax.at(i)-fParMin.at(i)));
251  graph_intervals->SetPointError(npoints,0.5, 0.5,
252  (xlocalmaxpos - xmin)/(fParMax.at(i)-fParMin.at(i)),
253  (xmax - xlocalmaxpos)/(fParMax.at(i)-fParMin.at(i)));
254  }
255  }
256  }
257 
258  // create legend
259  TLegend * legend = new TLegend(0.15, 0.88, 0.85, 0.99);
260  legend->SetBorderSize(0);
261  legend->SetFillColor(0);
262 
263  // create latex
264  TLatex * latex = new TLatex();
265  latex->SetTextSize(0.02);
266 
267  // create lines
268  TLine * line_top = new TLine(-0.5, 1.0, npar-0.5, 1.0);
269  line_top->SetLineColor(kBlack);
270  line_top->SetLineStyle(1);
271  line_top->SetLineWidth(2);
272 
273  TLine * line_bot = new TLine(-0.5, 0.0, npar-0.5, 0.0);
274  line_bot->SetLineColor(kBlack);
275  line_bot->SetLineStyle(1);
276  line_bot->SetLineWidth(2);
277 
278  // print to file
279  TCanvas * c_par = new TCanvas(TString::Format("c_par_%d",getNextIndex()));
280  c_par->cd();
281  hist_axes->Draw();
282  line_top->Draw();
283  line_bot->Draw();
284  if (fFlagInfoMarg) {
285  graph_intervals->DrawClone("SAME2");
286  for (int i = 0; i < graph_intervals->GetN(); ++i)
287  graph_intervals->SetPointError(i, 0.5, 0.5, 0.0, 0.0);
288  graph_intervals->Draw("SAMEPZ");
289  graph_quantiles->Draw("SAMEPZ");
290  graph_mean->Draw("SAMEP");
291  legend->AddEntry(graph_quantiles, "Quantiles (5%, 10%, 16%, 50%, 84%, 90%, 95%)", "L");
292  legend->AddEntry(graph_mean, "Mean and RMS", "LEP");
293  legend->AddEntry(graph_intervals, "Smallest 68% intervals and local modes", "FL");
294  }
295  if (fFlagInfoOpt) {
296  graph_mode->Draw("SAMEP");
297  legend->AddEntry(graph_mode, "Global mode", "P");
298  }
299  for (int i = 0; i < npar;++i) {
300  // latex->DrawLatex(double(i)-0.1, 0.010, Form("%+3.3f", fParMin.at(i)));
301  // latex->DrawLatex(double(i)-0.1, 0.965, Form("%+3.3f", fParMax.at(i)));
302  latex->DrawLatex(double(i)-0.1, 0.010-0.07, Form("%+3.3g", fParMin.at(i)));
303  latex->DrawLatex(double(i)-0.1, 0.965+0.07, Form("%+3.3g", fParMax.at(i)));
304  }
305  latex->SetNDC();
306  latex->DrawLatex(0.9, 0.175, "Par. min.");
307  latex->DrawLatex(0.9, 0.83, "Par. max.");
308  legend->Draw("SAME");
309  gPad->RedrawAxis();
310  c_par->Print(filename);
311 
312  // no error
313  return 1;
314 }
315 
316 // ---------------------------------------------------------
317 int BCSummaryTool::PrintCorrelationMatrix(const char * filename)
318 {
319  // copy summary data
320  if (!CopySummaryData())
321  return 0;
322 
323  // check if marginalized information is there
324  if (!fFlagInfoMarg)
325  return 0;
326 
327  // get number of parameters
328  int npar = fModel->GetNParameters();
329 
330  // create histogram
331  TH2D * hist_corr = new TH2D(
332  TString::Format("hist_corr_%d",getNextIndex()),
333  ";;",npar, -0.5, npar-0.5,npar, -0.5, npar-0.5);
334  hist_corr->SetStats(kFALSE);
335  hist_corr->GetXaxis()->SetTickLength(0.0);
336  hist_corr->GetYaxis()->SetTickLength(0.0);
337  hist_corr->GetZaxis()->SetRangeUser(-1.0, 1.0);
338 
339  for (int i = 0; i < npar; ++i) {
340  hist_corr->GetXaxis()->SetLabelSize(0.06);
341  hist_corr->GetYaxis()->SetLabelSize(0.06);
342  if (npar < 5) {
343  hist_corr->GetXaxis()->SetBinLabel( i+1, fParName.at(i).c_str() );
344  hist_corr->GetYaxis()->SetBinLabel( npar-i, fParName.at(i).c_str() );
345  }
346  else {
347  hist_corr->GetXaxis()->SetBinLabel( i+1, TString::Format("%d",i) );
348  hist_corr->GetYaxis()->SetBinLabel( npar-i, TString::Format("%d",i) );
349  }
350  }
351 
352  // fill plot
353  for (int i = 0; i < npar; ++i)
354  for (int j = 0; j < npar; ++j) {
355  int index = i * npar + j;
356  double corr = fCorrCoeff.at(index);
357  hist_corr->SetBinContent(i+1, npar-j, corr);
358  }
359 
360  // print to file
361  TCanvas * c_corr = new TCanvas(TString::Format("c_corr_matrix_%d",getNextIndex()));
362  c_corr->cd();
363  hist_corr->Draw("colz text");
364 
365  TF1 * f = new TF1("fUp","x",-0.5,npar-0.5);
366  TGaxis * A1 = new TGaxis(-0.5,npar-0.5,npar-0.5,npar-0.5,"fUp",100,"-");
367  A1->ImportAxisAttributes(hist_corr->GetXaxis());
368  A1->Draw();
369 
370  // redraw the histogram to overlay thetop axis tick marks since
371  // we don't know how to make them disappear
372  hist_corr->GetXaxis()->SetLabelSize(0.);
373  hist_corr->Draw("colz text same");
374 
375  for (int i = 0; i < npar; ++i)
376  for (int j = 0; j < npar; ++j) {
377  BCH2D * bch2d_temp = fModel->GetMarginalized(fModel->GetParameter(i),fModel->GetParameter(j));
378  if ( bch2d_temp || i==j )
379  continue;
380 
381  TBox * bempty = new TBox(
382  hist_corr->GetXaxis()->GetBinLowEdge(i+1),
383  hist_corr->GetYaxis()->GetBinLowEdge(npar-j),
384  hist_corr->GetXaxis()->GetBinLowEdge(i+2),
385  hist_corr->GetYaxis()->GetBinLowEdge(npar-j+1)
386  );
387  bempty->SetLineStyle(0);
388  bempty->SetLineWidth(0);
389  bempty->SetFillColor(kWhite);
390  bempty->Draw();
391  }
392 
393  // redraw top and right axes
394  TLine * lA1 = new TLine(-0.5,npar-0.5,npar-0.5,npar-0.5);
395  lA1->Draw("same");
396  TLine * lA2 = new TLine(npar-0.5,npar-0.5,npar-0.5,-0.5);
397  lA2->Draw("same");
398 
399  gPad->RedrawAxis();
400  c_corr->Print(filename);
401 
402  delete f;
403  delete A1;
404  delete lA1;
405  delete lA2;
406  delete hist_corr;
407  delete c_corr;
408 
409  // no error
410  return 1;
411 }
412 
413 // ---------------------------------------------------------
414 int BCSummaryTool::PrintCorrelationPlot(const char * filename)
415 {
416  // get number of parameters
417  int npar = fModel->GetNParameters();
418 
419  TCanvas * c = new TCanvas(TString::Format("c_corr_%d",getNextIndex()));
420  c->cd();
421 
422  double margin = .1;
423  double padsize = (1. - 2.*margin) / (double)npar;
424 
425  // array with pads holding the histograms
426  TPad ** pad = new TPad*[npar*npar];
427 
428  // position of pads
429  double xlow, xup, ylow, yup;
430  double marginleft, marginright, margintop, marginbottom;
431  marginleft=marginright=margintop=marginbottom=0.01;
432 
433  // rotate label to avoid overlap with long parameter names
434  size_t maxlength = 0;
435  for (int i = 0; i < npar ; ++i) {
436  maxlength = std::max(maxlength, fModel->GetParameter(i)->GetName().length());
437  }
438 
439  double rotation = 0;
440  short xalignment = 22;
441  short yalignment = 22;
442 
443  if (maxlength > 20) {
444  rotation = 10;
445  xalignment = 12;
446  yalignment = 32;
447  }
448 
449  // drawing all histograms
450  for (int i=npar-1;i>=0;--i) {
451  xlow = (double)i*padsize + margin;
452  xup = xlow+padsize;
453 
454  for (int j=i;j<=npar-1;j++) {
455  yup = 1. - (double)j*padsize - margin;
456  ylow = yup-padsize;
457 
458  // preparing the pad
459  int ipad=i*npar+j;
460  pad[ipad]=new TPad(TString::Format("pad_%d_%d",ipad,getNextIndex()), "", xlow, ylow, xup, yup);
461  pad[ipad]->SetTopMargin(margintop);
462  pad[ipad]->SetBottomMargin(marginbottom);
463  pad[ipad]->SetLeftMargin(marginleft);
464  pad[ipad]->SetRightMargin(marginright);
465  pad[ipad]->SetFillColor(kWhite);
466 
467  pad[ipad]->Draw();
468  pad[ipad]->cd();
469 
470  // get the histogram
471  TH1 * hh = 0;
472  BCH1D * bh1 = 0;
473  BCH2D * bh2 = 0;
474  if(i==j) {
475  bh1 = fModel->GetMarginalized(fModel->GetParameter(i));
476  if (bh1)
477  hh = (TH1D*)bh1->GetHistogram()->Clone();
478  }
479  else {
480  bh2 = fModel->GetMarginalized(fModel->GetParameter(i),fModel->GetParameter(j));
481  if (bh2)
482  hh = (TH2D*)bh2->GetHistogram()->Clone();
483  }
484 
485  // if the histogram is not available, draw N/A
486  if (!hh) {
487  pad[ipad]->SetFillColor(kGray);
488  TBox * box = new TBox(marginleft,marginbottom,1.-marginright,1.-margintop);
489  box->SetLineWidth(1);
490  box->SetLineColor(kGray+1);
491  box->SetFillColor(kWhite);
492  box->Draw();
493 
494  TText * text_na = new TText(.5,.5,"N/A");
495  text_na->SetTextFont(42);
496  text_na->SetTextAlign(22);
497  text_na->SetTextSize(.8/(double)npar);
498  text_na->SetTextColor(kGray+1);
499  text_na->Draw();
500  }
501  // otherwise draw the histogram
502  else {
503  if(i==j) {
504  bh1->Draw("BTsiB3CS1D0");
505  }
506  else {
507  bh2->Draw("BTfB3CS1nL");
508  }
509 
510  hh->GetXaxis()->SetLabelOffset(5500);
511  hh->GetYaxis()->SetLabelOffset(5500);
512  hh->GetXaxis()->SetTitleSize(10.00);
513  hh->GetYaxis()->SetTitleSize(10.00);
514  }
515 
516  c->cd();
517 
518  // draw axis label
519  double labelsize = .8/(double)npar/10.;
520  double xtext, ytext;
521 
522  // y axis
523  if(i==0) {
524  TLatex * label = new TLatex();
525  label->SetTextFont(62);
526  label->SetTextSize(labelsize);
527  label->SetTextAlign(yalignment);
528  label->SetNDC();
529 
530  label->SetTextAngle(90 - rotation);
531 
532  xtext = margin * (1. - 8. * labelsize);
533  ytext = yup - padsize / 2.;
534 
535  label->DrawLatex(xtext,ytext,fModel->GetParameter(j)->GetLatexName().c_str());
536  }
537 
538  // x axis
539  if(j==npar-1) {
540  TLatex * label = new TLatex();
541  label->SetTextFont(62);
542  label->SetTextSize(labelsize);
543  label->SetTextAlign(xalignment);
544  label->SetNDC();
545 
546  label->SetTextAngle(360 - rotation);
547 
548  xtext = xlow + padsize / 2.;
549  ytext = margin * (1. - 6. * labelsize);
550 
551  label->DrawLatex(xtext,ytext, fModel->GetParameter(i)->GetLatexName().c_str());
552  }
553  }
554  }
555 
556  gPad->RedrawAxis();
557  c->Print(filename);
558 
559  return 1;
560 }
561 
562 // ---------------------------------------------------------
563 int BCSummaryTool::PrintKnowledgeUpdatePlot1D(int index, const char * filename, std::string options_post, std::string options_prior)
564 {
565  // perform analysis
566  CalculatePriorModel();
567 
568  // create canvas
569  TCanvas * c = new TCanvas();
570  c->cd();
571 
572  // draw
573  DrawKnowledgeUpdatePlot1D(index, options_post, options_prior);
574 
575  // print
576  c->Print(filename);
577 
578  // no error
579  return 1;
580 }
581 
582 // ---------------------------------------------------------
583 int BCSummaryTool::DrawKnowledgeUpdatePlot1D(int index, std::string options_post, std::string options_prior)
584 {
585  // option flags
586  bool flag_slice_post = false;
587  bool flag_slice_prior = false;
588 
589  // check content of options string
590  if (options_post.find("slice") < options_post.size()) {
591  flag_slice_post = true;
592  }
593  if (options_prior.find("slice") < options_prior.size()) {
594  flag_slice_prior = true;
595  }
596 
597  // prepare legend
598  TLegend* legend = new TLegend();
599  legend->SetBorderSize(0);
600  legend->SetFillColor(kWhite);
601  legend->SetTextAlign(12);
602  legend->SetTextFont(62);
603  legend->SetTextSize(0.03);
604 
605  // get histograms;
606  const BCParameter * par = fModel->GetParameter(index);
607  BCH1D* hist_prior = fPriorModel->GetMarginalized(par);
608  BCH1D* hist_posterior = 0;
609 
610  if (flag_slice_prior && fPriorModel->GetNParameters()==2) {
611  if (index == 0) {
612  TH1D* hist = fPriorModel->GetSlice(fPriorModel->GetParameter(0),fPriorModel->GetParameter(1))->GetHistogram()->ProjectionX(Form("projx_%i",BCLog::GetHIndex()));
613  hist->Scale(1.0/hist->Integral("width"));
614  for (int i = 1; i <= hist_prior->GetHistogram()->GetNbinsX(); ++i)
615  hist_prior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
616  }
617  else {
618  TH1D* hist = fPriorModel->GetSlice(fPriorModel->GetParameter(0),fPriorModel->GetParameter(1))->GetHistogram()->ProjectionY(Form("projy_%i",BCLog::GetHIndex()));
619  hist->Scale(1.0/hist->Integral("width"));
620  for (int i = 1; i <= hist_prior->GetHistogram()->GetNbinsX(); ++i)
621  hist_prior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
622  }
623  hist_prior->GetHistogram()->SetStats(kFALSE);
624  }
625  else if (flag_slice_prior && fPriorModel->GetNParameters()==1) {
626  hist_prior = fPriorModel->GetSlice(par);
627  hist_prior->GetHistogram()->SetStats(kFALSE);
628  }
629  // if marginal doesn't exist, skip ahead
630  if ( !hist_prior)
631  return 0;
632 
633  hist_prior->GetHistogram()->SetLineColor(kRed);
634 
635  hist_posterior = fModel->GetMarginalized(par);
636  if (flag_slice_post && fModel->GetNParameters()==2) {
637  if (index == 0) {
638  TH1D* hist = fModel->GetSlice(fModel->GetParameter(0),fModel->GetParameter(1))->GetHistogram()->ProjectionX(Form("projx_%i",BCLog::GetHIndex()));
639  hist->Scale(1.0/hist->Integral("width"));
640  for (int i = 1; i <= hist_posterior->GetHistogram()->GetNbinsX(); ++i)
641  hist_posterior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
642  }
643  else {
644  TH1D* hist = fModel->GetSlice(fModel->GetParameter(0),fModel->GetParameter(1))->GetHistogram()->ProjectionY(Form("projy_%i",BCLog::GetHIndex()));
645  hist->Scale(1.0/hist->Integral("width"));
646  for (int i = 1; i <= hist_posterior->GetHistogram()->GetNbinsX(); ++i)
647  hist_posterior->GetHistogram()->SetBinContent(i, hist->GetBinContent(i));
648  }
649  hist_posterior->GetHistogram()->SetStats(kFALSE);
650  }
651  else if (flag_slice_post && fModel->GetNParameters()==1) {
652  hist_posterior = fModel->GetSlice(par);
653  hist_posterior->GetHistogram()->SetStats(kFALSE);
654  }
655 
656  // if marginal doesn't exist, skip ahead
657  if ( !hist_posterior)
658  return 0;
659 
660  legend->AddEntry(hist_prior->GetHistogram(), "prior", "L");
661  legend->AddEntry(hist_posterior->GetHistogram(), "posterior", "L");
662 
663  // scale histograms
664  hist_posterior->GetHistogram()->Scale(1./hist_posterior->GetHistogram()->Integral("width"));
665  hist_prior->GetHistogram()->Scale(1.0/hist_prior->GetHistogram()->Integral("width"));
666 
667  // get maximum
668  double max_prior = hist_prior->GetHistogram()->GetMaximum();
669  double max_posterior = hist_posterior->GetHistogram()->GetMaximum();
670  double maxy = 1.1 * TMath::Max(max_prior, max_posterior);
671 
672  double height = 0.03*legend->GetNRows();
673 
674  // plot
675  hist_prior->GetHistogram()->GetXaxis()->SetNdivisions(508);
676  hist_posterior->GetHistogram()->GetXaxis()->SetNdivisions(508);
677 
678  hist_prior->Draw(options_prior);
679  hist_posterior->Draw(std::string(options_post+"same").c_str());
680 
681  // scale axes
682  hist_prior->GetHistogram()->GetYaxis()->SetRangeUser(0.0, maxy);
683  hist_posterior->GetHistogram()->GetYaxis()->SetRangeUser(0.0, maxy);
684 
685  gPad->SetTopMargin(0.02);
686  double xlegend1 = gPad->GetLeftMargin() + 0.10 * (1.0 - gPad->GetRightMargin() - gPad->GetLeftMargin());
687  double xlegend2 = 1.0-gPad->GetRightMargin();
688  double ylegend1 = 1.-gPad->GetTopMargin()-height;
689  double ylegend2 = 1.-gPad->GetTopMargin();
690 
691  // place legend on top of histogram
692  legend->SetX1NDC(xlegend1);
693  legend->SetX2NDC(xlegend2);
694  legend->SetY1NDC(ylegend1);
695  legend->SetY2NDC(ylegend2);
696 
697  // draw legend
698  legend->Draw();
699 
700  // rescale top margin
701  gPad->SetTopMargin(1.-ylegend1+0.01);
702 
703  gPad->RedrawAxis();
704 
705  return 1;
706 }
707 
708 // ---------------------------------------------------------
709 int BCSummaryTool::PrintKnowledgeUpdatePlots(const char * filename, std::string options)
710 {
711  // perform analysis
712  CalculatePriorModel();
713 
714  // option flags
715  bool flag_slice = false;
716 
717  // check content of options string
718  if (options.find("slice") < options.size()) {
719  flag_slice = true;
720  }
721 
722  std::string file(filename);
723 
724  // check if file extension is pdf
725  if ( (file.find_last_of(".") != std::string::npos) &&
726  (file.substr(file.find_last_of(".")+1) == "pdf") ) {
727  ; // it's a PDF file
728 
729  }
730  else if ( (file.find_last_of(".") != std::string::npos) &&
731  (file.substr(file.find_last_of(".")+1) == "ps") ) {
732  ; // it's a PS file
733  }
734  else {
735  ; // make it a PDF file
736  file += ".pdf";
737  }
738 
739  // create canvas and prepare postscript
740  TCanvas * c = new TCanvas(TString::Format("c_%d",getNextIndex()));
741  c->cd();
742 
743  // loop over all parameters and draw 1D plots
744  int npar = fModel->GetNParameters();
745  c->Print(std::string(file + "[").c_str());
746  for (int i = 0; i < npar; ++i) {
747  if ( !DrawKnowledgeUpdatePlot1D(i, options, options))
748  continue;
749  c->Print(file.c_str());
750  }
751 
752  // create legend
753  TLegend * legend2d = new TLegend();
754  legend2d->SetBorderSize(0);
755  legend2d->SetFillColor(0);
756  legend2d->SetTextAlign(12);
757  legend2d->SetTextFont(62);
758  legend2d->SetTextSize(0.03);
759 
760  // create markers and arrows
761  TMarker * marker_prior = new TMarker();
762  marker_prior->SetMarkerStyle(24);
763  marker_prior->SetMarkerColor(kRed);
764 
765  TMarker * marker_posterior = new TMarker();
766  marker_posterior->SetMarkerStyle(24);
767  marker_posterior->SetMarkerColor(kBlack);
768 
769  TArrow * arrow = new TArrow();
770  arrow->SetArrowSize(0.02);
771  arrow->SetLineColor(kBlue);
772  // arrow->SetLineStyle(2);
773 
774  // loop over all parameters
775  for (int i = 0; i < npar; ++i) {
776  for (int j = 0; j < i; ++j) {
777  c->cd();
778 
779  // get parameters
780  const BCParameter * par1 = fModel->GetParameter(i);
781  const BCParameter * par2 = fModel->GetParameter(j);
782 
783  // get 2-d histograms
784  BCH2D* bch2d_2dprior = 0;
785  BCH2D* bch2d_2dposterior = 0;
786  if (flag_slice && npar == 2) {
787  bch2d_2dprior = fPriorModel->GetSlice(par2, par1);
788  bch2d_2dposterior = fModel->GetSlice(par2, par1);
789  }
790  else {
791  bch2d_2dprior = fPriorModel->GetMarginalized(par1, par2);
792  bch2d_2dposterior = fModel->GetMarginalized(par1, par2);
793  }
794 
795  // can't draw anything
796  if ( !bch2d_2dprior || !bch2d_2dposterior)
797  continue;
798 
799  // get histograms
800  TH2D* hist_2dprior = bch2d_2dprior->GetHistogram();
801  hist_2dprior->SetLineColor(kRed);
802  TH2D* hist_2dposterior = bch2d_2dposterior->GetHistogram();
803 
804  // scale histograms
805  hist_2dprior->Scale(1.0/hist_2dprior->Integral("width"));
806  hist_2dposterior->Scale(1.0/hist_2dposterior->Integral("width"));
807 
808  // calculate contours
809  bch2d_2dprior->CalculateIntegratedHistogram();
810  bch2d_2dposterior->CalculateIntegratedHistogram();
811 
812  double level[1] = {bch2d_2dprior->GetLevel(0.32)};
813  hist_2dprior->SetContour(1, level);
814  hist_2dprior->Draw("CONT3");
815  level[0] = bch2d_2dposterior->GetLevel(0.32);
816  hist_2dposterior->SetContour(1, level);
817  hist_2dposterior->Draw("CONT3 SAME");
818 
819  std::vector<double> mode_prior = fPriorModel->GetBestFitParameters();
820  std::vector<double> mode_posterior = fModel->GetBestFitParameters();
821 
822  marker_prior->DrawMarker(mode_prior.at(j), mode_prior.at(i));
823  marker_posterior->DrawMarker(mode_posterior.at(j), mode_posterior.at(i));
824  arrow->DrawArrow(mode_prior.at(j), mode_prior.at(i), mode_posterior.at(j), mode_posterior.at(i));
825 
826  if (i==1 && j == 0) {
827  legend2d->AddEntry(hist_2dprior, "smallest 68% interval(s) of prior", "L");
828  legend2d->AddEntry(hist_2dposterior, "smallest 68% interval(s) of posterior", "L");
829  legend2d->AddEntry(marker_prior, "prior mode", "P");
830  legend2d->AddEntry(marker_posterior, "posterior mode", "P");
831  legend2d->AddEntry(arrow, "change in mode", "L");
832  }
833  gPad->SetTopMargin(0.02);
834 
835  double height = 0.03*legend2d->GetNRows();
836 
837  double xlegend1 = gPad->GetLeftMargin();
838  double xlegend2 = 1.0-gPad->GetRightMargin();
839  double ylegend1 = 1.-gPad->GetTopMargin()-height;
840  double ylegend2 = 1.-gPad->GetTopMargin();
841 
842  // place legend on top of histogram
843  legend2d->SetX1NDC(xlegend1);
844  legend2d->SetX2NDC(xlegend2);
845  legend2d->SetY1NDC(ylegend1);
846  legend2d->SetY2NDC(ylegend2);
847 
848  legend2d->Draw();
849 
850  gPad->SetTopMargin(1.-ylegend1+0.01);
851 
852  gPad->RedrawAxis();
853  c->Print(file.c_str());
854  }
855  }
856 
857  // close output
858  c->Print(std::string(file + "]").c_str());
859  c->Update();
860 
861  // free memory
862  delete legend2d;
863  delete marker_prior;
864  delete marker_posterior;
865  delete arrow;
866  delete c;
867 
868  // no error
869  return 1;
870 }
871 
872 // ---------------------------------------------------------
873 int BCSummaryTool::PrintParameterLatex(const char * filename)
874 {
875  // open file
876  std::ofstream ofi(filename);
877  ofi.precision(3);
878 
879  // check if file is open
880  if(!ofi.is_open()) {
881  std::cerr << "Couldn't open file " << filename <<std::endl;
882  return 0;
883  }
884 
885  // get number of parameters and quantiles
886  int npar = fModel->GetNParameters();
887 
888  // print table
889  ofi
890  << "\\documentclass[11pt, a4paper]{article}" << std::endl
891  << std::endl
892  << "\\begin{document}" << std::endl
893  << std::endl
894  << "\\begin{table}[ht!]" << std::endl
895  << "\\begin{center}" << std::endl
896  <<"\\begin{tabular}{llllllll}" << std::endl
897  << "\\hline" << std::endl
898  << "Parameter & Mean & RMS & Gl. mode & Mode & Median & 16\\% quant. & 84\\% quant. \\\\" << std::endl
899  << "\\hline" << std::endl;
900 
901  for (int i = 0; i < npar; ++i) {
902  const BCParameter * par = fModel->GetParameter(i);
903  BCH1D * bch1d = fModel->GetMarginalized(par);
904  ofi
905  << par->GetName() << " & "
906  << bch1d->GetMean() << " & "
907  << bch1d->GetRMS() << " & "
908  << fModel->GetBestFitParameters().at(i) << " & "
909  << bch1d->GetMode() << " & "
910  << bch1d->GetMedian() << " & "
911  << bch1d->GetQuantile(0.16) << " & "
912  << bch1d->GetQuantile(0.84) << " \\\\" << std::endl;
913  }
914  ofi
915  << "\\hline" << std::endl
916  << "\\end{tabular}" << std::endl
917  << "\\caption{Summary of the parameter estimates.}" << std::endl
918  << "\\end{center}" << std::endl
919  << "\\end{table}" << std::endl
920  << std::endl
921  << "\\end{document}" << std::endl;
922 
923  // close file
924  ofi.close();
925 
926  // no error
927  return 1;
928 }
929 
930 // ---------------------------------------------------------
931 int BCSummaryTool::CalculatePriorModel()
932 {
933  // create new prior model
934  delete fPriorModel;
935 
936  fPriorModel = new BCSummaryPriorModel();
937 
938  // set model
939  fPriorModel->SetModel(fModel);
940 
941  // perform marginalization
942  fPriorModel->MarginalizeAll();
943 
944  // perform minimization
945  fPriorModel->FindMode( fPriorModel->GetBestFitParameters() );
946 
947  // no error
948  return 1;
949 }
BCSummaryPriorModel
A helper class for the BCSummaryTool.
Definition: BCSummaryPriorModel.h:29
BCH2D::Draw
void Draw(std::string options="BTfB3CS1meangmodelmode", std::vector< double > intervals=std::vector< double >(0))
Definition: BCH2D.cxx:171
BCSummaryTool::PrintCorrelationMatrix
int PrintCorrelationMatrix(const char *filename="matrix.pdf")
Definition: BCSummaryTool.cxx:317
BCH1D::GetRMS
double GetRMS()
Definition: BCH1D.h:87
BCSummaryTool::DrawKnowledgeUpdatePlot1D
int DrawKnowledgeUpdatePlot1D(int index, std::string options_post="", std::string options_prior="")
Definition: BCSummaryTool.cxx:583
BCModel
The base class for all user-defined models.
Definition: BCModel.h:50
BCH2D
A class for handling 2D distributions.
Definition: BCH2D.h:37
BCH2D::GetLevel
double GetLevel(double p)
Definition: BCH2D.cxx:599
BCH1D::GetMean
double GetMean()
Definition: BCH1D.h:58
BCSummaryTool::PrintKnowledgeUpdatePlots
int PrintKnowledgeUpdatePlots(const char *filename="update.pdf", std::string options="")
Definition: BCSummaryTool.cxx:709
BCSummaryTool::CalculatePriorModel
int CalculatePriorModel()
Definition: BCSummaryTool.cxx:931
BCH2D::GetHistogram
TH2D * GetHistogram()
Definition: BCH2D.h:60
BCSummaryPriorModel::SetModel
void SetModel(BCModel *model)
Definition: BCSummaryPriorModel.cxx:38
BCLog::GetHIndex
static int GetHIndex()
Definition: BCLog.h:164
BCH1D::GetMedian
double GetMedian()
Definition: BCH1D.h:67
BCSummaryTool::PrintParameterPlot
int PrintParameterPlot(const char *filename="parameters.pdf")
Definition: BCSummaryTool.cxx:166
BCParameter
A class representing a parameter of a model.
Definition: BCParameter.h:28
BCH1D::GetMode
double GetMode()
Definition: BCH1D.cxx:52
BCSummaryTool::PrintCorrelationPlot
int PrintCorrelationPlot(const char *filename="correlation.pdf")
Definition: BCSummaryTool.cxx:414
BCSummaryTool::PrintParameterLatex
int PrintParameterLatex(const char *filename)
Definition: BCSummaryTool.cxx:873
BCH1D::GetHistogram
TH1D * GetHistogram()
Definition: BCH1D.h:53
BCH1D
A class for handling 1D distributions.
Definition: BCH1D.h:31
BCH1D::GetSmallestIntervals
std::vector< double > GetSmallestIntervals(double content=0.68)
Definition: BCH1D.cxx:1008
BCSummaryTool::PrintKnowledgeUpdatePlot1D
int PrintKnowledgeUpdatePlot1D(int index, const char *filename, std::string options_post="", std::string options_prior="")
Definition: BCSummaryTool.cxx:563
BCH1D::GetQuantile
double GetQuantile(double probablity)
Definition: BCH1D.cxx:58
BCParameter::GetName
const std::string & GetName() const
Definition: BCParameter.h:54
BCH2D::CalculateIntegratedHistogram
void CalculateIntegratedHistogram()
Definition: BCH2D.cxx:564
BCModel::GetName
const std::string & GetName() const
Definition: BCModel.h:85
BCH1D::Draw
void Draw(std::string options="BTsiB3CS1D0Lmeanmode", std::vector< double > intervals=std::vector< double >(0))
Definition: BCH1D.cxx:231