FrequentistSurface.cxx

Go to the documentation of this file.

#include "CAFAna/Fit/FrequentistSurface.h"

#include <iostream>
#include <functional>
#include <limits>

#include "TCanvas.h"
#include "TGraph.h"
#include "TH2F.h"
#include "TObjArray.h"
#include "TPad.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TKey.h"
#include "TVectorD.h"
#include "TObjString.h"
//#include "TCollection.h"

#include "CAFAna/Core/Var.h"
#include "CAFAna/Experiment/IExperiment.h"
#include "CAFAna/Fit/MinuitFitter.h"
#include "CAFAna/Core/LoadFromFile.h"
#include "CAFAna/Core/IFitVar.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/ThreadPool.h"
#include "CAFAna/Core/Utilities.h"

#include "OscLib/IOscCalc.h"

namespace ana
{
  //----------------------------------------------------------------------
  FrequentistSurface::FrequentistSurface(const IExperiment* expt,
                                         osc::IOscCalcAdjustable* calc,
                                         const FitAxis& xax,
                                         const FitAxis& yax,
                                         const std::vector<const IFitVar*>& profVars,
                                         const std::vector<const ISyst*>& profSysts,
                                         const SeedList& seedPts,
                                         const std::vector<SystShifts>& systSeedPts,
                                         bool parallel,
                                         MinuitFitter::FitOpts opts)
    : fParallel(parallel), fFitOpts(opts)
  {
    fLogX = xax.islog;
    fLogY = yax.islog;

    CreateHistograms(xax, yax, profVars, profSysts);

    for(const IFitVar* v: profVars) fSeedValues.push_back(v->GetValue( calc));

    FillSurface(expt, calc, xax.var, yax.var, profVars, profSysts, seedPts, systSeedPts);

    FindMinimum(expt, calc, xax.var, yax.var, profVars, profSysts, seedPts, systSeedPts);
  }

  //---------------------------------------------------------------------
  FrequentistSurface::~FrequentistSurface()
  {
  }

  //---------------------------------------------------------------------
  void FrequentistSurface::
  CreateHistograms(const FitAxis& xax, const FitAxis& yax,
                   const std::vector<const IFitVar*>& profVars,
                   const std::vector<const ISyst*>& profSysts)
  {
    fHist = ExpandedHistogram(";"+xax.var->LatexName()+";"+yax.var->LatexName(),
                              xax.nbins, xax.min, xax.max, xax.islog,
                              yax.nbins, yax.min, yax.max, yax.islog);

    for(unsigned int i = 0; i < profVars.size()+profSysts.size(); ++i){
      std::string title;
      if(i < profVars.size())
        title = profVars[i]->LatexName();
      else
        title = profSysts[i-profVars.size()]->LatexName();

      fProfHists.push_back(ExpandedHistogram(title+";"+xax.var->LatexName()+";"+yax.var->LatexName(),
                                             xax.nbins, xax.min, xax.max, xax.islog,
                                             yax.nbins, yax.min, yax.max, yax.islog));
    }

  }

  //---------------------------------------------------------------------
  std::string FrequentistSurface::
  ProgressBarTitle(const IFitVar* xvar, const IFitVar* yvar,
                   const std::vector<const IFitVar*>& profVars,
                   const std::vector<const ISyst*>& profSysts) const
  {
    std::string title = "Filling surface for "+yvar->ShortName()+" vs "+xvar->ShortName();

    if(!profVars.empty() || !profSysts.empty()){
      title += " (profiling ";

      for(const IFitVar* v: profVars) title += v->ShortName() + ", ";
      for(const ISyst* s: profSysts)  title += s->ShortName() + ", ";

      // Always have one superfluous ", " at the end
      title.resize(title.size()-2);
      title += ")";
    }

    return title;
  }

