testTimeDistribution.cc

Go to the documentation of this file.

#include <tst/Verify.h>

#include <utl/AugerUnits.h>

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

#include <utl/TimeDistribution.h>
#include <utl/MultiTimeDistribution.h>

#include <utl/TimeDistributionAlgorithm.h>

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


class testTimeDistribution : public CppUnit::TestFixture {

  CPPUNIT_TEST_SUITE(testTimeDistribution);

  CPPUNIT_TEST(testInitialize);
  CPPUNIT_TEST(testBinSizes);
  CPPUNIT_TEST(testOperatorBracket);
  CPPUNIT_TEST(testSparsness);
  CPPUNIT_TEST(testAssignment);
  CPPUNIT_TEST_EXCEPTION(testSumException, AugerException);
  CPPUNIT_TEST(testEqual);
  CPPUNIT_TEST(testOperatorArithmetic);

  //CPPUNIT_TEST(testSlotIterator);

  CPPUNIT_TEST(testSparseIterator);
  CPPUNIT_TEST(testMultiTimeDistribution); 

  CPPUNIT_TEST(testAlgorithms);
  
  CPPUNIT_TEST_SUITE_END();

private:
  TimeDistributionD* td1;
  TimeDistributionD* td2;
  TimeDistributionD* td2weighted;
  TimeDistributionD* td3;
  TimeDistributionI* td4;
  TimeDistributionD* td5;
  TimeDistributionD* td6;

public:
  void
  setUp()
  {
    // add data by slot number
    td1 = new TimeDistributionD(1*nanosecond);    
    for (int i = -10 ; i <= 10 ; ++i)
      td1->AddTime(i, double(i));

    // add data by specifying time.  This should put 3 entries
    // in each of 10 (1 ns) time slots
    td2 = new TimeDistributionD(1*nanosecond);
    for (int i = 0 ; i < 30 ; ++i) {
      const double time = (1/3.0)*i - 5 + 1e-6;
      td2->AddTime(time);
    }

    td2weighted = new TimeDistributionD(1*nanosecond);
    for (int i = 0 ; i < 30 ; ++i) {
      const double time = (1/3.0)*i - 5 + 1e-6;
      td2weighted->AddTime(time, 2.0);
    }

    // fill a distribution from -10,0 and 11,20.  leave 1-10 empty.
    // This distribution cand then be used to check Iterator and 
    // SparseIterator
    td3 = new TimeDistributionD(1*nanosecond);
    for (int i = -10 ; i <=  0 ; ++i)  // note i=0 does not get stored
      td3->AddTime(i, double(i));
    for (int i =  11 ; i <= 20 ; ++i)
      td3->AddTime(i, double(i));


    // time distribution containing ints
    td4 = new TimeDistributionI(1*nanosecond);    
    for (int i = -10 ; i <= 10 ; ++i)
      td4->AddTime(i, i);

    // fill a distribution from 3-5, and 25-30
    td5 = new TimeDistributionD(1*nanosecond);
    for (int i = 3 ; i <= 5 ; ++i)
      td5->AddTime(i, double(i));
    for (int i = 25; i <= 30; ++i)
      td5->AddTime(i, double(i));

    // use a different binning
    td6 = new TimeDistributionD(2*nanosecond);
    for (int i = 3 ; i <= 5 ; ++i)
      td6->AddTime(i, double(i));
  }

  void
  tearDown()
  {
    delete td1;
    delete td2;
    delete td2weighted;
    delete td3;
    delete td4;
    delete td5;
    delete td6;
  }

  void
  testBinSizes()
  {
    CPPUNIT_ASSERT(Verify<CloseTo>(td1->GetBinning(), 1*nanosecond));
  }

  void
  testInitialize()
  {  
    for (TimeDistributionD::SparseIterator tdi = td1->SparseBegin();
         tdi != td1->SparseEnd(); ++tdi)
      CPPUNIT_ASSERT(Verify<CloseTo>(double(tdi->get<0>()), tdi->get<1>()));
  }

