NCSpectrumInterpolator.cxx

Go to the documentation of this file.

//$Id: NCSpectrumInterpolator.cxx,v 1.18 2009/11/18 18:58:24 bckhouse Exp $
//
// NCSpectrumInterpolator.cxx
//
// Given a set of reference histograms, interpolate between them.
//
// C. Backhouse 10/2008

#include "NCSpectrumInterpolator.h"

#include "MessageService/MsgService.h"

#include "NCUtils/NCUtility.h"
using namespace NC::Utility;

#include "TCanvas.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TMath.h"
#include "TMatrixD.h"

#include <cassert>
#include <cmath>

CVSID("$Id: NCSpectrumInterpolator.cxx,v 1.18 2009/11/18 18:58:24 bckhouse Exp $");

using std::vector;

bool NC::ISpectrumInterpolator::ShiftInfo::IsNominal() const
{
  for(unsigned int n = 0; n < shift.size(); ++n)
    if(shift[n] != 0) return false;
  return true;
}


NC::ISpectrumInterpolator::~ISpectrumInterpolator()
{
  for(unsigned int n = 0; n < shifted.size(); ++n)
    for(unsigned int m = 0; m < shifted[n].spectra.size(); ++m)
      delete shifted[n].spectra[m];
}


void NC::ISpectrumInterpolator::AddInputSpectra(std::vector<double> shift,
                                                std::vector<TH1*> spectra)
{
  for(unsigned int n = 0; n < shifted.size(); ++n){
    if(shifted[n].shift == shift){
      assert(0 && "We already added spectra at this shift");
      //      shifted[n].spectra.insert(shifted[n].spectra.end(),
      //                                spectra.begin(), spectra.end());
      //      return;
    }
  }

  // Didn't find a matching one so make a new one
  shifted.push_back(ShiftInfo(shift, spectra));
}


vector<TH1*> NC::ISpectrumInterpolator::
GetInputSpectra(vector<double> shift) const
{
  vector<TH1*> ret;

  for(unsigned int n = 0; n < shifted.size(); ++n){
    const ShiftInfo& sn = shifted[n];
    if(sn.shift == shift){
      ret.reserve(sn.spectra.size());
      for(unsigned int n = 0; n < sn.spectra.size(); ++n)
        ret.push_back(CloneFast(sn.spectra[n]));
      return ret;
    }
  }

  return ret; // empty
}


vector<TH1*> NC::ISpectrumInterpolator::
InterpolateHists(const vector<ShiftInfo>& hs,
                 const vector<double>& shift) const
{
  assert(!hs.empty());

  vector<vector<double> > xs;
  xs.reserve(hs.size());
  for(unsigned int n = 0; n < hs.size(); ++n) xs.push_back(hs[n].shift);

  const TMatrixD coeffs = FindCoeffs(xs, shift);

  const unsigned int I = hs[0].spectra.size();
  vector<TH1*> ret;
  ret.reserve(I);
  for(unsigned int i = 0; i < I; ++i){
    // Get the correct binning
    TH1* h = CloneFast(hs[0].spectra[i]);
    h->Reset();

    const unsigned int J = hs.size();

    for(unsigned int j = 0; j < J; ++j)
      AddFast(h, hs[j].spectra[i], coeffs(0, j));

    ret.push_back(h);
  } // end for i

  return ret;
}


int NC::SpectrumInterpolatorFancy::intintpow(int a, int b) const
{
  int ret = 1;
  while(b--) ret *= a;
  return ret;
}


double NC::SpectrumInterpolatorFancy::doubleintpow(double a, int b) const
{
  double ret = 1;
  while(b--) ret *= a;
  return ret;
}

