SurfaceKrige.cxx

Go to the documentation of this file.

#include "CAFAna/Fit/SurfaceKrige.h"

#include "CAFAna/Core/Kriger.h"
#include "CAFAna/Core/IFitVar.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/Utilities.h"

#include "Utilities/func/MathUtil.h"

#include "TH2.h"

namespace ana
{
  //----------------------------------------------------------------------
  KrigeKernel::KrigeKernel()
    : fGrad(1), fGradX(1), fGradY(1)
  {
  }

  //----------------------------------------------------------------------
  double KrigeKernel::operator()(double x1, double y1, double x2, double y2) const
  {
    // TODO - the 2D model doesn't seem to actually work
    if(true){//fPts.size() < 10){
      // Use a very simple function when we don't have much data
      return fGrad*util::pythag(x1-x2, y1-y2);
    }
    else{
      // Treat x and y separately later
      return fGradX*fabs(x1-x2) + fGradY*fabs(y1-y2);
    }
  }

  //----------------------------------------------------------------------
  void KrigeKernel::AddPoint(double x, double y, double z)
  {
    // Learn the best function of dx and dy to predict the variance

    Pt p = {x, y, z};
    fPts.push_back(p);

    double R2 = 0, ZR = 0;
    double X2 = 0, XY = 0, Y2 = 0, ZX = 0, ZY = 0;
    for(unsigned int ai = 0; ai < fPts.size(); ++ai){
      for(unsigned int bi = ai+1; bi < fPts.size(); ++bi){
        const Pt a = fPts[ai];
        const Pt b = fPts[bi];

        const double dz2 = util::sqr(a.z-b.z);

        const double dx = fabs(a.x-b.x);
        const double dy = fabs(a.y-b.y);

        const double dr = util::pythag(a.x-b.x, a.y-b.y);

        R2 += dr*dr;
        ZR += dz2*dr;

        X2 += dx*dx;
        XY += dx*dy;
        Y2 += dy*dy;
        ZX += dz2 * dx;
        ZY += dz2 * dy;
      }
    }

    if(R2) fGrad = ZR / R2;

    // Only compute 2D coefficients once we have some stats
    if(fPts.size() < 10) return;

    TMatrixD M(2, 2);              TVectorD v(2);
    M(0, 0) = X2; M(0, 1) = XY;    v(0) = ZX;
    M(1, 0) = XY; M(1, 1) = Y2;    v(1) = ZY;

    M.Invert();
    TVectorD g = M*v;

    fGradX = fabs(g(0));
    fGradY = fabs(g(1));
  }

  //----------------------------------------------------------------------
  SurfaceKrige::
  SurfaceKrige(const IExperiment* expt,
               osc::IOscCalcAdjustable* calc,
               const IFitVar* xvar, int nbinsx, double xmin, double xmax,
               const IFitVar* yvar, int nbinsy, double ymin, double ymax,
               const std::vector<const IFitVar*>& profVars,
               const std::vector<const ISyst*>& profSysts,
               const SeedList& seedPts,
               const std::vector<SystShifts>& systSeedPts,
               MinuitFitter::FitOpts opts)
  {
    fLogX = fLogY = false;

    fParallel = false;
    fFitOpts = opts;

    CreateHistograms(FitAxis(xvar, nbinsx, xmin, xmax, fLogX),
                     FitAxis(yvar, nbinsy, ymin, ymax, fLogY),
                     profVars, profSysts);

    fErr = ExpandedHistogram(";"+xvar->LatexName()+";"+yvar->LatexName(),
                             nbinsx, xmin, xmax, fLogX,
                             nbinsy, ymin, ymax, fLogY);

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

    FillSurface(expt, calc, xvar, yvar, profVars, profSysts, seedPts, systSeedPts);

    FindMinimum(expt, calc, xvar, yvar, profVars, profSysts, seedPts, systSeedPts);
  }

  //---------------------------------------------------------------------
  SurfaceKrige::~SurfaceKrige()
  {
    delete fErr;
  }

  //---------------------------------------------------------------------
  void SurfaceKrige::
  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)
  {
    KrigeKernel kern;
    // Need this layer of indirection to capture the kernel by reference
    Kriger kriger([&kern](double x1, double y1, double x2, double y2)
                  {
                    return kern(x1, y1, x2, y2);
                  });

    // Start off with one point in the centre
    int targetx = fHist->GetNbinsX()/2;
    int targety = fHist->GetNbinsY()/2;

    Progress prog(ProgressBarTitle(xvar, yvar, profVars, profSysts));

    // TODO: Arbitrary number of passes, need real stopping condition...
    const int maxPass = 100;//fHist->GetNbinsX()*fHist->GetNbinsY()/10;

    for(int pass = 0; pass < maxPass; ++pass){

      // Evaluate the function at the target point and add the info into the
      // Kriger
      const double targetz = FillSurfacePoint(expt, calc,
                                              xvar, BinCenterX(targetx),
                                              yvar, BinCenterY(targety),
                                              profVars, profSysts,
                                              seedPts, systSeedPts);

      kriger.AddPoint(targetx, targety, targetz);

      kern.AddPoint(targetx, targety, targetz);

      prog.SetProgress(double(pass)/maxPass);

      // Search through the surface for the next most uncertain point
      double maxerr = -1;

      for(int ix = 1; ix <= fHist->GetNbinsX(); ++ix){
        for(int iy = 1; iy <= fHist->GetNbinsY(); ++iy){

          double err;
          const double z = kriger.Predict(ix, iy, &err);

          // Keep surface and error surface up-to-date as we go
          fHist->SetBinContent(ix, iy, z);
          fErr->SetBinContent(ix, iy, err);

          // New candidate
          if(err > maxerr){
            maxerr = err;
            targetx = ix;
            targety = iy;
          }
        } // end for iy
      } // end for ix

    } // end for pass
  }
}