  //---------------------------------------------------------------------
  void FrequentistSurface::FillSurface(const IExperiment *expt,
                                       osc::IOscCalcAdjustable *calc,
                                       const IFitVar *xvar, const IFitVar *yvar,
                                       const std::vector<const IFitVar *> &profVars,
                                       const std::vector<const ISyst *> &profSysts,
                                       const SeedList& seedPts,
                                       const std::vector<SystShifts> &systSeedPts)
  {
    // Nothing created during surface filling belongs in a
    // directory. Unfortunately the local guards in Spectrum etc are racey when
    // run in parallel. But this should cover the whole lot safely.
    DontAddDirectory guard;

    const std::string progTitle = ProgressBarTitle(xvar, yvar, profVars, profSysts);

    Progress *prog = 0;
    // Difficult to keep a progress bar properly up to date when threaded
    if (!fParallel) prog = new Progress(progTitle);
    ThreadPool *pool = 0;

    if(fParallel){
      // A hack/workaround needed for parallel running:
      //
      // Give all the constituents of the Prediction a chance to do lazy
      // initialization, before they race themselves trying to do it in
      // parallel.
      expt->ChiSq(calc);

      pool = new ThreadPool;
      pool->ShowProgress(progTitle + TString::Format(" using %d threads", pool->NThreads()).Data());
    }

    const int Nx = fHist->GetNbinsX();
    const int Ny = fHist->GetNbinsY();

    // Fill bins in "random" order so that the progress bar is accurate
    // straight away instead of possibly being misled by whatever atypical
    // points we start with. This step is a prime which guarantees we get every
    // cell.
    int step = 7919;
    // Very unlikely (Nx or Ny is a multiple of step), but just to be safe.
    if ((Nx * Ny) % step == 0) step = 1;

    int bin = 0;
    int neval = 0;

    // Allow the surface to be parallelised across multiple jobs by splitting up
    // the full surface into patches, and only running bins that fall inside a
    // certain patch
    int first = 0, last = Nx * Ny;
    if (RunningOnGrid() && NumJobs()>1) {
      double stride = double(Nx*Ny) / double(NumJobs());
      int job = JobNumber();
      first = job * stride;
      last = (job+1) * stride;
    }

    do{
      if (neval < first) {
        ++neval;
        bin = (bin + step) % (Nx * Ny);
        continue;
      } else if (neval >= last) {
        break;
      }

      const int x = bin % Nx + 1;
      const int y = bin / Nx + 1;

      const double xv = BinCenterX(x);
      const double yv = BinCenterY(y);

      if (xvar->Penalty(xv, calc) > 1e-10)
      {
        std::cerr << "Warning! " << xvar->ShortName() << " = " << xv
                  << " has penalty of " << xvar->Penalty(xv, calc)
                  << " that could have been applied in surface. "
                  << "This should never happen." << std::endl;
      }
      if (yvar->Penalty(yv, calc) > 1e-10)
      {
        std::cerr << "Warning! " << yvar->ShortName() << " = " << yv
                  << " has penalty of " << yvar->Penalty(yv, calc)
                  << " that could have been applied in surface. "
                  << "This should never happen." << std::endl;
      }

      ThreadPool::func_t task = [=](){
        FillSurfacePoint(expt, calc,
                         xvar, xv, yvar, yv,
                         profVars, profSysts, seedPts, systSeedPts);
      };

      ++neval;
      if(fParallel){
        pool->AddTask(task);
      }
      else{
        task(); // Just do it straight away
        prog->SetProgress(double(neval-first) / double(last-first));
      }

      if (RunningOnGrid() && NumJobs()>1) fBinMask.push_back(bin);

      bin = (bin + step) % (Nx * Ny);
    } while (bin != 0);


    if(fParallel){
      pool->Finish();
      delete pool;
    }
    else{
      prog->Done();
      delete prog;
    }
  }

  //----------------------------------------------------------------------
  double FrequentistSurface::FillSurfacePoint(const IExperiment* expt,
                                              osc::IOscCalcAdjustable* calc,
                                              const IFitVar* xvar, double x,
                                              const IFitVar* yvar, double y,
                                              const std::vector<const IFitVar*>& profVars,
                                              const std::vector<const ISyst*>& profSysts,
                                              const SeedList& seedPts,
                                              const std::vector<SystShifts>& systSeedPts)
  {
    if(fParallel){
      // Need to take our own copy so that we don't get overwritten by someone
      // else's changes.
      calc = calc->Copy();
    }

    xvar->SetValue(calc, x);
    yvar->SetValue(calc, y);

    //Make sure that the profiled values of fitvars do not persist between steps.
    for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );

    expt->Reset();

    double chi;
    if(profVars.empty() && profSysts.empty()){
      chi = expt->ChiSq(calc);
    }
    else{
      MinuitFitter fitter(expt, profVars, profSysts);
      fitter.SetFitOpts(fFitOpts);
      SystShifts bestSysts;
      chi = fitter.Fit(calc, bestSysts, seedPts, systSeedPts, MinuitFitter::kQuiet)->EvalMetricVal();

      for(unsigned int i = 0; i < profVars.size(); ++i){
        fProfHists[i]->Fill(x, y, profVars[i]->GetValue(calc));
      }
      for(unsigned int j = 0; j < profSysts.size(); ++j){
        fProfHists[j+profVars.size()]->Fill(x, y, bestSysts.GetShift(profSysts[j]));
      }
    }

    fHist->Fill(x, y, chi);

    if(fParallel) delete calc;

