KS: You validate stability of posterior covariance matrix, you set burn calc cov matrix increase burn calc again and compare. By performing such operation several hundred times we can check when matrix becomes stable.
435{
436
438
439 auto Canvas = std::make_unique<TCanvas>("Canvas", "Canvas", 0, 0, 1024, 1024);
440 Canvas->SetGrid();
441 gStyle->SetOptStat(0);
442 gStyle->SetOptTitle(0);
443 Canvas->SetTickx();
444 Canvas->SetTicky();
445 Canvas->SetBottomMargin(0.1f);
446 Canvas->SetTopMargin(0.05f);
447 Canvas->SetRightMargin(0.15f);
448 Canvas->SetLeftMargin(0.10f);
449
450
451 const int NRGBs = 10;
452 TColor::InitializeColors();
453 Double_t stops[NRGBs] = { 0.00, 0.10, 0.25, 0.35, 0.50, 0.60, 0.65, 0.75, 0.90, 1.00 };
454 Double_t red[NRGBs] = { 0.50, 1.00, 1.00, 0.25, 0.00, 0.10, 0.50, 1.00, 0.75, 0.55 };
455 Double_t green[NRGBs] = { 0.00, 0.25, 1.00, 0.25, 0.00, 0.60, 0.90, 1.00, 0.75, 0.75 };
456 Double_t blue[NRGBs] = { 0.00, 0.25, 1.00, 1.00, 0.50, 0.60, 0.90, 1.00, 0.05, 0.05 };
457 TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, 255);
458 gStyle->SetNumberContours(255);
459
460 std::string OutName = inputFile;
461 OutName = OutName.substr(0, OutName.find(".root"));
462 Canvas->Print(Form("Correlation_%s.pdf[", OutName.c_str()), "pdf");
463 Canvas->Print(Form("Covariance_%s.pdf[", OutName.c_str()), "pdf");
464
466 YAML::Node Settings = card_yaml["ProcessMCMC"];
467
468 const int entries = int(Processor->
GetnSteps());
469 const int NIntervals = GetFromManager<int>(Settings["NIntervals"], 5);
470 const int IntervalsSize = entries/NIntervals;
471
472 int BurnIn = 0;
473 MACH3LOG_INFO(
"Diagnosing matrices with entries={}, NIntervals={} and IntervalsSize={}", entries, NIntervals, IntervalsSize);
474
475 TMatrixDSym *Covariance = nullptr;
476 TMatrixDSym *Correlation = nullptr;
477
478 TH2D *CovariancePreviousHist = nullptr;
479 TH2D *CorrelationPreviousHist = nullptr;
480
481 TH2D *CovarianceHist = nullptr;
482 TH2D *CorrelationHist = nullptr;
483
484
487
490
491 delete Covariance;
492 Covariance = nullptr;
493 delete Correlation;
494 Correlation = nullptr;
495
496
497 for(int k = 1; k < NIntervals; ++k)
498 {
499 BurnIn = k*IntervalsSize;
503
506
507 TH2D *CovarianceDiff = static_cast<TH2D*>(CovarianceHist->Clone("Covariance_Ratio"));
508 TH2D *CorrelationDiff = static_cast<TH2D*>(CorrelationHist->Clone("Correlation_Ratio"));
509
510
511 #ifdef MULTITHREAD
512 #pragma omp parallel for
513 #endif
514 for (int j = 1; j < CovarianceDiff->GetXaxis()->GetNbins()+1; ++j)
515 {
516 for (int i = 1; i < CovarianceDiff->GetYaxis()->GetNbins()+1; ++i)
517 {
518 if( std::fabs (CovarianceDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CovariancePreviousHist->GetBinContent(j, i)) < 1.e-5)
519 {
520 CovarianceDiff->SetBinContent(j, i,
_UNDEF_);
521 CovariancePreviousHist->SetBinContent(j, i,
_UNDEF_);
522 }
523 if( std::fabs (CorrelationDiff->GetBinContent(j, i)) < 1.e-5 && std::fabs (CorrelationPreviousHist->GetBinContent(j, i)) < 1.e-5)
524 {
525 CorrelationDiff->SetBinContent(j, i,
_UNDEF_);
526 CorrelationPreviousHist->SetBinContent(j, i,
_UNDEF_);
527 }
528 }
529 }
530
531 CovarianceDiff->Divide(CovariancePreviousHist);
532 CorrelationDiff->Divide(CorrelationPreviousHist);
533
534
535 for (int j = 0; j < CovarianceDiff->GetXaxis()->GetNbins(); ++j)
536 {
537 TString Title = "";
538 double Prior = 1.0;
539 double PriorError = 1.0;
540
542
543 CovarianceDiff->GetXaxis()->SetBinLabel(j+1, Title);
544 CovarianceDiff->GetYaxis()->SetBinLabel(j+1, Title);
545 CorrelationDiff->GetXaxis()->SetBinLabel(j+1, Title);
546 CorrelationDiff->GetYaxis()->SetBinLabel(j+1, Title);
547 }
548 CovarianceDiff->GetXaxis()->SetLabelSize(0.015f);
549 CovarianceDiff->GetYaxis()->SetLabelSize(0.015f);
550 CorrelationDiff->GetXaxis()->SetLabelSize(0.015f);
551 CorrelationDiff->GetYaxis()->SetLabelSize(0.015f);
552
553 std::stringstream ss;
554 ss << "BCut_";
555 ss << BurnIn;
556 ss << "/";
557 ss << "BCut_";
558 ss << (k-1)*IntervalsSize;
559 std::string str = ss.str();
560
561 TString Title = "Cov " + str;
562 CovarianceDiff->GetZaxis()->SetTitle( Title );
563 Title = "Corr " + str;
564 CorrelationDiff->GetZaxis()->SetTitle(Title);
565
566 CovarianceDiff->SetMinimum(-2);
567 CovarianceDiff->SetMaximum(2);
568 CorrelationDiff->SetMinimum(-2);
569 CorrelationDiff->SetMaximum(2);
570
571 Canvas->cd();
572 CovarianceDiff->Draw("colz");
573 Canvas->Print(Form("Covariance_%s.pdf", OutName.c_str()), "pdf");
574
575 CorrelationDiff->Draw("colz");
576 Canvas->Print(Form("Correlation_%s.pdf", OutName.c_str()), "pdf");
577
578
579 delete CovariancePreviousHist;
580 CovariancePreviousHist = static_cast<TH2D*>(CovarianceHist->Clone());
581 delete CorrelationPreviousHist;
582 CorrelationPreviousHist = static_cast<TH2D*>(CorrelationHist->Clone());
583
584 delete CovarianceHist;
585 CovarianceHist = nullptr;
586 delete CorrelationHist;
587 CorrelationHist = nullptr;
588
589 delete CovarianceDiff;
590 delete CorrelationDiff;
591 delete Covariance;
592 Covariance = nullptr;
593 delete Correlation;
594 Correlation = nullptr;
595 }
596 Canvas->cd();
597 Canvas->Print(Form("Covariance_%s.pdf]", OutName.c_str()), "pdf");
598 Canvas->Print(Form("Correlation_%s.pdf]", OutName.c_str()), "pdf");
599
601 if(Covariance != nullptr) delete Covariance;
602 if(Correlation != nullptr) delete Correlation;
603 if(CovariancePreviousHist != nullptr) delete CovariancePreviousHist;
604 if(CorrelationPreviousHist != nullptr) delete CorrelationPreviousHist;
605 if(CovarianceHist != nullptr) delete CovarianceHist;
606 if(CorrelationHist != nullptr) delete CorrelationHist;
607}
TH2D * TMatrixIntoTH2D(TMatrixDSym *Matrix, const std::string &title)
KS: Convert TMatrix to TH2D, mostly useful for making fancy plots.
void GetNthParameter(const int param, double &Prior, double &PriorError, TString &Title) const
Get properties of parameter by passing it number.
void ResetHistograms()
Reset 2D posteriors, in case we would like to calculate in again with different BurnInCut.
void SetPrintToPDF(const bool PlotOrNot)
Long64_t GetnSteps()
Get Number of Steps that Chain has, for merged chains will not be the same nEntries.
void GetCovariance(TMatrixDSym *&Cov, TMatrixDSym *&Corr)
Get the post-fit covariances and correlations.
void SetStepCut(const std::string &Cuts)
Set the step cutting by string.