testErrorPropagator.cc

Go to the documentation of this file.

#include <tst/Verify.h>

#include <utl/AugerUnits.h>
#include <utl/AugerException.h>

#include <iostream>
#include <cmath>
#include <vector>

#include <cppunit/extensions/HelperMacros.h>
#include <boost/tuple/tuple.hpp>

#include <utl/NumericalErrorPropagator.h>
#include <utl/CovarianceMatrix.h>
#include <utl/AugerException.h>

using namespace std;
using namespace utl;
using namespace tst;

#define ASSERT_EQUAL(x, y) CPPUNIT_ASSERT(x == y)

class testErrorPropagator : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(testErrorPropagator);
  CPPUNIT_TEST(TestCovarianceMatrix);
  CPPUNIT_TEST(TestErrorPropagatorFaster);
  CPPUNIT_TEST(TestErrorPropagatorBetter);
  CPPUNIT_TEST_SUITE_END();
  
private:
  
  template<class Stream>
  void
  OutputMatrix(Stream& s, const CovarianceMatrix& c)
  {
    for (unsigned int i=0;i<c.GetExtent();++i)
      for (unsigned int j=0;j<c.GetExtent();++j)
        s << c(i,j) << (j<c.GetExtent()-1 ? " " : "\n");
  }

  // output 3, input 4
  class TestFunctor : public NumericalErrorPropagator
  {
  public:
    double m[3][4];

    void
    Transform(std::vector<double>& out, const std::vector<double>& inp)
      const
    {
      if (out.size()!=3) out = std::vector<double>(3);
      out[0] = m[0][0]*inp[0]+m[0][1]*inp[1]+m[0][2]*inp[2]+m[0][3]*inp[3] + 10.;
      out[1] = m[1][0]*inp[0]+m[1][1]*inp[1]+m[1][2]*inp[2]+m[1][3]*inp[3] + 20.;
      out[2] = m[2][0]*inp[0]+m[2][1]*inp[1]+m[2][2]*inp[2]+m[2][3]*inp[3] + 30.;
    }

  };

public:

  void
  TestCovarianceMatrix()
  {
    { // test with nPar = 2
      CovarianceMatrix c(2);
      c(0,0) = 0.1;
      c(0,1) = 0.2;
      c(1,1) = 0.3;

      ostringstream s;
      OutputMatrix(s,c);
      
      const string truth =
        "0.1 0.2\n"
        "0.2 0.3\n";

      ASSERT_EQUAL(s.str(),truth);
    }
    
    { // test with nPar = 6
      CovarianceMatrix c(6);
      c[0] = 0;
      c[1] = 1;
      c[2] = 2;
      c[3] = 3;
      c[4] = 4;
      c[5] = 5;
      c(0,1) = 0.1;
      c(0,2) = 0.2;
      c(0,3) = 0.3;
      c(0,4) = 0.4;
      c(0,5) = 0.5;
      c(1,2) = 1.2;
      c(1,3) = 1.3;
      c(1,4) = 1.4;
      c(1,5) = 1.5;
      c(2,3) = 2.3;
      c(2,4) = 2.4;
      c(2,5) = 2.5;
      c(3,4) = 3.4;
      c(3,5) = 3.5;
      c(4,5) = 4.5;

      ostringstream s;
      OutputMatrix(s,c);
      
      const string truth =
        "0 0.1 0.2 0.3 0.4 0.5\n"
        "0.1 1 1.2 1.3 1.4 1.5\n"
        "0.2 1.2 2 2.3 2.4 2.5\n"
        "0.3 1.3 2.3 3 3.4 3.5\n"
        "0.4 1.4 2.4 3.4 4 4.5\n"
        "0.5 1.5 2.5 3.5 4.5 5\n";

      ASSERT_EQUAL(s.str(),truth);
    }
  }

  void
  TestErrorPropagatorFaster()
  {
    TestErrorPropagator(NumericalErrorPropagator::eFaster);
  }

  void
  TestErrorPropagatorBetter()
  {
    TestErrorPropagator(NumericalErrorPropagator::eBetter);
  }

  void
  TestErrorPropagator(NumericalErrorPropagator::Quality quality)
  {
    vector<double> iVec(4);
    CovarianceMatrix iCov(4);

    iVec[0] = 1.23;
    iVec[1] = 4.56;
    iVec[2] = 7.89;
    iVec[3] = 2.01;

    iCov[0] = 0.321;
    iCov[1] = 0.654;
    iCov[2] = 0.987;
    iCov[3] = 0.234;
    iCov(0,1) = 0.11;
    iCov(0,2) = -0.12;
    iCov(0,3) = -0.13;
    iCov(1,2) = 0.22;
    iCov(1,3) = -0.23;
    iCov(2,3) = -0.34;
    
    TestFunctor fcn; // does y = m*x + b
    fcn.m[0][0] =  1.1;
    fcn.m[0][1] = -2.2;
    fcn.m[0][2] =  3.3;
    fcn.m[0][3] =  4.4;
    fcn.m[1][0] = -4.5;
    fcn.m[1][1] = -5.6;
    fcn.m[1][2] = -6.7;
    fcn.m[1][3] =  7.8;
    fcn.m[2][0] =  0.1;
    fcn.m[2][1] =  0.2;
    fcn.m[2][2] = -0.3;
    fcn.m[2][3] =  0.5;

    vector<double> oVec(3);
    CovarianceMatrix oCov(3);
    fcn(oVec,oCov,iVec,iCov,quality);

    // V_out = J V_in J^T, mit J = Jacobi matrix, here for TestFunctor J = m

    double tmp[3][4];
    for (int i=0;i<3;++i)
      for (int j=0;j<4;++j)
    {
      tmp[i][j] = 0.0;
      for (int k=0;k<4;++k)
        tmp[i][j] += fcn.m[i][k]*iCov(k,j);
    }

    ostringstream s1;
    for (int i=0;i<3;++i)
      for (int j=0;j<3;++j)
    {
      double vout = 0.0;
      for (int k=0;k<4;++k)
        vout += tmp[i][k]*fcn.m[j][k]; // fcn.m[j][k] is transposed of fcn.m[k][j]
      s1 << vout << (j<2 ? " " : "\n");
    }
    
    ostringstream s2;
    for (int i=0;i<3;++i)
      for (int j=0;j<3;++j)
    {
      s2 << oCov(i,j) << (j<2 ? " " : "\n");
    }

    ASSERT_EQUAL(s1.str(),s2.str());
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(testErrorPropagator);