source: trunk/test/averager_test.cc @ 499

// $Id: averager_test.cc 499 2006-01-18 16:26:09Z peter $

#include <c++_tools/statistics/Averager.h>
#include <c++_tools/statistics/AveragerPair.h>
#include <c++_tools/statistics/AveragerPairWeighted.h>
#include <c++_tools/statistics/AveragerWeighted.h>
#include <c++_tools/gslapi/vector.h>

#include <fstream>
#include <limits>
#include <iostream>

using namespace theplu::statistics;

//Forward declarations
bool equal(const Averager&, const Averager&);
bool equal(const AveragerWeighted&, const AveragerWeighted&, const double, 
           std::ostream* error);
bool equal(const Averager&, const AveragerWeighted&, const double, 
           std::ostream* error);
bool equal(const AveragerPair&, const AveragerPair&, 
           const double, std::ostream* error);
bool equal(const AveragerPair&, const AveragerPairWeighted&, 
           const double, std::ostream* error);
bool equal(const AveragerPairWeighted&, const AveragerPairWeighted&, 
           const double, std::ostream* error);


int main(const int argc,const char* argv[])
{  

  std::ostream* error;
  if (argc>1 && argv[1]==std::string("-v"))
    error = &std::cerr;
  else {
    error = new std::ofstream("/dev/null");
    if (argc>1)
      std::cout << "averager_test -v : for printing extra information\n";
  }
  bool ok = true;

  // Testing Averager
  *error << "testing Averager" << std::endl;
  Averager a;
  a.add(1);
  a.add(3);
  a.add(5);
  if (a.n()!=3 || a.mean()!=3 || a.sum_xx()!=35){
    ok=false;
    *error << "error: add\n";
  }

  Averager* a1 = new Averager(1+3+5,1+9+25,3);
  if (!equal(a,*a1)){
    ok=false;
    *error << "error: Averager(const double x,const double xx,const long n)\n";
  }
  delete a1;

  a1 = new Averager(a);
  if (!equal(a,*a1)){
    ok=false;
    *error << "error: Copy constructor\n";
  }
  delete a1;

  a.add(3,5);
  if (! a.standard_error()==sqrt(a.variance()/a.n())){
    ok=false;
    *error << "error: standard_error\n";
  }
      

  if ( fabs(a.variance() - a.std()*a.std())>
       std::numeric_limits<double>().round_error() ){
    ok=false;
    *error << "error: std squared should be variance" << std::endl;
    *error << "std2: " << a.std()*a.std() << std::endl;
    *error << "variance: " << a.variance() << std::endl;
    *error << "difference is: " << a.std()*a.std()-a.variance() << std::endl;
  }
  
  if ( a.variance() != a.variance(a.mean()) ){
    ok=false;
    *error << "error: variance incorrect\n" << std::endl;
    *error << "variance: " << a.variance() << std::endl;
    *error << "mean: " << a.mean() << std::endl;
    *error << "variance(mean) " << a.variance(a.mean()) << std::endl;
  }

  // Testing AveragerWeighted
  *error << "testing AveragerWeighted" << std::endl;
  theplu::gslapi::vector x(3,0);
  x(0)=0;
  x(1)=1;
  x(2)=2;
  theplu::gslapi::vector w(3,1);
  theplu::statistics::AveragerWeighted aw;
  aw.add(x,w);
  a.reset();
  a.add(x);
  const double tol=std::numeric_limits<double>().round_error();
  if (!equal(a,aw,tol,error)){
    *error << "error: AveragerWeighted with unitary weights should " 
           << "be equal to Averager" << std::endl;
    ok=false;
  }

  AveragerWeighted* aw2 = new AveragerWeighted(aw);
  if (!equal(aw,*aw2,tol,error)){
    *error << "error: AveragerWeighted copy constructor " << std::endl;
    ok=false;
  }
    
  aw2->add(12,0);
  if (!equal(aw,*aw2,tol,error)){
    *error << "error: AveragerWeighted adding a data point with weight=0 " 
           << "should make no change " << std::endl;
    ok=false;
  }
      
  aw2->reset();
  w.scale(17);
  aw2->add(x,w);
  if (!equal(aw,*aw2,tol,error)){
    *error << "error: AveragerWeighted rescaling weights " 
           << "should make no change " << std::endl;
    ok=false;
  }
  delete aw2;
  

  *error << "testing AveragerPair" << std::endl;
  AveragerPair ap;
  for (int i=0; i<10; i++)
    ap.add(i,i);
  if (fabs(ap.correlation()-1)>tol){
    ok=false;
    *error << "correlation: " << ap.correlation() << std::endl;
    *error << "error: correlation between identical vectors should be unity" 
           << std::endl;
  }
  if (ap.x_averager().variance()!=ap.covariance()){
    ok=false;
    *error << "error: covariance of identical vectors should equal to variance" 
           << std::endl;
  }
  AveragerPair* ap2 = new AveragerPair(ap);
  delete ap2;

  *error << "testing AveragerPairWeighted" << std::endl;
  AveragerPairWeighted apw;
  x(0)=0; x(1)=1; x(2)=2;
  theplu::gslapi::vector y(3,0);
  x(0)=0; x(1)=0; x(2)=2;
  apw.add(x,y,w,w);
  ap.reset();
  ap.add(x,y);
  if (!equal(ap,apw,tol,error)){
    *error << "error: AveragerPairWeighted with unitary weights should " 
           << "be equal to AveragerPair" << std::endl;
    ok=false;
  }

  AveragerPairWeighted* apw2 = new AveragerPairWeighted(apw);
  if (!equal(apw,*apw2,tol,error)){
    *error << "error: AveragerPairWeighted copy constructor " << std::endl;
    ok=false;
  }
    
  apw2->add(12,23222.03,32.3,0);
  if (!equal(apw,*apw2,tol,error)){
    *error << "error: AveragerWeighted adding a data point with weight=0 " 
           << "should make no change " << std::endl;
    ok=false;
  }
      
  apw2->reset();
  w.scale(17);
  apw2->add(x,y,w,w);
  if (!equal(apw,*apw2,tol,error)){
    *error << "error: AveragerWeighted rescaling weights " 
           << "should make no change " << std::endl;
    ok=false;
  }
  delete apw2;

  if (error!=&std::cerr)
    delete error;

  if (!ok)
    return -1;
  return 0;
}

