source: trunk/test/averager_test.cc @ 1120

Last change on this file since 1120 was 1120, checked in by Peter, 14 years ago

vector is now Vector

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 10.0 KB
Line 
1// $Id: averager_test.cc 1120 2008-02-21 23:18:41Z peter $
2
3/*
4  Copyright (C) 2005 Peter Johansson
5  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
6  Copyright (C) 2007 Jari Häkkinen
7
8  This file is part of the yat library, http://trac.thep.lu.se/yat
9
10  The yat library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU General Public License as
12  published by the Free Software Foundation; either version 2 of the
13  License, or (at your option) any later version.
14
15  The yat library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  General Public License for more details.
19
20  You should have received a copy of the GNU General Public License
21  along with this program; if not, write to the Free Software
22  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
23  02111-1307, USA.
24*/
25
26#include "yat/statistics/Averager.h"
27#include "yat/statistics/AveragerPair.h"
28#include "yat/statistics/AveragerPairWeighted.h"
29#include "yat/statistics/AveragerWeighted.h"
30#include "yat/utility/Vector.h"
31
32#include <cmath>
33#include <fstream>
34#include <limits>
35#include <iostream>
36#include <set>
37
38using namespace theplu::yat::statistics;
39
40//Forward declarations
41bool equal(const Averager&, const Averager&);
42bool equal(const AveragerWeighted&, const AveragerWeighted&, const double, 
43           std::ostream* error);
44bool equal(const Averager&, const AveragerWeighted&, const double, 
45           std::ostream* error);
46bool equal(const AveragerPair&, const AveragerPair&, 
47           const double, std::ostream* error);
48bool equal(const AveragerPair&, const AveragerPairWeighted&, 
49           const double, std::ostream* error);
50bool equal(const AveragerPairWeighted&, const AveragerPairWeighted&, 
51           const double, std::ostream* error);
52
53
54int main(const int argc,const char* argv[])
55{ 
56
57  std::ostream* error;
58  if (argc>1 && argv[1]==std::string("-v"))
59    error = &std::cerr;
60  else {
61    error = new std::ofstream("/dev/null");
62    if (argc>1)
63      std::cout << "averager_test -v : for printing extra information\n";
64  }
65  bool ok = true;
66
67  // Testing Averager
68  *error << "testing Averager" << std::endl;
69  Averager a;
70  a.add(1);
71  a.add(3);
72  a.add(5);
73  if (a.n()!=3 || a.mean()!=3 || a.sum_xx()!=35){
74    ok=false;
75    *error << "error: add\n";
76  }
77
78
79  Averager* a1 = new Averager(1.0+3+5,1.0+9+25,3);
80  if (!equal(a,*a1)){
81    ok=false;
82    std::cout.precision(25);
83    std::cout << equal(a,*a1) << std::endl;
84    std::cout << a.sum_x() << '\t' << a1->sum_x() << std::endl;
85    std::cout << a.sum_xx() << '\t' << a1->sum_xx() << std::endl;
86    std::cout << a.n() << '\t' << a1->n() << std::endl;
87    std::cout << a.variance() << '\t' << a1->variance() << std::endl;
88    std::cout << a.mean() << '\t' << a1->mean() << std::endl;
89    std::cout << a.mean() - a1->mean() << std::endl;
90    std::cout << a.variance() - a1->variance() << std::endl;
91    *error << "error: Averager(double x, double xx, u_long n)\n";
92  }
93  delete a1;
94
95  a1 = new Averager(a);
96  if (!equal(a,*a1)){
97    ok=false;
98    *error << "error: Copy constructor\n";
99  }
100  delete a1;
101
102  a.add(3,5);
103  if (! a.standard_error()==sqrt(a.variance()/a.n())){
104    ok=false;
105    *error << "error: standard_error\n";
106  }
107     
108
109  if ( fabs(a.variance() - a.std()*a.std())>
110       std::numeric_limits<double>().round_error() ){
111    ok=false;
112    *error << "error: std squared should be variance" << std::endl;
113    *error << "std2: " << a.std()*a.std() << std::endl;
114    *error << "variance: " << a.variance() << std::endl;
115    *error << "difference is: " << a.std()*a.std()-a.variance() << std::endl;
116  }
117 
118  if ( a.variance() != a.variance(a.mean()) ){
119    ok=false;
120    *error << "error: variance incorrect\n" << std::endl;
121    *error << "variance: " << a.variance() << std::endl;
122    *error << "mean: " << a.mean() << std::endl;
123    *error << "variance(mean) " << a.variance(a.mean()) << std::endl;
124  }
125  theplu::yat::utility::Vector* tmp_vec = new theplu::yat::utility::Vector(10);
126  add(a, tmp_vec->begin(), tmp_vec->end());
127  delete tmp_vec;
128  std::set<double>* tmp_set = new std::set<double>;
129  tmp_set->insert(1.0);
130  tmp_set->insert(2.0);
131  add(a, tmp_set->begin(), tmp_set->end());
132  delete tmp_set;
133
134
135  // Testing AveragerWeighted
136  *error << "testing AveragerWeighted" << std::endl;
137  theplu::yat::utility::Vector x(3,0);
138  x(0)=0;
139  x(1)=1;
140  x(2)=2;
141  theplu::yat::utility::Vector w(3,1);
142  theplu::yat::statistics::AveragerWeighted aw;
143  add(aw, x.begin(), x.end(), w.begin());
144  a.reset();
145  add(a, x.begin(), x.end());
146  const double tol=std::numeric_limits<double>().round_error();
147  if (!equal(a,aw,tol,error)){
148    *error << "error: AveragerWeighted with unitary weights should " 
149           << "be equal to Averager" << std::endl;
150    ok=false;
151  }
152
153  AveragerWeighted* aw2 = new AveragerWeighted(aw);
154  if (!equal(aw,*aw2,tol,error)){
155    *error << "error: AveragerWeighted copy constructor " << std::endl;
156    ok=false;
157  }
158   
159  aw2->add(12,0);
160  if (!equal(aw,*aw2,tol,error)){
161    *error << "error: AveragerWeighted adding a data point with weight=0 " 
162           << "should make no change " << std::endl;
163    ok=false;
164  }
165     
166  aw2->reset();
167  w*=17;
168  add(*aw2, x.begin(), x.end(), w.begin());
169  if (!equal(aw,*aw2,tol,error)){
170    *error << "error: AveragerWeighted rescaling weights " 
171           << "should make no change " << std::endl;
172    ok=false;
173  }
174  delete aw2;
175  {
176    theplu::yat::utility::Vector tmp(10);
177    add(aw, tmp.begin(), tmp.end());
178  }
179  {
180    std::set<double> tmp;
181    tmp.insert(1.0);
182    tmp.insert(2.0);
183    add(aw, tmp.begin(), tmp.end());
184  }
185 
186
187  *error << "testing AveragerPair" << std::endl;
188  AveragerPair ap;
189  for (int i=0; i<10; i++)
190    ap.add(static_cast<double>(i),i);
191  if (fabs(ap.correlation()-1)>tol){
192    ok=false;
193    *error << "correlation: " << ap.correlation() << std::endl;
194    *error << "error: correlation between identical vectors should be unity" 
195           << std::endl;
196  }
197  if (ap.x_averager().variance()!=ap.covariance()){
198    ok=false;
199    *error << "error: covariance of identical vectors should equal to variance" 
200           << std::endl;
201  }
202  AveragerPair* ap2 = new AveragerPair(ap);
203  delete ap2;
204
205  *error << "testing AveragerPairWeighted" << std::endl;
206  AveragerPairWeighted apw;
207  x(0)=0; x(1)=1; x(2)=2;
208  theplu::yat::utility::Vector y(3,0);
209  x(0)=0; x(1)=0; x(2)=2;
210  add(apw, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
211  ap.reset();
212  add(ap, x.begin(), x.end(), y.begin());
213  if (!equal(ap,apw,tol,error)){
214    *error << "error: AveragerPairWeighted with unitary weights should " 
215           << "be equal to AveragerPair" << std::endl;
216    ok=false;
217  }
218
219  AveragerPairWeighted* apw2 = new AveragerPairWeighted(apw);
220  if (!equal(apw,*apw2,tol,error)){
221    *error << "error: AveragerPairWeighted copy constructor " << std::endl;
222    ok=false;
223  }
224   
225  apw2->add(12,23222.03,32.3,0);
226  if (!equal(apw,*apw2,tol,error)){
227    *error << "error: AveragerWeighted adding a data point with weight=0 " 
228           << "should make no change " << std::endl;
229    ok=false;
230  }
231     
232  apw2->reset();
233  w*=17;
234  add(*apw2, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
235  if (!equal(apw,*apw2,tol,error)){
236    *error << "error: AveragerWeighted rescaling weights " 
237           << "should make no change " << std::endl;
238    ok=false;
239  }
240  delete apw2;
241
242  if (error!=&std::cerr)
243    delete error;
244
245  if (!ok)
246    return -1;
247  return 0;
248}
249
250bool equal(const Averager& a, const Averager& b)
251{
252//  std::cout << (a.n()==b.n()) << std::endl;
253//  std::cout << (a.mean()==b.mean()) << std::endl;
254//  std::cout << (fabs(a.variance()-b.variance()<1e-15)) << std::endl;
255  return (a.n()==b.n() && a.mean()==b.mean() &&
256          fabs(a.variance()-b.variance()<1e-15));
257}
258
259bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
260           const double tol, std::ostream* error)
261{
262  bool equal = true;
263  if  ( fabs(a.mean()-b.mean())>tol){
264    equal=false;
265    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
266  }
267  if ( fabs(a.variance()-b.variance())>tol ) {
268    equal=false;
269    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
270           << std::endl;
271  } 
272  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
273    equal =false;
274    *error << "error for standard error:\t" << std::endl;
275  }
276  return equal;
277}
278
279bool equal(const Averager& a, const AveragerWeighted& b, const double tol, 
280           std::ostream* error)
281{
282  bool equal = true;
283  if  ( fabs(a.mean()-b.mean())>tol){
284    equal=false;
285    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
286  }
287  if ( fabs(a.variance()-b.variance())>tol ) {
288    equal=false;
289    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
290           << std::endl;
291  } 
292  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
293    equal =false;
294    *error << "error for standard error:\t" << std::endl;
295  }
296  return equal;
297}
298
299bool equal(const AveragerPair& a, const AveragerPair& b, 
300           const double tol, std::ostream* error)
301{
302  bool ok = true;
303  if  ( fabs(a.covariance()-b.covariance())>tol){
304    ok=false;
305    *error << "error covariance: " << a.covariance() << "\t" 
306           << b.covariance() << std::endl;
307  }
308  if ( fabs(a.correlation()-b.correlation())>tol ) {
309    ok=false;
310    *error << "error correlation" << std::endl;
311  } 
312  return ok;
313}
314
315bool equal(const AveragerPair& a, const AveragerPairWeighted& b, 
316           const double tol, std::ostream* error)
317{
318  bool ok = true;
319  if  ( fabs(a.covariance()-b.covariance())>tol){
320    ok=false;
321    *error << "error covariance: " << a.covariance() << "\t" 
322           << b.covariance() << std::endl;
323  }
324  if ( fabs(a.correlation()-b.correlation())>tol ) {
325    ok=false;
326    *error << "error correlation" << std::endl;
327  } 
328  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
329    ok =false;
330    *error << "error for x_averager():\t" << std::endl;
331  }
332  return ok;
333}
334bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
335           const double tol, std::ostream* error)
336{
337  bool ok = true;
338  if  ( fabs(a.covariance()-b.covariance())>tol){
339    ok=false;
340    *error << "error covariance: " << a.covariance() << "\t" 
341           << b.covariance() << std::endl;
342  }
343  if ( fabs(a.correlation()-b.correlation())>tol ) {
344    ok=false;
345    *error << "error correlation" << std::endl;
346  } 
347  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
348    ok =false;
349    *error << "error for x_averager():\t" << std::endl;
350  }
351  return ok;
352}
353
354
355
Note: See TracBrowser for help on using the repository browser.