/*
double NC::SpectrumInterpolator::
InterpolateMultiPoints(const vector<vector<double> >& xs,
                       const vector<double>& ys,
                       const vector<double>& x) const
{
  assert(xs.size() == ys.size());

  assert(xs[0].size() == x.size());

  // 2 to fit a general quadratic, 1 to fit a general linear function, etc.
  const int fitOrder = 2;

  // The dimensionality of the space we are interpolating in
  const int D = x.size();

  // The size of the matrices, which is the number of coefficients we
  // need to determine, which is the number of different products of the
  // components of the x coordinate each raised to a power up to fitOrder
  const int K = intintpow(fitOrder+1, D);

  assert(xs.size() >= (unsigned int)K);

  TMatrixD mat(K, K);
  TMatrixD vec(K, 1);

  // What powers to raise the xs to for each row in the matrix
  vector<vector<int> > pows(K);
  for(int n = 0; n < K; ++n) pows[n] = vector<int>(xs.size());

  // Calculate every combination of powers
  // This produces a matrix like |012012012012012012012012012|
  // (where k is the col number, |000111222000111222000111222|
  // and d the row number)       |000000000111111111222222222|
  for(int k = 0; k < K; ++k)
    for(int d = 0; d < D; ++d)
      pows[k][d] = (k%intintpow(fitOrder+1, d+1))/intintpow(fitOrder+1, d);

  // The number of data points we are fitting
  const int N = xs.size();

  // How close is really close?
  const double eps = 1e-12;
  // If evaluating really close to a known point, just use that point
  // since when  x-x_n -> 0 we get a divide-by-zero
  for(int n = 0; n < N; ++n){
    double dist = 0;
    for(int d = 0; d < D; ++d){
      dist += SQR(xs[n][d]-x[d]);
    }
    if(dist < eps) return ys[n];
  }

  // The equation we are solving is mat*a = vec

  // Loop over all the data points
  for(int n = 0; n < N; ++n){
    // The weight by which this point is discounted in the fit
    // denom=1 corresponds to the usual unweighted fit
    // TODO - only taking weight in one dimension
    double denom = 0;
    for(int d = 0; d < D; ++d) denom += SQR(x[d]-xs[n][d]);

    // Loop over the rows in the vector
    for(int j = 0; j < K; ++j){
      double vecadd = ys[n]/denom;
      for(int d = 0; d < D; ++d)
        vecadd *= doubleintpow(xs[n][d], pows[j][d]);
      vec(j, 0) += vecadd;
    }

    // Loop over the cells of the matrix
    for(int j = 0; j < K; ++j){
      for(int i = 0; i < K; ++i){
        double matadd = 1/denom;
        for(int d = 0; d < D; ++d){
          matadd *= doubleintpow(xs[n][d], pows[j][d]);
          matadd *= doubleintpow(xs[n][d], pows[i][d]);
        }
        mat(j, i) += matadd;
      }
    }
  }

  // This part solves the matrix equation
  mat.Invert();

  const TMatrixD coeffs = mat*vec;

  // And now calculate the value of the polynomial at x
  double ret = 0;
  for(int k = 0; k < K; ++k){
    double retadd = coeffs(k, 0);
    for(int d = 0; d < D; ++d){
      retadd *= doubleintpow(x[d], pows[k][d]);
    }
    ret += retadd;
  }
  return ret;
}
*/


void NC::ISpectrumInterpolator::WriteResources() const
{
  MSG("NCISpectrumInterpolator",
      Msg::kInfo) << "NC::ISpectrumInterpolator::WriteResources()" << endl;

  const int N = int(shifted.size());
  if(N == 0) return;
  const int M = int(shifted[0].spectra.size());

  // Go through all the plots, find the maximum range and extend by 10%
  // then clamp to the range 0.5-2.0
  double min = FLT_MAX;
  double max = FLT_MIN;
  for(int n = 0; n < N; ++n){
    for(int m = 0; m < M; ++m){
      TH1* s = shifted[n].spectra[m];
      if(s->GetMinimum(0) < min) min = s->GetMinimum(0);
      if(s->GetMaximum() > max) max = s->GetMaximum();
    }
  }
  max = 1+(max-1)*1.1;
  min = 1-(1-min)*1.1;
  if(max > 2) max = 2;
  if(min < 0.5) min = 0.5;

  for(int n = 0; n < N; ++n){
    for(int m = 0; m < M; ++m){
      TH1* s = shifted[n].spectra[m];
      s->SetMinimum(min);
      s->SetMaximum(max);
      s->Write();
    }
  }

  TCanvas canv("all_spectra_canv");
  // Try dividing as a square. If this isn't big enough then expand
  // horizontally. If still not then also expand vertically.
  int X = int(TMath::Sqrt(N*M));
  int Y = int(TMath::Sqrt(N*M));
  if(X*Y < N*M) ++X;
  if(X*Y < N*M) ++Y;
  canv.Divide(X, Y);

  for(int n = 0; n < N; ++n){
    for(int m = 0; m < M; ++m){
      canv.cd(M*n+m+1);
      TH1* s = shifted[n].spectra[m];
      s->SetTitle(s->GetName());
      if(s->GetXaxis()->GetXmin() == 0 && s->GetXaxis()->GetXmax() == 120)
        s->GetXaxis()->SetRangeUser(0, 19.5);
      if(s->GetDimension() == 2) s->Draw("colz"); else s->Draw();
    }
  }
  canv.Write();

  MSG("NCISpectrumInterpolator",
      Msg::kInfo) << "Done writing interpolator resources" << endl;
}


