source: trunk/test/averager_test.cc @ 2119

// $Id: averager_test.cc 2119 2009-12-12 23:11:43Z peter $

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

  This file is part of the yat library, http://dev.thep.lu.se/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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include "Suite.h"

#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;
using namespace theplu::yat::statistics;

using theplu::yat::test::Suite;

//Forward declarations
bool equal(const Averager&, const Averager&, unsigned int, Suite&);
bool equal(const AveragerWeighted&, const AveragerWeighted&, unsigned int, 
           Suite& suite);
bool equal(const Averager&, const AveragerWeighted&, unsigned int, Suite& suite);
bool equal(const AveragerPair&, const AveragerPair&, unsigned int, Suite& suite);
bool equal(const AveragerPair&, const AveragerPairWeighted&, unsigned int, 
           Suite& suite);
bool equal(const AveragerPairWeighted&, const AveragerPairWeighted&,
           unsigned int, Suite& suite);


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

  Suite suite(argc, argv);

  // Testing Averager
  suite.out() << "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){
    suite.add(false);
    suite.err() << "error: add\n";
  }
  Averager b(a);
  b.add(5,-1);
  if (b.n()!=2 || b.mean()!=2 || b.sum_xx()!=10){
    suite.add(false);
    suite.err() << "error: add with negative n\n";
  }

  Averager* a1 = new Averager(1.0+3+5,1.0+9+25,3);
  unsigned int tol = 10;
  if (!equal(a,*a1, tol, suite)){
    suite.add(false);
    suite.err().precision(25);
    suite.err() << a.sum_x() << '\t' << a1->sum_x() << std::endl;
    suite.err() << a.sum_xx() << '\t' << a1->sum_xx() << std::endl;
    suite.err() << a.n() << '\t' << a1->n() << std::endl;
    suite.err() << a.variance() << '\t' << a1->variance() << std::endl;
    suite.err() << a.mean() << '\t' << a1->mean() << std::endl;
    suite.err() << a.mean() - a1->mean() << std::endl;
    suite.err() << a.variance() - a1->variance() << std::endl;
    suite.err() << "error: Averager(double x, double xx, long n)\n";
  }
  delete a1;

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

  a.add(3,5);
  if (!suite.equal_sqrt(a.standard_error(), sqrt(a.variance()/a.n()),1)) {
    suite.add(false);
    suite.err() << "error: standard_error\n";
  }
      
  if ( !suite.equal(a.variance(), a.std()*a.std()) ) {
    suite.add(false);
    suite.err() << "error: std squared should be variance" << std::endl;
    suite.err() << "std2: " << a.std()*a.std() << std::endl;
    suite.err() << "variance: " << a.variance() << std::endl;
    suite.err() << "difference is: " << a.std()*a.std()-a.variance() << std::endl;
  }
  
  if ( a.variance() != a.variance(a.mean()) ){
    suite.add(false);
    suite.err() << "error: variance incorrect\n" << std::endl;
    suite.err() << "variance: " << a.variance() << std::endl;
    suite.err() << "mean: " << a.mean() << std::endl;
    suite.err() << "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
  suite.err() << "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());
  if (!equal(a,aw,tol,suite)){
    suite.err() << "error: AveragerWeighted with unitary weights should " 
           << "be equal to Averager" << std::endl;
    suite.add(false);
  }

  AveragerWeighted* aw2 = new AveragerWeighted(aw);
  if (!equal(aw,*aw2,tol,suite)){
    suite.err() << "error: AveragerWeighted copy constructor " << std::endl;
    suite.add(false);
  }
    
  aw2->add(12,0);
  if (!equal(aw,*aw2,tol,suite)){
    suite.err() << "error: AveragerWeighted adding a data point with weight=0 " 
           << "should make no change " << std::endl;
    suite.add(false);
  }
      
  aw2->reset();
  w*=17;
  add(*aw2, x.begin(), x.end(), w.begin());
  if (!equal(aw,*aw2,tol,suite)){
    suite.err() << "error: AveragerWeighted rescaling weights " 
           << "should make no change " << std::endl;
    suite.add(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());
  }
  

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

  for (int i=0; i<8; i++)
    ap.add(static_cast<double>(i),i,-1);
  if (!suite.equal(ap.correlation(),1, tol)) {
    suite.add(false);
    suite.err() << "correlation after removal of data: " << ap.correlation()
                << std::endl;
    suite.err() << "error: correlation between identical vectors is unity" 
                << std::endl;
  }


  suite.out() << "testing AveragerPairWeighted" << std::endl;
  AveragerPairWeighted apw;
  x(0)=0; x(1)=1; x(2)=2;
  theplu::yat::utility::Vector y(3,0);
  y(0)=0; y(1)=0; y(2)=2;
  add(apw, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
  ap.reset();
  add(ap, x.begin(), x.end(), y.begin());
  if (!equal(ap,apw,tol,suite)){
    suite.err() << "error: AveragerPairWeighted with unitary weights should " 
           << "be equal to AveragerPair" << std::endl;
    suite.add(false);
  }

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

  return suite.return_value();
}

bool equal(const Averager& a, const Averager& b, unsigned int tol, 
           Suite& suite)
{
  return (a.n()==b.n() && suite.equal(a.mean(),b.mean(),tol) &&
          suite.equal(a.variance(),b.variance(),tol));
}

bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
           const unsigned int tol, Suite& suite)
{
  bool equal = true;
  if  ( !suite.equal(a.mean(),b.mean(),tol)){
    equal=false;
    suite.err() << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
    suite.err() << "difference:\t" << a.mean()-b.mean() << std::endl;
  }
  if ( !suite.equal(a.variance(),b.variance(),tol) ) {
    equal=false;
    suite.err() << "error for variance:\t" << a.variance() << " " 
                << b.variance() << std::endl;
  } 
  if ( !suite.equal(a.standard_error(),b.standard_error(),tol) ) {
    equal =false;
    suite.err() << "error for standard error:\t" << std::endl;
  }
  return equal;
}