  void
  testAssignment()
  {
    TimeDistributionD copyTd(35*nanosecond);
    copyTd = *td1;
    CPPUNIT_ASSERT(Verify<CloseTo>(copyTd.GetBinning(), 1*nanosecond));
    for (int i = -10 ; i <= 10 ; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>(copyTd[i], double(i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(copyTd.At(i), double(i)));
    }
  }

  void
  testOperatorBracket()
  {
    // defined slots
    for (int i = -10 ; i <= 10 ; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>((*td1)[i], double(i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(td1->At(i), double(i)));
      CPPUNIT_ASSERT(Verify<Equal>((*td4)[i], i));
      CPPUNIT_ASSERT(Verify<Equal>(td4->At(i), i));

      (*td1)[i] += 2.0;
      CPPUNIT_ASSERT(Verify<CloseTo>((*td1)[i], double(i) + 2));
      CPPUNIT_ASSERT(Verify<CloseTo>(td1->At(i), double(i) + 2));

      (*td4)[i] += 2;
      CPPUNIT_ASSERT(Verify<Equal>((*td4)[i], i + 2));
      CPPUNIT_ASSERT(Verify<Equal>(td4->At(i), i + 2));
    }

    // undefined slots return 0
    for (int i = -20 ; i < -10 ; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>((*td1)[i], 0.0));      
      CPPUNIT_ASSERT(Verify<CloseTo>(td1->At(i), 0.0));      
      CPPUNIT_ASSERT(Verify<Equal>((*td4)[i], 0));
      CPPUNIT_ASSERT(Verify<Equal>(td4->At(i), 0));
      
      (*td1)[i] += 10.0;
      CPPUNIT_ASSERT(Verify<CloseTo>((*td1)[i], 10.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(td1->At(i), 10.0));
      (*td4)[i] += 10;
      CPPUNIT_ASSERT(Verify<Equal>((*td4)[i], 10));
      CPPUNIT_ASSERT(Verify<Equal>(td4->At(i), 10));
    }
    for (int i = 11 ; i < 20 ; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>((*td1)[i], 0.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(td1->At(i), 0.0));
      CPPUNIT_ASSERT(Verify<Equal>((*td4)[i], 0));
      CPPUNIT_ASSERT(Verify<Equal>(td4->At(i), 0));
    }

    // histos were filled by times (rather than slot indices) with 3 time entries per slot
    for (int i = -5 ; i < 4 ; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>((*td2)[i], 3.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(td2->At(i), 3.0));
    }

    // histos were filled by times (rather than slot indices) with 3 time entries per slot
    // and with weights of 2 on each entry
    for (int i = -5 ; i < 4 ; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>((*td2weighted)[i], 6.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(td2weighted->At(i), 6.0));
    }

    // look up histos by times (rather than slot indices)
    for (int i = 0 ; i < 30 ; ++i) {
      const double time = (1/3.0)*i - 5 + 1e-6;
      CPPUNIT_ASSERT(Verify<CloseTo>((*td2)[time], 3.0));
      CPPUNIT_ASSERT(Verify<CloseTo>(td2->At(time), 3.0));
    }

    // GetStart() and GetStop() methods
    CPPUNIT_ASSERT(Verify<Equal>(td3->GetStart(), -10));
    CPPUNIT_ASSERT(Verify<Equal>(td3->GetNumSlots(), 31));
    CPPUNIT_ASSERT(Verify<Equal>(td3->GetStop(), 20));

    TimeDistributionD tdEmpty(1*nanosecond);
    CPPUNIT_ASSERT(Verify<Equal>(tdEmpty.GetNumSlots(), 0));
    CPPUNIT_ASSERT(Verify<Equal>(tdEmpty.GetStart(), 0));
    CPPUNIT_ASSERT(Verify<Equal>(tdEmpty.GetStop(), -1));
  }

  void
  testSparsness()
  {
    {
      TimeDistributionI td(1*nanosecond);
      for (int i = 0; i < 10; ++i)
        td.SetTime(i, i);
      for (int i = 0; i < 10; ++i)
        CPPUNIT_ASSERT(Verify<Equal>(td.At(i), i));
      TimeDistributionI::SparseIterator first = td.SparseBegin();
      // 0 should not exist
      CPPUNIT_ASSERT(Verify<Equal>(first->get<0>(), 1));
      CPPUNIT_ASSERT(Verify<Equal>(first->get<1>(), 1));
      td.SetTime(1, 0);
      first = td.SparseBegin();
      // set 1 to 0 should remove the entry
      CPPUNIT_ASSERT(Verify<Equal>(first->get<0>(), 2));
      CPPUNIT_ASSERT(Verify<Equal>(first->get<1>(), 2));
      td.AddTime(2, -2);
      first = td.SparseBegin();
      // if the add cancels the entry should dissapear
      CPPUNIT_ASSERT(Verify<Equal>(first->get<0>(), 3));
      CPPUNIT_ASSERT(Verify<Equal>(first->get<1>(), 3));
      // setting all to zero
      for (int i = 0; i < 10; ++i)
        td.SetTime(i, 0);
      // no entries
      CPPUNIT_ASSERT((td.SparseBegin() == td.SparseEnd()));
    }

    { // repeat with double
      TimeDistributionD td(1*nanosecond);
      for (int i = 0; i < 10; ++i)
        td.SetTime(i, double(i));
      for (int i = 0; i < 10; ++i)
        CPPUNIT_ASSERT(Verify<CloseTo>(td.At(i), double(i)));
      TimeDistributionD::SparseIterator first = td.SparseBegin();
      CPPUNIT_ASSERT(Verify<CloseTo>(first->get<0>(), 1));
      CPPUNIT_ASSERT(Verify<CloseTo>(first->get<1>(), double(1)));
      td.SetTime(1, 0);
      first = td.SparseBegin();
      CPPUNIT_ASSERT(Verify<CloseTo>(first->get<0>(), 2));
      CPPUNIT_ASSERT(Verify<CloseTo>(first->get<1>(), double(2)));
      td.AddTime(2, -2);
      first = td.SparseBegin();
      CPPUNIT_ASSERT(Verify<CloseTo>(first->get<0>(), 3));
      CPPUNIT_ASSERT(Verify<CloseTo>(first->get<1>(), double(3)));
      for (int i = 0; i < 10; ++i)
        td.SetTime(i, 0);
      CPPUNIT_ASSERT((td.SparseBegin() == td.SparseEnd()));
    }
  }

  void
  testSumException()
  {
    tst::Expected();
    // two distributions of unequal size
    TimeDistributionD tdBarf = *td3 + *td6;
  }

  void
  testOperatorArithmetic()
  {
    // operator+
    TimeDistributionD tdSum = *td3 + *td5;

    // check that only bins with data in them were affected
    int sum = 0;
    for (TimeDistributionD::SparseIterator sit = tdSum.SparseBegin();
         sit != tdSum.SparseEnd(); ++sit, ++sum)
      CPPUNIT_ASSERT(Verify<Equal>(sit->get<0>(), int(sit->get<1>())));

    CPPUNIT_ASSERT(Verify<Equal>(sum, 30-1));  // confirm empty bins were skipped

    // check that start and num slots properly set
    for (int i = tdSum.GetStart(); i != tdSum.GetStop(); ++i)
      if (i <= 0 || (i >= 11 && i <= 20) || (i >= 3 && i <= 5) ||
          (i >= 30 && i <= 35)) {
        CPPUNIT_ASSERT(Verify<Equal>(i, int(tdSum[i])));
        CPPUNIT_ASSERT(Verify<Equal>(i, int(tdSum.At(i))));
      }

    // add distribution to itself + check the sum
    TimeDistributionD tdSum2 = *td6 + *td6;
    for (int i = tdSum2.GetStart(); i != tdSum2.GetStop(); ++i) {
      CPPUNIT_ASSERT(Verify<Equal>(2*i, int(tdSum2[i])));
      CPPUNIT_ASSERT(Verify<Equal>(2*i, int(tdSum2.At(i))));
    }
  }

  void
  testSparseIterator()
  {
    // dereferencing and preincrement
    int i = -10;
    for (TimeDistributionD::SparseIterator sit = td3->SparseBegin();
         sit != td3->SparseEnd(); ++sit, ++i) {
      if (i == 0) // this distribution should skip bins 0-10
        i = 11;
      CPPUNIT_ASSERT((*sit == boost::make_tuple<const int, const double>(i, i)));
    }

    // operator-> and and postincrement
    i = -10;
    for (TimeDistributionD::SparseIterator sit = td3->SparseBegin();
         sit != td3->SparseEnd(); ++sit, ++i) {
      if (i == 0) // this distribution should skip bins 0-10
        i = 11;
      CPPUNIT_ASSERT(sit->get<0>() == i);
      CPPUNIT_ASSERT(sit->get<1>() == i);
    }
    
    // assignment
    i = -10;
    for (TimeDistributionD::SparseIterator sit = td3->SparseBegin();
         sit != td3->SparseEnd(); ++sit, ++i) {
      if (i == 0) // this distribution should skip bins 0-10
        i = 11;
      TimeDistributionD::SparseIterator copyIt = sit;
      CPPUNIT_ASSERT((*copyIt == boost::make_tuple<const int, const double>(i, i)));
    }
  }

  void
  testMultiTimeDistribution()
  {  
    MultiTimeDistributionD mtd;
    mtd.AddTimeDistribution(1*nanosecond, 0);
    TimeDistributionD& td = mtd.GetTimeDistribution(0);
    for (int i = 0; i < 3; ++i)
      td[i] = double(i);
    for (int i = 0; i < 3; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>(td[i], double(i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(td.At(i), double(i)));
    }
    for (int i = 0; i < 3; ++i)
      td.SetTime(i, double(i));
    for (int i = 0; i < 3; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>(td[i], double(i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(td.At(i), double(i)));
    }
    td.Clear();
    for (int i = 0; i < 3; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>(td[i], 0.));
      CPPUNIT_ASSERT(Verify<CloseTo>(td.At(i), 0.));
    }
    for (int i = 0; i < 3; ++i)
      td.SetTime(i, double(i));
    for (int i = 0; i < 3; ++i) {
      CPPUNIT_ASSERT(Verify<CloseTo>(td[i], double(i)));
      CPPUNIT_ASSERT(Verify<CloseTo>(td.At(i), double(i)));
    }

    // HasLabel
    CPPUNIT_ASSERT(Verify<Equal>(mtd.HasLabel(0), true));
    CPPUNIT_ASSERT(Verify<Equal>(mtd.HasLabel(1), false));

    // Add more time distributions
    mtd.AddTimeDistribution(*td1, 1);
    mtd.AddTimeDistribution(*td1, 2);
    CPPUNIT_ASSERT(Verify<Equal>(mtd.GetNLabels(), 3U));

    // multiTimeDistribution accessor
    TimeDistributionD tdSource = mtd.GetTimeDistribution(1);
    CPPUNIT_ASSERT(Verify<CloseTo>(tdSource[-5], -5.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(tdSource[5], 5.0));

    // multiTimeDistribution iterator
    int i = 0;
    for (MultiTimeDistributionD::Iterator iter = mtd.Begin();
         iter != mtd.End(); ++iter, ++i)
      CPPUNIT_ASSERT(Verify<Equal>(iter->GetLabel(), i));

    // copy ctor
    MultiTimeDistributionD mcopy = mtd;
    i = 0;
    for (MultiTimeDistributionD::Iterator iter = mcopy.Begin();
         iter != mcopy.End(); ++iter, ++i) {
      CPPUNIT_ASSERT(Verify<Equal>(iter->GetLabel(), i));

      TimeDistributionD ta = iter->GetTimeDistribution();
      for (int j = 0 ; j < 2 ; ++j)
        CPPUNIT_ASSERT(Verify<Equal>(ta[j], double(j)));
    }

    // test a MultiTimeDistribution for ints
    MultiTimeDistributionI mtdi;
    mtdi.AddTimeDistribution(1*nanosecond, 0);
    TimeDistributionI& tdi = mtdi.GetTimeDistribution(0);
    for (int i = 0; i < 3; ++i)
      tdi[i] = i;

    // copy ctor
    TimeDistributionI tdiSource = mtdi.GetTimeDistribution(0);
    CPPUNIT_ASSERT(Verify<Equal>(tdiSource[0], 0));
    CPPUNIT_ASSERT(Verify<Equal>(tdiSource[2], 2));
  }
  
  void
  testAlgorithms()
  {
    TimeDistributionD td(1*nanosecond);
    
    td[-1] = 1.0; 
    td[0] = 4.0;
    // skip bin 1
    td[2] = 4.0;
    td[3] = 5.0;
    td[4] = 6.0;
    // skip 5
    td[6] = 0.;  // explicitly set 6
    td[7] = 10.;

    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, -1, 7), 30.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, -100, 100), 30.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, 2, 4), 15.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, 4, 2), 15.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, 100, -100), 30.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, 50, 100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, 100, 50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, -100, -50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, -50, -100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(td, 2, 2), 4.0));

    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, -1, 2), 9.0/4.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, -100, 100), 30.0/201.));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, 100, -100), 30.0/201.));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, -100, -50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, 50, 100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, 100, 50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, -100, -50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, -50, -100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Mean(td, 2, 2), 4.0));

    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, 2, 4), sqrt(2.0/3.0)));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, 4, 2), sqrt(2.0/3.0)));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, 0, 2), sqrt(32./9.)));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, 50, 100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, 100, 50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, -100, -50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, -50, -100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::RMS(td, 2, 2), 0.0));

    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, -1, 2), 7.0/9.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, -100, 2), 7.0/9.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, 50, 100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, 100, 50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, -100, -50), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, -50, -100), 0.0));
    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Centroid(td, 2, 2), 2.0));

    TimeDistributionI tdi(1*nanosecond);
    
    tdi[-1] = 1; 
    tdi[0] = 4;
    // skip bin 1
    tdi[2] = 4;
    tdi[3] = 5;
    tdi[4] = 6;
    // skip 5
    tdi[6] = 0;  // explicitly set 6
    tdi[7] = 10;

    CPPUNIT_ASSERT(Verify<CloseTo>(TimeDistributionAlgorithm::Sum(tdi, -1, 7), 30.0));
    CPPUNIT_ASSERT(Verify<CloseTo> (TimeDistributionAlgorithm::Mean(tdi, -1, 2), 9.0/4.0));
    CPPUNIT_ASSERT(Verify<CloseTo> (TimeDistributionAlgorithm::RMS(td, 2, 4), sqrt(2.0/3.0)));
    CPPUNIT_ASSERT(Verify<CloseTo> (TimeDistributionAlgorithm::Centroid(td, -1, 2), 7.0/9.0));
  }

  void
  testEqual()
  {
    TimeDistributionI tdi(*td4);
    CPPUNIT_ASSERT(tdi == *td4);
    ++tdi[0];
    CPPUNIT_ASSERT(tdi != *td4);

    TimeDistributionD tdd(*td1);
    CPPUNIT_ASSERT(tdd == *td1);
    ++tdd[1];
    CPPUNIT_ASSERT(tdd != *td1);
  }
  
};


CPPUNIT_TEST_SUITE_REGISTRATION(testTimeDistribution);