source: trunk/test/averager_test.cc @ 916

// $Id: averager_test.cc 916 2007-09-30 00:50:10Z peter $

/*
  Copyright (C) 2005 Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
  Copyright (C) 2007 Jari Häkkinen

  This file is part of the yat library, http://trac.thep.lu.se/trac/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "yat/statistics/Averager.h"
#include "yat/statistics/AveragerPair.h"
#include "yat/statistics/AveragerPairWeighted.h"
#include "yat/statistics/AveragerWeighted.h"
#include "yat/utility/vector.h"

#include <cmath>
#include <fstream>
#include <limits>
#include <iostream>
#include <set>

using namespace theplu::yat::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.0+3+5,1.0+9+25,3);
  if (!equal(a,*a1)){
    ok=false;
    std::cout.precision(25);
    std::cout << equal(a,*a1) << std::endl;
    std::cout << a.sum_x() << '\t' << a1->sum_x() << std::endl;
    std::cout << a.sum_xx() << '\t' << a1->sum_xx() << std::endl;
    std::cout << a.n() << '\t' << a1->n() << std::endl;
    std::cout << a.variance() << '\t' << a1->variance() << std::endl;
    std::cout << a.mean() << '\t' << a1->mean() << std::endl;
    std::cout << a.mean() - a1->mean() << std::endl;
    std::cout << a.variance() - a1->variance() << std::endl;
    *error << "error: Averager(double x, double xx, u_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;
  }
  theplu::yat::utility::vector* tmp_vec = new theplu::yat::utility::vector(10);
  add(a, tmp_vec->begin(), tmp_vec->end());
  delete tmp_vec;
  std::set<double>* tmp_set = new std::set<double>;
  tmp_set->insert(1.0);
  tmp_set->insert(2.0);
  add(a, tmp_set->begin(), tmp_set->end());
  delete tmp_set;


  // Testing AveragerWeighted
  *error << "testing AveragerWeighted" << std::endl;
  theplu::yat::utility::vector x(3,0);
  x(0)=0;
  x(1)=1;
  x(2)=2;
  theplu::yat::utility::vector w(3,1);
  theplu::yat::statistics::AveragerWeighted aw;
  add(aw, x.begin(), x.end(), w.begin());
  a.reset();
  add(a, x.begin(), x.end());
  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*=17;
  add(*aw2, x.begin(), x.end(), w.begin());
  if (!equal(aw,*aw2,tol,error)){
    *error << "error: AveragerWeighted rescaling weights " 
           << "should make no change " << std::endl;
    ok=false;
  }
  delete aw2;
  {
    theplu::yat::utility::vector tmp(10);
    add(aw, tmp.begin(), tmp.end());
  }
  {
    std::set<double> tmp;
    tmp.insert(1.0);
    tmp.insert(2.0);
    add(aw, tmp.begin(), tmp.end());
  }
  

  *error << "testing AveragerPair" << std::endl;
  AveragerPair ap;
  for (int i=0; i<10; i++)
    ap.add(static_cast<double>(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::yat::utility::vector y(3,0);
  x(0)=0; x(1)=0; x(2)=2;
  apw.add_values(x,y,w,w);
  ap.reset();
  ap.add_values(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*=17;
  apw2->add_values(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)
{
//  std::cout << (a.n()==b.n()) << std::endl;
//  std::cout << (a.mean()==b.mean()) << std::endl;
//  std::cout << (fabs(a.variance()-b.variance()<1e-15)) << std::endl;
  return (a.n()==b.n() && a.mean()==b.mean() &&
          fabs(a.variance()-b.variance()<1e-15));
}

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;
}