testErrorPropagation.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>

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

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

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

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

class testErrorPropagation : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(testErrorPropagation);
  CPPUNIT_TEST(TestCorrelationMatrix);
  CPPUNIT_TEST(TestErrorPropagation);
  CPPUNIT_TEST_SUITE_END();
  
private:
  
  template<class Stream>
  void
  OutputMatrix(Stream& s, const CorrelationMatrix& c)
  {
    for (unsigned int i=0;i<c.GetNParameter();++i)
    {
      for (unsigned int j=0;j<c.GetNParameter()-1;++j)
        s << c.Get(i,j) << " ";
      s << c.Get(i,c.GetNParameter()-1) << "\n";
    }
  }

  // input 4, output 3
  class TestFunctor : public FCNCalculator {
    public:

    double m[3][4];

    std::vector<double> operator()(const std::vector<double>& inp) const
    {
      std::vector<double> out;
      out.push_back(m[0][0]*inp[0]+m[0][1]*inp[1]+m[0][2]*inp[2]+m[0][3]*inp[3] + 10.);
      out.push_back(m[1][0]*inp[0]+m[1][1]*inp[1]+m[1][2]*inp[2]+m[1][3]*inp[3] + 20.);
      out.push_back(m[2][0]*inp[0]+m[2][1]*inp[1]+m[2][2]*inp[2]+m[2][3]*inp[3] + 30.);
      return out;
    }

  };

public:

  void
  TestCorrelationMatrix()
  {
    { // test with nPar = 2
      CorrelationMatrix c(2);
      c.Set(0,1,0.3);

      ostringstream s;
      OutputMatrix(s,c);
      
      const string truth =
        "1 0.3\n"
        "0.3 1\n";
      ASSERT_EQUAL(s.str(),truth);
    }
    
    { // test with nPar = 6
      CorrelationMatrix c(6);
      c.Set(0,1,0.1);
      c.Set(2,0,0.2);
      c.Set(0,3,0.3);
      c.Set(4,0,0.4);
      c.Set(0,5,0.5);
      c.Set(1,2,-0.11);
      c.Set(1,3,-0.21);
      c.Set(4,1,-0.31);
      c.Set(1,5,-0.41);
      c.Set(2,3,0.12);
      c.Set(2,4,0.22);
      c.Set(2,5,0.32);
      c.Set(4,3,-0.13);
      c.Set(3,5,-0.23);
      c.Set(5,4,0.111);

      ostringstream s;
      OutputMatrix(s,c);
      
      const string truth =
        "1 0.1 0.2 0.3 0.4 0.5\n"
        "0.1 1 -0.11 -0.21 -0.31 -0.41\n"
        "0.2 -0.11 1 0.12 0.22 0.32\n"
        "0.3 -0.21 0.12 1 -0.13 -0.23\n"
        "0.4 -0.31 0.22 -0.13 1 0.111\n"
        "0.5 -0.41 0.32 -0.23 0.111 1\n";
      ASSERT_EQUAL(s.str(),truth);
    }
  }

  void
  TestErrorPropagation()
  {
    vector<Parameter> input;
    input.push_back(Parameter(1.23,0.321));
    input.push_back(Parameter(4.56,0.654));
    input.push_back(Parameter(7.89,0.987));
    input.push_back(Parameter(2.01,0.234));

    CorrelationMatrix inputCorr(4);
    inputCorr.Set(0,1, 0.11);
    inputCorr.Set(0,2,-0.12);
    inputCorr.Set(0,3,-0.13);
    inputCorr.Set(1,2, 0.22);
    inputCorr.Set(1,3,-0.23);
    inputCorr.Set(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;
    
    NumericalErrorPropagation result(fcn,input,inputCorr);

    // 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]*inputCorr.Get(k,j)*input[k].error*input[j].error;
    }

    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 << result.GetPropagatedCorrelations().Get(i,j)*
            result.GetPropagatedErrors()[i]*
            result.GetPropagatedErrors()[j] << (j<2 ? " " : "\n");
    }
    
    ASSERT_EQUAL(s1.str(),s2.str());
  }

};


CPPUNIT_TEST_SUITE_REGISTRATION(testErrorPropagation);