    return chi;
  }


  //---------------------------------------------------------------------
  void FrequentistSurface::FindMinimum(const IExperiment* expt,
                                       osc::IOscCalcAdjustable* calc,
                                       const IFitVar* xvar, const IFitVar* yvar,
                                       const std::vector<const IFitVar*>& profVars,
                                       const std::vector<const ISyst*>& profSysts,
                                       const SeedList& seedPts,
                                       const std::vector<SystShifts>& systSeedPts)
  {
    // Location of the best minimum found from filled surface
    double minchi = 1e10;
    int minx = fHist->GetNbinsX()/2;
    int miny = fHist->GetNbinsY()/2;
    for(int x = 1; x <= fHist->GetNbinsX(); ++x){
      for(int y = 1; y <= fHist->GetNbinsY(); ++y){
        int bin = ((y-1)*fHist->GetNbinsX())+(x-1);
        const double chi = fHist->GetBinContent(x, y);
        if(RunningOnGrid() && NumJobs()>1
          && !std::count(fBinMask.begin(), fBinMask.end(), bin)) continue;
        if (chi < minchi){
          minchi = chi;
          minx = x;
          miny = y;
        }
      }
    }

    std::vector<const IFitVar*> allVars = {xvar, yvar};
    allVars.insert(allVars.end(), profVars.begin(), profVars.end());
    MinuitFitter fit(expt, allVars, profSysts);
    fit.SetFitOpts(fFitOpts);
    expt->Reset();
    // Seed from best grid point
    xvar->SetValue(calc, BinCenterX(minx));
    yvar->SetValue(calc, BinCenterY(miny));
    for(int i = 0; i < (int)fSeedValues.size(); ++i) profVars[i]->SetValue( calc, fSeedValues[i] );
    SystShifts systSeed = SystShifts::Nominal();
    fBestLikelihood = fit.Fit(calc, systSeed, seedPts)->EvalMetricVal();
    fBestFitX = xvar->GetValue(calc);
    fBestFitY = yvar->GetValue(calc);

    for(int x = 0; x < fHist->GetNbinsX()+2; ++x){
      for(int y = 0; y < fHist->GetNbinsY()+2; ++y){
        fHist->SetBinContent(x, y, fHist->GetBinContent(x, y)-fBestLikelihood);
      }
    }

    fHist->SetMinimum(0);
  }

  //----------------------------------------------------------------------
  void FrequentistSurface::SaveTo(TDirectory* dir, const std::string& name) const
  {
    TDirectory *tmp = gDirectory;

    dir = dir->mkdir(name.c_str()); // switch to subdir
    dir->cd();

    TObjString("FrequentistSurface").Write("type");

    ISurface::SaveToHelper(dir);

    TDirectory *profDir = dir->mkdir("profHists");
    int idx = 0;
    for (auto it: FrequentistSurface::fProfHists)
    {
      profDir->cd();
      it->Write(TString::Format("hist%d", idx++));
    }

    dir->Write();
    delete dir;

    tmp->cd();
  }

  //----------------------------------------------------------------------
  std::unique_ptr<FrequentistSurface> FrequentistSurface::LoadFrom(TDirectory* dir, const std::string& name)
  {
    dir = dir->GetDirectory(name.c_str()); // switch to subdir
    assert(dir);

    DontAddDirectory guard;

    TObjString *tag = (TObjString *) dir->Get("type");
    assert(tag);
    assert(tag->GetString() == "FrequentistSurface" || tag->GetString() == "Surface");
    delete tag;

    std::unique_ptr<FrequentistSurface> surf(new FrequentistSurface);
    ISurface::FillSurfObj(*surf, dir);

    for (std::size_t idx = 0; idx < surf->fSeedValues.size(); ++idx)
    {
      // Search for old "marg" name here too for backwards compatibility
      TH2 *h = (TH2 *) dir->Get(TString::Format("profHists/hist%lu", idx));
      if (h)
        surf->fProfHists.push_back(h);
      else
        surf->fProfHists.push_back((TH2 *) dir->Get(TString::Format("margHists/hist%lu", idx)));
    }

    const TObjString* logx = (TObjString*)dir->Get("logx");
    const TObjString* logy = (TObjString*)dir->Get("logy");
    // Tolerate missing log tags for backwards compatibility
    surf->fLogX = (logx && logx->GetString() == "yes");
    surf->fLogY = (logy && logy->GetString() == "yes");

    delete dir;

    return surf;
  }

  //----------------------------------------------------------------------
  std::unique_ptr<FrequentistSurface> FrequentistSurface::LoadFromMulti(
    std::vector<TFile*> files, std::string label)
  {
    std::vector<std::unique_ptr<FrequentistSurface>> surfs;
    for (TFile* f : files) {
      surfs.push_back(FrequentistSurface::LoadFrom(f, label));
    }

    int Nx = surfs[0]->fHist->GetNbinsX();
    int Ny = surfs[0]->fHist->GetNbinsY();
    size_t nbins = Nx * Ny;
    std::vector<int> binMask;

    // Loop over the surfaces to check all bins are present
    for (size_t i = 0; i < surfs.size(); ++i) {
      for (int bin : surfs[i]->fBinMask) {
        if (std::count(binMask.begin(), binMask.end(), bin))
          assert(false && "Repeated bin found in surfaces being merged.");
        binMask.push_back(bin);
      }
    }
    if (binMask.size() != nbins) {
      std::cout << "Missing bins found in surfaces being merged. "
        << "Are you sure you included all files for this surface?"
        << std::endl;
      assert(false);
    }

    DontAddDirectory guard;

    // Create return surface and initialise with first in list
    std::unique_ptr<FrequentistSurface> ret(new FrequentistSurface);
    ret->fHist = (TH2F*)surfs[0]->ToTH2(0)->Clone();
    for (TH2* h : surfs[0]->fProfHists)
      ret->fProfHists.push_back((TH2F*)h->Clone());

    ret->fLogX = surfs[0]->fLogX;
    ret->fLogY = surfs[0]->fLogY;
    ret->fSeedValues = surfs[0]->fSeedValues;

    // Find job that contains minimum bin (and therefore true best fit)
    int best = -1;
    double minChi = std::numeric_limits<double>::max();
    for (size_t i = 0; i < surfs.size(); ++i) {
      TH2F* hist = (TH2F*)surfs[i]->ToTH2(0);
      for (size_t bin : surfs[i]->fBinMask) {
        const int x = bin % Nx + 1;
        const int y = bin / Nx + 1;
        double val = hist->GetBinContent(x, y);
        if (val < minChi) {
          minChi = val;
          best = i;
        }
      }
    }

    ret->fBestLikelihood = surfs[best]->fBestLikelihood;
    ret->fBestFitX = surfs[best]->fBestFitX;
    ret->fBestFitY = surfs[best]->fBestFitY;

    // Loop over other surfaces and add them in
    for (size_t i = 1; i < surfs.size(); ++i) {
      ret->fHist->Add(surfs[i]->ToTH2(0));
      for (size_t j = 0; j < ret->fProfHists.size(); ++j)
        ret->fProfHists[j]->Add(surfs[i]->fProfHists[j]);
    }

    // Scale hist by global minimum
    for(int x = 0; x < ret->fHist->GetNbinsX()+2; ++x){
      for(int y = 0; y < ret->fHist->GetNbinsY()+2; ++y){
        ret->fHist->SetBinContent(x, y, ret->fHist->GetBinContent(x, y)-ret->fBestLikelihood);
      }
    }

    ret->fHist->SetMinimum(0);

    return ret;
  }

  //----------------------------------------------------------------------
  std::unique_ptr<FrequentistSurface> FrequentistSurface::LoadFromMulti(
    std::string wildcard, std::string label)
  {
    std::vector<std::string> fileNames = Wildcard(wildcard);
    std::vector<TFile*> files;
    for (std::string fileName : fileNames) {
      files.push_back(TFile::Open(fileName.c_str(), "read"));
    }
    std::unique_ptr<FrequentistSurface> ret = LoadFromMulti(files, label);
    for(TFile* f: files) delete f;
    return ret;
  }

  // See eg the statistics section of the PDG
  TH2* Gaussian68Percent2D(const FrequentistSurface& s){return Flat(2.30, s);}
  TH2* Gaussian90Percent2D(const FrequentistSurface& s){return Flat(4.61, s);}
  TH2* Gaussian95Percent2D(const FrequentistSurface& s){return Flat(5.99, s);}
  TH2* Gaussian2Sigma2D   (const FrequentistSurface& s){return Flat(6.18, s);}
  TH2* Gaussian99Percent2D(const FrequentistSurface& s){return Flat(9.21, s);}
  TH2* Gaussian3Sigma2D   (const FrequentistSurface& s){return Flat(11.83, s);}

  TH2* Gaussian68Percent1D(const FrequentistSurface& s){return Flat(1.00, s);}
  TH2* Gaussian90Percent1D(const FrequentistSurface& s){return Flat(2.71, s);}
  TH2* Gaussian95Percent1D(const FrequentistSurface& s){return Flat(3.84, s);}
  TH2* Gaussian2Sigma1D   (const FrequentistSurface& s){return Flat(4.00, s);}
  TH2* Gaussian99Percent1D(const FrequentistSurface& s){return Flat(6.63, s);}
  TH2* Gaussian3Sigma1D   (const FrequentistSurface& s){return Flat(9.00, s);}


} // namespace