regression_averagedevs.cxx

Go to the documentation of this file.

#include <limits>

#include "TMVA/regression_averagedevs.h"

#include "TLatex.h"
#include "TGraphErrors.h"
#include "TFrame.h"

/*
  this macro plots the quadratic deviation of the estimated from the target value, averaged over the first nevt events in test sample (all if Nevt=-1)
  a) normal average
  b) truncated average, using best 90%
  created January 2009, Eckhard von Toerne, University of Bonn, Germany
*/

void TMVA::regression_averagedevs(TString dataset,TString fin, Int_t Nevt, Bool_t useTMVAStyle )
{
   bool debug=false;
   if (Nevt <0)  Nevt=1000000;
   TMVAGlob::Initialize( useTMVAStyle );
   // checks if file with name "fin" is already open, and if not opens one
   TFile* file = TMVAGlob::OpenFile( fin );
   TList jobDirList;
   TMVAGlob::GetListOfJobs((TFile*)file->GetDirectory(dataset.Data()),jobDirList);
   if (jobDirList.GetSize()==0) {
      cout << "error could not find jobs" << endl;
      return;
   }

   Bool_t __PLOT_LOGO__  = kTRUE;
   Bool_t __SAVE_IMAGE__ = kTRUE;

   TDirectory* dir0 = (TDirectory*) (jobDirList.At(0));
   //TDirectory* dir0 = (TDirectory*) (file->Get("InputVariables_Id"));
   Int_t nTargets = TMVAGlob::GetNumberOfTargets( dir0);

   if (debug) cout << "found targets " << nTargets<<endl;
   TCanvas* c=0;
   for (Int_t itrgt = 0 ; itrgt < nTargets; itrgt++){
      if (debug) cout << "loop targets " << itrgt<<endl;
      TString xtit = "Method";
      TString ytit = "Average Quadratic Deviation";
      TString ftit = ytit + " versus " + xtit + Form(" for target %d",itrgt);
      c = new TCanvas( Form("c%d",itrgt), ftit , 50+20*itrgt, 10*itrgt, 750, 650 );

      // global style settings
      c->SetGrid();
      c->SetTickx(1);
      c->SetTicky(0);
      c->SetTopMargin(0.28);
      c->SetBottomMargin(0.1);

      TString hNameRef(Form("regression_average_devs_target%d",itrgt));

      const Int_t maxMethods = 100;
      //     const Int_t maxTargets = 100;
      Float_t m[4][maxMethods]; // h0 train-all, h1 train-90%, h2 test-all, h3 test-90%
      Float_t em[4][maxMethods];
      Float_t x[4][maxMethods];
      Float_t ex[4][maxMethods];

      TIter next(&jobDirList);
      Float_t mymax=0., mymin=std::numeric_limits<float>::max();
      TString mvaNames[maxMethods];
      TDirectory *jobDir;
      Int_t nMethods = 0;
      // loop over all methods
      while ( (jobDir = (TDirectory*)next()) ) {
         TString methodTitle;
         TMVAGlob::GetMethodTitle(methodTitle,jobDir);
         mvaNames[nMethods]=methodTitle;
         if (debug) cout << "--- Found directory for method: " << methodTitle << endl;
         TIter keyIt(jobDir->GetListOfKeys());
         TKey *histKey;
         while ( (histKey = (TKey*)keyIt()) ) {
            if (histKey->ReadObj()->InheritsFrom("TH1F") ){
               TString s(histKey->ReadObj()->GetName());
               if( !s.Contains("Quadr_Dev") ) continue;
               if( !s.Contains(Form("target_%d_",itrgt))) continue;
               Int_t ihist = 0 ;
               if( !s.Contains("best90perc") && s.Contains("train")) ihist=0;
               if( s.Contains("best90perc") && s.Contains("train")) ihist=1;
               if( !s.Contains("best90perc") && s.Contains("test")) ihist=2;
               if( s.Contains("best90perc") && s.Contains("test")) ihist=3;
               if (debug) cout <<"using histogram" << s << ", ihist="<<ihist<<endl;
               TH1F* h = (TH1F*) (histKey->ReadObj());
               m[ihist][nMethods] = sqrt(h->GetMean());
               em[ihist][nMethods] = h->GetRMS()/(sqrt(h->GetEntries())*2.*h->GetMean());
               x[ihist][nMethods] = nMethods+0.44+0.12*ihist;
               ex[ihist][nMethods] = 0.001;
               mymax=  m[ihist][nMethods] > mymax ? m[ihist][nMethods] : mymax;
               mymin=  m[ihist][nMethods] < mymin ? m[ihist][nMethods] : mymin;
               if (debug) cout << "m"<< ihist << "="<<m[ihist][nMethods]<<endl;
            }
         }
         nMethods++;
      }
      TH1F* haveragedevs= new TH1F(Form("haveragedevs%d",itrgt),ftit,nMethods,0.,nMethods);
      for (int i=0;i<nMethods;i++) haveragedevs->GetXaxis()->SetBinLabel(i+1, mvaNames[i]);
      haveragedevs->SetStats(0);
      TGraphErrors* graphTrainAv= new TGraphErrors(nMethods,x[0],m[0],ex[0],em[0]);
      TGraphErrors* graphTruncTrainAv= new TGraphErrors(nMethods,x[1],m[1],ex[1],em[1]);
      TGraphErrors* graphTestAv= new TGraphErrors(nMethods,x[2],m[2],ex[2],em[2]);
      TGraphErrors* graphTruncTestAv= new TGraphErrors(nMethods,x[3],m[3],ex[3],em[3]);

      Double_t xmax = 1.2 * mymax;
      Double_t xmin = 0.8 * mymin - (mymax - mymin)*0.05;
      Double_t xheader = 0.2;
      Double_t yheader = xmax*0.92;
      xmin = xmin > 0.? xmin : 0.;
      if (mymin > 1.e-20 && log10(mymax/mymin)>1.5){
         c->SetLogy();
         cout << "--- result differ significantly using log scale for display of regression results"<< endl;
         xmax = 1.5 * xmax;
         xmin = 0.75 * mymin;
         yheader = xmax*0.78;
      }
      Float_t x0L = 0.03,     y0H = 0.91;
      Float_t dxL = 0.457-x0L, dyH = 0.14;
      //     TLegend *legend = new TLegend( x0L, y0H-dyH, x0L+dxL, y0H , "Average Deviation = (#sum_{evts} (f_{MVA} - f_{target})^{2} )^{1/2}");
      TLegend *legend = new TLegend( x0L, y0H-dyH, x0L+dxL, y0H );
      legend->SetTextSize( 0.035 );
      legend->SetTextAlign(12);
      legend->SetMargin( 0.1 );

      TH1F *hr = c->DrawFrame(-1.,0.,nMethods+1, xmax);
      cout << endl;
      cout << "Training: Average Deviation between target " << itrgt <<" and estimate" << endl;
      cout << Form("%-15s%-15s%-15s", "Method","Average Dev.","trunc. Aver.(90%)") <<endl;
      for (int i=0;i<nMethods;i++){
         cout << Form("%-15s:%#10.3g%#10.3g",
                      (const char*)mvaNames[i], m[0][i],m[1][i])<<endl;
         //       cout << mvaNames[i] << "  " << m[0][i]<< "  "<< m[1][i]<<endl;
         hr->GetXaxis()->SetBinLabel(i+1," ");
      }
      cout << endl;
      cout << "Testing: Average Deviation between target " << itrgt <<" and estimate" << endl;
      cout << Form("%-15s%-15s%-15s", "Method","Average Dev.","trunc. Aver.(90%)") <<endl;
      for (int i=0;i<nMethods;i++){
         cout << Form("%-15s:%#10.3g%#10.3g",
                      (const char*)mvaNames[i], m[2][i],m[3][i])<<endl;
         //cout << mvaNames[i] << "  " << m[2][i]<< "  "<< m[3][i]<<endl;
      }

      haveragedevs->SetMinimum(xmin);
      haveragedevs->SetMaximum(xmax);
      haveragedevs->SetXTitle("Method");
      haveragedevs->SetYTitle("Deviation from target");
      haveragedevs->Draw();
      c->GetFrame()->SetFillColor(21);
      c->GetFrame()->SetBorderSize(12);
      graphTrainAv->SetMarkerSize(1.);
      graphTrainAv->SetMarkerColor(kBlue);
      graphTrainAv->SetMarkerStyle(25);
      graphTrainAv->Draw("P");

      graphTruncTrainAv->SetMarkerSize(1.);
      graphTruncTrainAv->SetMarkerColor(kBlack);
      graphTruncTrainAv->SetMarkerStyle(25);
      graphTruncTrainAv->Draw("P");

      graphTestAv->SetMarkerSize(1.);
      graphTestAv->SetMarkerColor(kBlue);
      graphTestAv->SetMarkerStyle(21);
      graphTestAv->Draw("P");

      graphTruncTestAv->SetMarkerSize(1.);
      graphTruncTestAv->SetMarkerColor(kBlack);
      graphTruncTestAv->SetMarkerStyle(21);
      graphTruncTestAv->Draw("P");
      legend->AddEntry(graphTrainAv,TString("Training Sample, Average Deviation"),"p");
      legend->AddEntry(graphTruncTrainAv,TString("Training Sample, truncated Average Dev. (best 90%)"),"p");
      legend->AddEntry(graphTestAv,TString("Test Sample, Average Deviation"),"p");
      legend->AddEntry(graphTruncTestAv,TString("Test Sample, truncated Average Dev. (best 90%)"),"p");

      legend->Draw();
      TLatex legHeader;
      legHeader.SetTextSize(0.035);
      legHeader.SetTextAlign(12);
      //legHeader.DrawLatex(x0L, y0H+0.01, "Average Deviation = (#sum (_{ } f_{MVA} - f_{target})^{2} )^{1/2}");
      legHeader.DrawLatex(xheader, yheader, "Average Deviation = (#sum (_{ } f_{MVA} - f_{target})^{2} )^{1/2}");
      // ============================================================

      if (__PLOT_LOGO__) TMVAGlob::plot_logo();
      // ============================================================

      c->Update();
      TString fname = dataset+"/plots/" + hNameRef;
      if (__SAVE_IMAGE__) TMVAGlob::imgconv( c, fname );
   } // end loop itrgt
   return;
}