Kriger.cxx

Go to the documentation of this file.

#include "CAFAna/Core/Kriger.h"

namespace ana
{
  //----------------------------------------------------------------------
  void Kriger::AddPoint(double x, double y, double z)
  {
    fPts.emplace_back(x, y, z);
  }

  //----------------------------------------------------------------------
  double Kriger::Predict(double x, double y, double* err)
  {
    double junk;
    if(!err) err = &junk;

    EnsureDecomp();

    Pt p(x, y, 0);

    TVectorD b(N()+1);
    for(unsigned int i = 0; i < N(); ++i){
      b(i) = fKern(p.x, p.y, fPts[i].x, fPts[i].y);
    }
    b(N()) = 1;

    bool ok;
    TVectorD lambda = fDecomp.Solve(b, ok);

    double z = 0;
    *err = 0;
    for(unsigned int i = 0; i < N(); ++i){
      z += lambda[i] * fPts[i].z;
      *err += lambda[i] * b[i];
    }
    *err = sqrt(fabs(*err));

    return z;
  }

  //----------------------------------------------------------------------
  void Kriger::EnsureMatrix()
  {
    if(unsigned(fMat.GetNrows()) == N()+1) return;
    fMat.ResizeTo(N()+1, N()+1);

    for(unsigned int i = 0; i < N(); ++i){
      fMat(i, i) = 0;
      for(unsigned int j = i+1; j < N(); ++j){
        fMat(i, j) = fMat(j, i) = fKern(fPts[i].x, fPts[i].y, fPts[j].x, fPts[j].y);
      }
    }

    // Sum of coefficients is 1
    for(unsigned int i = 0; i < N(); ++i) fMat(N(), i) = 1;
    // Include the lagrange multiplier in each
    for(unsigned int i = 0; i < N(); ++i) fMat(i, N()) = 1;
  }

  //----------------------------------------------------------------------
  void Kriger::EnsureDecomp()
  {
    EnsureMatrix();

    if(fDecomp.GetNrows() == fMat.GetNrows()) return;
    fDecomp.SetMatrix(fMat);
  }
}