bool equal(const Averager& a, const AveragerWeighted& b, const unsigned int tol, 
           Suite& suite)
{
  bool equal = true;
  if  ( !suite.equal(a.mean(),b.mean(),tol)){
    equal=false;
    suite.err() << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
  }
  if ( !suite.equal(a.variance(),b.variance(),tol) ) {
    equal=false;
    suite.err() << "error for variance:\t" << a.variance() << " " 
                << b.variance() << std::endl;
  } 
  if ( !suite.equal(a.standard_error(), b.standard_error(),tol) ) {
    equal =false;
    suite.err() << "error for standard error:\t" << std::endl;
  }
  return equal;
}

bool equal(const AveragerPair& a, const AveragerPair& b, 
           const unsigned int tol, Suite& suite)
{
  bool ok = true;
  if (!suite.equal(a.covariance(), b.covariance(), tol)) {
    suite.add(false);
    suite.err() << "error covariance: " << a.covariance() << "\t" 
           << b.covariance() << std::endl;
  }
  if ( !suite.equal(a.correlation(),b.correlation(), tol)) {
    suite.add(false);
    suite.err() << "error correlation" << std::endl;
  } 
  return ok;
}

bool equal(const AveragerPair& a, const AveragerPairWeighted& b, 
           const unsigned int tol, Suite& suite)
{
  bool ok = true;
  if (!suite.equal(a.covariance(), b.covariance(), tol) ) {
    suite.add(false);
    suite.err() << "error covariance: " << a.covariance() << "\t" 
           << b.covariance() << std::endl;
  }
  if ( !suite.equal(a.correlation(),b.correlation(), tol) ) {
    suite.add(false);
    suite.err() << "error correlation" << std::endl;
    suite.err() << "unweighted:" << a.correlation() << std::endl;
    suite.err() << "weighted:" << b.correlation() << std::endl;
    suite.err() << "difference:" << a.correlation()-b.correlation() 
                << std::endl;
  } 
  if ( !equal(a.x_averager(),b.x_averager(),tol,suite)) {
    ok =false;
    suite.err() << "error for x_averager():\t" << std::endl;
  }
  return ok;
}
bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
           const unsigned int tol, Suite& suite)
{
  bool ok = true;
  if  ( !suite.equal(a.covariance(),b.covariance(),tol) ){
    suite.add(false);
    suite.err() << "error covariance: " << a.covariance() << "\t" 
           << b.covariance() << std::endl;
  }
  if ( !suite.equal(a.correlation(), b.correlation(), tol) ) {
    suite.add(false);
    suite.err() << "error correlation" << std::endl;
    suite.err() << "a:" << a.correlation() << std::endl;
    suite.err() << "b:" << b.correlation() << std::endl;
    suite.err() << "difference:" << a.correlation()-b.correlation() 
                << std::endl;
    suite.err() << "tol:" << tol << std::endl;
  } 
  if ( !equal(a.x_averager(),b.x_averager(),tol,suite)) {
    ok =false;
    suite.err() << "error for x_averager():\t" << std::endl;
  }
  return ok;
}