bool equal(const Averager& a, const Averager& b)
{
  return (a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance());
}

bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
           const double tol, std::ostream* error)
{
  bool equal = true;
  if  ( fabs(a.mean()-b.mean())>tol){
    equal=false;
    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
  }
  if ( fabs(a.variance()-b.variance())>tol ) {
    equal=false;
    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
           << std::endl;
  } 
  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
    equal =false;
    *error << "error for standard error:\t" << std::endl;
  }
  return equal;
}

bool equal(const Averager& a, const AveragerWeighted& b, const double tol, 
           std::ostream* error)
{
  bool equal = true;
  if  ( fabs(a.mean()-b.mean())>tol){
    equal=false;
    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
  }
  if ( fabs(a.variance()-b.variance())>tol ) {
    equal=false;
    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
           << std::endl;
  } 
  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
    equal =false;
    *error << "error for standard error:\t" << std::endl;
  }
  return equal;
}

bool equal(const AveragerPair& a, const AveragerPair& b, 
           const double tol, std::ostream* error)
{
  bool ok = true;
  if  ( fabs(a.covariance()-b.covariance())>tol){
    ok=false;
    *error << "error covariance: " << a.covariance() << "\t" 
           << b.covariance() << std::endl;
  }
  if ( fabs(a.correlation()-b.correlation())>tol ) {
    ok=false;
    *error << "error correlation" << std::endl;
  } 
  return ok;
}

bool equal(const AveragerPair& a, const AveragerPairWeighted& b, 
           const double tol, std::ostream* error)
{
  bool ok = true;
  if  ( fabs(a.covariance()-b.covariance())>tol){
    ok=false;
    *error << "error covariance: " << a.covariance() << "\t" 
           << b.covariance() << std::endl;
  }
  if ( fabs(a.correlation()-b.correlation())>tol ) {
    ok=false;
    *error << "error correlation" << std::endl;
  } 
  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
    ok =false;
    *error << "error for x_averager():\t" << std::endl;
  }
  return ok;
}
bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
           const double tol, std::ostream* error)
{
  bool ok = true;
  if  ( fabs(a.covariance()-b.covariance())>tol){
    ok=false;
    *error << "error covariance: " << a.covariance() << "\t" 
           << b.covariance() << std::endl;
  }
  if ( fabs(a.correlation()-b.correlation())>tol ) {
    ok=false;
    *error << "error correlation" << std::endl;
  } 
  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
    ok =false;
    *error << "error for x_averager():\t" << std::endl;
  }
  return ok;
}