NumericalErrorPropagation.cc

Go to the documentation of this file.

#include <iostream>
#include <vector>

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/io.hpp>

#include <utl/ErrorLogger.h>
#include <utl/AugerException.h>
#include <utl/NumericalErrorPropagation.h>
#include <utl/CorrelationMatrix.h>

using namespace std;
using namespace boost::numeric;
using utl::CorrelationMatrix;
using utl::NumericalErrorPropagation;
using utl::FCNCalculator;
using utl::Parameter;


NumericalErrorPropagation::NumericalErrorPropagation(const FCNCalculator& fcn,
                                                     const vector<Parameter>& par,
                                                     const CorrelationMatrix& corr) :
  fFcn(fcn),
  fParameter(par),
  fCorrelation(corr)
{
  if (par.size() != corr.GetNParameter()) {
    ostringstream err;
    err << "ERROR "
           "#parameter: " << par.size() << " is different from size of "
           "correlation matrix " << corr.GetNParameter();
    ERROR(err);
    throw OutOfBoundException(err.str());
  }

  Propagate();
}


NumericalErrorPropagation::~NumericalErrorPropagation()
{
  delete fPropagatedCorrelation;
}


void
NumericalErrorPropagation::Propagate()
{
  /*
    calculate Jacobi-Matrix

    alpha_j: input parameter j
    q_i (alpha): output-parameter i for set of input parameters
    q_i (alpha_1, alpha_2, ... , alpha_n)

    J_ij = partial ( d(q_i) / d(alpha_j) )
  */

  const unsigned int J_cols = fParameter.size(); // alpha_j input-parameter

  vector<double> x(J_cols),xm(J_cols),xmm(J_cols),xp(J_cols),xpp(J_cols);

  for (unsigned int i = 0; i < J_cols; ++i)  // initialize
    x[i] = xm[i] = xmm[i] = xp[i] = xpp[i] = fParameter[i].value;

  fResult = fFcn(x);

  const unsigned int J_rows = fResult.size();  // q_i output-parameter

  ublas::matrix<double> J(J_rows, J_cols);

  for (unsigned col = 0; col < J_cols; ++col) {  // alpha_j

    if (col > 0) {
      // reset previous change to vectors, x is untouched
      xpp[col-1] = x[col-1];
      xp [col-1] = x[col-1];
      xm [col-1] = x[col-1];
      xmm[col-1] = x[col-1];
    }

    /* if error of alpha_j is zero, we can skip calculation of all
       d q_i / d alpha_j, they are multiplied with zero later -
       still guard against uninitialized values in J, though */
    if (fParameter[col].error) {
      /* try to guess dx which is small compared to curvature of fFcn, but large compared
         to machine precision, since we don't want to see numerical round off errors */
      double dx = fParameter[col].error * 1e-3;

      // derivative numerically calculated with five-point stencil (see e.g. Wikipedia)
      xpp[col] += 2*dx;
      xp[col] += dx;
      xm[col] -= dx;
      xmm[col] -= 2*dx;

      vector<double> fpp = fFcn(xpp);
      vector<double> fp = fFcn(xp);
      vector<double> fm = fFcn(xm);
      vector<double> fmm = fFcn(xmm);

      for (unsigned row = 0; row < J_rows; ++row)  // q_i
        J(row, col) = (-fpp[row] + 8*fp[row] - 8*fm[row] + fmm[row]) / (12*dx);

    } else {

      for (unsigned row = 0; row < J_rows; ++row)  // q_i
        J(row, col) = 0;

    }
  }

  /*
    calculate the covariance function of the output parameters q_i
    V' = J * V * J^T
  */

  ublas::matrix<double> V_input(J_cols, J_cols);

  for (unsigned row = 0; row < J_cols; ++row)
    for (unsigned col = 0; col < J_cols; ++col)
      V_input(row,col) =
        fCorrelation.Get(row, col) * fParameter[row].error * fParameter[col].error;

  const ublas::matrix<double> temp = ublas::prod(V_input, trans(J));
  ublas::matrix<double> V_output = ublas::prod(J, temp);

  /*
    now retrieve the output parameter errors and their correlations
  */

  fPropagatedErrors.clear();
  fPropagatedErrors.reserve(J_rows);
  for (unsigned i = 0; i < J_rows; ++i) {
    const double V = V_output(i, i);

    if (V < 0 || std::isnan(V)) {
      ostringstream err;
      err << "invalid variance " << V;
      ERROR(err);
      fPropagatedErrors.push_back(0);
      V_output(i, i) = 0;
    } else
      fPropagatedErrors.push_back(sqrt(V));
  }

  fPropagatedCorrelation = new CorrelationMatrix(J_rows);

  for (unsigned i = 0; i < J_rows; ++i) {
    for (unsigned j = i+1; j < J_rows; ++j) {

      const double tmp = V_output(i, i) * V_output(j, j);

      const double correlation = tmp > 0 ? V_output(i, j) / sqrt(tmp) : 0;

      if (correlation < -1 || correlation > 1 || std::isnan(correlation)) {
        ostringstream err;
        err << "correlation = " << correlation << " is invalid, converted to 0 (!)";
        ERROR(err);
        fPropagatedCorrelation->Set(i, j, 0);
      } else
        fPropagatedCorrelation->Set(i, j, correlation);
    }
  }
}