vector<TH1*> NC::ISpectrumInterpolator::GetNominal() const
{
  for(unsigned int n = 0; n < shifted.size(); ++n){
    if(shifted[n].IsNominal()) return shifted[n].spectra;
  }
  assert(0 && "No nominal spectrum exists");
}


ostream& operator<<(ostream& op, const TMatrixD& m)
{
  op << "| ";
  for(int n = 0; n < m.GetNcols(); ++n) op << m(0, n) << " ";
  op << "|";
  return op;
}

ostream& operator<<(ostream& op, const vector<double>& m)
{
  op << "[ ";
  for(unsigned int n = 0; n < m.size(); ++n) op << m[n] << " ";
  op << "]";
  return op;
}


TMatrixD NC::SpectrumInterpolatorFancy::
FindCoeffs(const vector<vector<double> >& xs, const vector<double>& x) const
{
  assert(xs[0].size() == x.size());

  // 2 to fit a general quadratic, 1 to fit a general linear function, etc.
  const int fitOrder = 1;

  // The dimensionality of the space we are interpolating in
  const int D = x.size();

  // The size of the matrices, which is the number of coefficients we
  // need to determine, which is the number of different products of the
  // components of the x coordinate each raised to a power up to fitOrder
  const int K = intintpow(fitOrder+1, D);

  assert(xs.size() >= (unsigned int)K);


  // What powers to raise the xs to for each row in the matrix
  vector<vector<int> > pows(K);
  for(int k = 0; k < K; ++k) pows[k] = vector<int>(xs.size());

  // Calculate every combination of powers
  // This produces a matrix like |012012012012012012012012012|
  // (where k is the col number, |000111222000111222000111222|
  // and d the row number)       |000000000111111111222222222|
  for(int k = 0; k < K; ++k)
    for(int d = 0; d < D; ++d)
      pows[k][d] = (k%intintpow(fitOrder+1, d+1))/intintpow(fitOrder+1, d);


  // Calculate the value of the different polynomial parts at x
  TMatrixD basis(1, K);
  for(int k = 0; k < K; ++k){
    double basisAdd = 1;
    for(int d = 0; d < D; ++d){
      basisAdd *= doubleintpow(x[d], pows[k][d]);
    }
    basis(0, k) = basisAdd;
  } // end for k


  // The number of data points we are fitting
  const int N = xs.size();

  // How close is really close?
  const double eps = 1e-12;
  // If evaluating really close to a known point, just use that point
  // since when  x-x_n -> 0 we get a divide-by-zero
  for(int n = 0; n < N; ++n){
    double dist = 0;
    for(int d = 0; d < D; ++d){
      dist += SQR(xs[n][d]-x[d]);
    }
    if(dist < eps){
      TMatrixD coeffs(1, N);
      coeffs(0, n) = 1;
      return coeffs;
    }
  } // end for n

  TMatrixD mat(K, K);

  TMatrixD vecPart(K, N);

  // Loop over all the data points
  for(int n = 0; n < N; ++n){
    // The weight by which this point is discounted in the fit
    // denom=1 corresponds to the usual unweighted fit
    // TODO - only taking weight in one dimension
    double denom = 0;
    for(int d = 0; d < D; ++d) denom += SQR(x[d]-xs[n][d]);

    // Loop over the rows in the vector
    for(int j = 0; j < K; ++j){
      double vecadd = 1/denom;
      for(int d = 0; d < D; ++d)
        vecadd *= doubleintpow(xs[n][d], pows[j][d]);
      vecPart(j, n) = vecadd;
    }

    // Loop over the cells of the matrix
    for(int j = 0; j < K; ++j){
      for(int i = 0; i < K; ++i){
        double matadd = 1/denom;
        for(int d = 0; d < D; ++d){
          matadd *= doubleintpow(xs[n][d], pows[j][d]);
          matadd *= doubleintpow(xs[n][d], pows[i][d]);
        }
        mat(j, i) += matadd;
      } // end for i
    } // end for j
  } // end for n
  mat.Invert();

  return basis*mat*vecPart;
}


TMatrixD NC::SpectrumInterpolatorSimple::
FindCoeffs(const vector<vector<double> >& xs, const vector<double>& x) const
{
  //  const unsigned int D = x.size(); // dimensionality of the space

  // Ideally these should match, but the case where the reference
  // points are in more dimensions than the input should still work.
  assert(xs[0].size() <= x.size());
  const unsigned int D = TMath::Min(x.size(), xs[0].size());

  const unsigned int N = xs.size(); // number of input points

  TMatrixD ret(1, N);

  for(unsigned int n = 0; n < N; ++n){ // for every point
    double weight = 1;
    for(unsigned int d = 0; d < D; ++d){ // for every dimension
      const bool ltOneSigma = x[d] > -1 && x[d] < +1;

      if(ltOneSigma){
        const double dist = fabs(xs[n][d]-x[d]);
        if(dist > 1)
          weight = 0;
        else
          weight *= (1-dist);
      }
      else{
        const bool sign = xs[n][d]/x[d] >= 0; // right sign

        if(sign){
          if(fabs(xs[n][d]) < 1)
            weight *= 1-fabs(x[d]);
          else
            weight *= fabs(x[d]);
        }
        else weight = 0;
      }
    } // end for d
    ret(0, n) = weight;
  } // end for n

  //  cout << x << " -> " << ret << endl;

  double norm = 0;
  for(unsigned int n = 0; n < N; ++n) norm += ret(0, n);
  assert(fabs(norm - 1) < 1e-3);


  return ret;
}


// Need every point to be shifted in one dimension only. Give weight from
// each point linear in how far 'x' is shifted. We are dealing with ratios
// not shifts here - so need to subtract off a sufficient amount of the nominal
// spectrum.
TMatrixD NC::SpectrumInterpolatorSimplest::
FindCoeffs(const vector<vector<double> >& xs, const vector<double>& x) const
{
  const unsigned int N = xs.size(); // number of input points
  TMatrixD ret(1, N);

  const unsigned int D = TMath::Min(x.size(), xs[0].size());

  // If there is more than one point, and no dimensions, then they'll all be
  // picked as nominal. Go find out why your systematic shift is empty.
  // (perhaps because we don't support fitting the systematic you're trying
  // to use).
  assert((N == 1 || D > 0) && "Empty shift");

  int nominalIdx = -1;

  for(unsigned int n = 0; n < N; ++n){ // for every point
    bool isNominal = true;
    double weight = 0;
    for(unsigned int d = 0; d < D; ++d){ // for every dimension
      if(xs[n][d] != 0) isNominal = false;
      if((xs[n][d] < 0 && x[d] < 0) ||
         (xs[n][d] > 0 && x[d] > 0)){
        assert(weight == 0); // No point should contribute more than once
        weight = x[d] / xs[n][d];
      } // end if
    } // end for d
    ret(0, n) = weight;
    if(isNominal){
      // Should only be one nominal point
      assert(nominalIdx == -1 && "Multiple nominal spectra");
      nominalIdx = n;
    }
  } // end for n
  assert(nominalIdx >= 0);

  assert(ret(0, nominalIdx) == 0); // Nominal shouldn't have got weight
  double norm = 0;
  for(unsigned int n = 0; n < N; ++n) norm += ret(0, n);
  ret(0, nominalIdx) = 1-norm; // Make the total up to one

  //  cout << x << " -> " << ret << endl;

  return ret;
}

---