source: branches/0.4-stable/test/averager_test.cc @ 1291

Last change on this file since 1291 was 1291, checked in by Jari Häkkinen, 15 years ago

Fixes #361
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.6 KB
Line 
1// $Id: averager_test.cc 1291 2008-05-09 19:42:31Z jari $
2
3/*
4  Copyright (C) 2005 Peter Johansson
5  Copyright (C) 2006 Jari Häkkinen, Peter Johansson, Markus Ringnér
6  Copyright (C) 2007, 2008 Jari Häkkinen, Peter Johansson
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 "Suite.h"
27
28#include "yat/statistics/Averager.h"
29#include "yat/statistics/AveragerPair.h"
30#include "yat/statistics/AveragerPairWeighted.h"
31#include "yat/statistics/AveragerWeighted.h"
32#include "yat/utility/Vector.h"
33
34#include <cmath>
35#include <fstream>
36#include <limits>
37#include <iostream>
38#include <set>
39
40using namespace theplu::yat;
41using namespace theplu::yat::statistics;
42
43using theplu::yat::test::Suite;
44
45//Forward declarations
46bool equal(const Averager&, const Averager&, unsigned int, Suite&);
47bool equal(const AveragerWeighted&, const AveragerWeighted&, unsigned int, 
48           Suite& suite);
49bool equal(const Averager&, const AveragerWeighted&, unsigned int, Suite& suite);
50bool equal(const AveragerPair&, const AveragerPair&, unsigned int, Suite& suite);
51bool equal(const AveragerPair&, const AveragerPairWeighted&, unsigned int, 
52           Suite& suite);
53bool equal(const AveragerPairWeighted&, const AveragerPairWeighted&,
54           unsigned int, Suite& suite);
55
56
57int main(int argc, char* argv[])
58{ 
59
60  Suite suite(argc, argv);
61
62  // Testing Averager
63  suite.out() << "testing Averager" << std::endl;
64  Averager a;
65  a.add(1);
66  a.add(3);
67  a.add(5);
68  if (a.n()!=3 || a.mean()!=3 || a.sum_xx()!=35){
69    suite.add(false);
70    suite.err() << "error: add\n";
71  }
72  Averager b(a);
73  b.add(5,-1);
74  if (b.n()!=2 || b.mean()!=2 || b.sum_xx()!=10){
75    suite.add(false);
76    suite.err() << "error: add with negative n\n";
77  }
78  b.add(5,-4);
79  if (b.n()!=-2 || b.mean()!=8 || b.sum_xx()!=-90){
80    suite.add(false);
81    suite.err() << "error: add with negatibe n\n";
82  }
83
84
85  Averager* a1 = new Averager(1.0+3+5,1.0+9+25,3);
86  unsigned int tol = 10;
87  if (!equal(a,*a1, tol, suite)){
88    suite.add(false);
89    suite.err().precision(25);
90    suite.err() << a.sum_x() << '\t' << a1->sum_x() << std::endl;
91    suite.err() << a.sum_xx() << '\t' << a1->sum_xx() << std::endl;
92    suite.err() << a.n() << '\t' << a1->n() << std::endl;
93    suite.err() << a.variance() << '\t' << a1->variance() << std::endl;
94    suite.err() << a.mean() << '\t' << a1->mean() << std::endl;
95    suite.err() << a.mean() - a1->mean() << std::endl;
96    suite.err() << a.variance() - a1->variance() << std::endl;
97    suite.err() << "error: Averager(double x, double xx, unsigned long n)\n";
98  }
99  delete a1;
100
101  a1 = new Averager(a);
102  if (!equal(a,*a1, tol, suite)){
103    suite.add(false);
104    suite.err() << "error: Copy constructor\n";
105  }
106  delete a1;
107
108  a.add(3,5);
109  if (std::abs(a.standard_error()-sqrt(a.variance()/a.n()))>
110      std::numeric_limits<double>().round_error() ){
111    suite.add(false);
112    suite.err() << "error: standard_error\n";
113  }
114     
115
116  if ( std::abs(a.variance() - a.std()*a.std())>
117       std::numeric_limits<double>().round_error() ){
118    suite.add(false);
119    suite.err() << "error: std squared should be variance" << std::endl;
120    suite.err() << "std2: " << a.std()*a.std() << std::endl;
121    suite.err() << "variance: " << a.variance() << std::endl;
122    suite.err() << "difference is: " << a.std()*a.std()-a.variance() << std::endl;
123  }
124 
125  if ( a.variance() != a.variance(a.mean()) ){
126    suite.add(false);
127    suite.err() << "error: variance incorrect\n" << std::endl;
128    suite.err() << "variance: " << a.variance() << std::endl;
129    suite.err() << "mean: " << a.mean() << std::endl;
130    suite.err() << "variance(mean) " << a.variance(a.mean()) << std::endl;
131  }
132  theplu::yat::utility::Vector* tmp_vec = new theplu::yat::utility::Vector(10);
133  add(a, tmp_vec->begin(), tmp_vec->end());
134  delete tmp_vec;
135  std::set<double>* tmp_set = new std::set<double>;
136  tmp_set->insert(1.0);
137  tmp_set->insert(2.0);
138  add(a, tmp_set->begin(), tmp_set->end());
139  delete tmp_set;
140
141
142  // Testing AveragerWeighted
143  suite.err() << "testing AveragerWeighted" << std::endl;
144  theplu::yat::utility::Vector x(3,0);
145  x(0)=0;
146  x(1)=1;
147  x(2)=2;
148  theplu::yat::utility::Vector w(3,1);
149  theplu::yat::statistics::AveragerWeighted aw;
150  add(aw, x.begin(), x.end(), w.begin());
151  a.reset();
152  add(a, x.begin(), x.end());
153  if (!equal(a,aw,tol,suite)){
154    suite.err() << "error: AveragerWeighted with unitary weights should " 
155           << "be equal to Averager" << std::endl;
156    suite.add(false);
157  }
158
159  AveragerWeighted* aw2 = new AveragerWeighted(aw);
160  if (!equal(aw,*aw2,tol,suite)){
161    suite.err() << "error: AveragerWeighted copy constructor " << std::endl;
162    suite.add(false);
163  }
164   
165  aw2->add(12,0);
166  if (!equal(aw,*aw2,tol,suite)){
167    suite.err() << "error: AveragerWeighted adding a data point with weight=0 " 
168           << "should make no change " << std::endl;
169    suite.add(false);
170  }
171     
172  aw2->reset();
173  w*=17;
174  add(*aw2, x.begin(), x.end(), w.begin());
175  if (!equal(aw,*aw2,tol,suite)){
176    suite.err() << "error: AveragerWeighted rescaling weights " 
177           << "should make no change " << std::endl;
178    suite.add(false);
179  }
180  delete aw2;
181  {
182    theplu::yat::utility::Vector tmp(10);
183    add(aw, tmp.begin(), tmp.end());
184  }
185  {
186    std::set<double> tmp;
187    tmp.insert(1.0);
188    tmp.insert(2.0);
189    add(aw, tmp.begin(), tmp.end());
190  }
191 
192
193  suite.out() << "testing AveragerPair" << std::endl;
194  AveragerPair ap;
195  for (int i=0; i<10; i++)
196    ap.add(static_cast<double>(i),i);
197  if (std::abs(ap.correlation()-1)>tol){
198    suite.add(false);
199    suite.err() << "correlation: " << ap.correlation() << std::endl;
200    suite.err() << "error: correlation between identical vectors should be unity" 
201           << std::endl;
202  }
203  if (ap.x_averager().variance()!=ap.covariance()){
204    suite.add(false);
205    suite.err() << "error: covariance of identical vectors should equal to variance" 
206           << std::endl;
207  }
208  AveragerPair* ap2 = new AveragerPair(ap);
209  delete ap2;
210
211  for (int i=0; i<8; i++)
212    ap.add(static_cast<double>(i),i,-1);
213  if (std::abs(ap.correlation()-1)>tol) {
214    suite.add(false);
215    suite.err() << "correlation after removal of data: " << ap.correlation()
216                << std::endl;
217    suite.err() << "error: correlation between identical vectors is unity" 
218                << std::endl;
219  }
220  ap.add(static_cast<double>(1),1,-4);
221  if (std::abs(ap.correlation()+1)>tol || ap.n()!=-2) {
222    suite.add(false);
223    suite.err() << "AveragerPair error: add with negatibe n\n";
224  }
225
226
227  suite.out() << "testing AveragerPairWeighted" << std::endl;
228  AveragerPairWeighted apw;
229  x(0)=0; x(1)=1; x(2)=2;
230  theplu::yat::utility::Vector y(3,0);
231  y(0)=0; y(1)=0; y(2)=2;
232  add(apw, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
233  ap.reset();
234  add(ap, x.begin(), x.end(), y.begin());
235  if (!equal(ap,apw,tol,suite)){
236    suite.err() << "error: AveragerPairWeighted with unitary weights should " 
237           << "be equal to AveragerPair" << std::endl;
238    suite.add(false);
239  }
240
241  AveragerPairWeighted* apw2 = new AveragerPairWeighted(apw);
242  if (!equal(apw,*apw2,tol,suite)){
243    suite.err() << "error: AveragerPairWeighted copy constructor " << std::endl;
244    suite.add(false);
245  }
246   
247  apw2->add(12,23222.03,32.3,0);
248  if (!equal(apw,*apw2,tol,suite)){
249    suite.err() << "error: AveragerWeighted adding a data point with weight=0 " 
250           << "should make no change " << std::endl;
251    suite.add(false);
252  }
253     
254  apw2->reset();
255  w*=17;
256  add(*apw2, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
257  if (!equal(apw,*apw2,tol,suite)){
258    suite.err() << "error: AveragerWeighted rescaling weights " 
259           << "should make no change " << std::endl;
260    suite.add(false);
261  }
262  delete apw2;
263
264  return suite.return_value();
265}
266
267bool equal(const Averager& a, const Averager& b, unsigned int tol, 
268           Suite& suite)
269{
270//  std::cout << (a.n()==b.n()) << std::endl;
271//  std::cout << (a.mean()==b.mean()) << std::endl;
272//  std::cout << (std::abs(a.variance()-b.variance()<1e-15)) << std::endl;
273  return (a.n()==b.n() && suite.equal(a.mean(),b.mean(),tol) &&
274          suite.equal(a.variance(),b.variance(),tol));
275}
276
277bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
278           const unsigned int tol, Suite& suite)
279{
280  bool equal = true;
281  if  ( !suite.equal(a.mean(),b.mean(),tol)){
282    equal=false;
283    suite.err() << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
284    suite.err() << "difference:\t" << a.mean()-b.mean() << std::endl;
285  }
286  if ( !suite.equal(a.variance(),b.variance(),tol) ) {
287    equal=false;
288    suite.err() << "error for variance:\t" << a.variance() << " " 
289                << b.variance() << std::endl;
290  } 
291  if ( !suite.equal(a.standard_error(),b.standard_error(),tol) ) {
292    equal =false;
293    suite.err() << "error for standard error:\t" << std::endl;
294  }
295  return equal;
296}
297
298bool equal(const Averager& a, const AveragerWeighted& b, const unsigned int tol, 
299           Suite& suite)
300{
301  bool equal = true;
302  if  ( !suite.equal(a.mean(),b.mean(),tol)){
303    equal=false;
304    suite.err() << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
305  }
306  if ( !suite.equal(a.variance(),b.variance(),tol) ) {
307    equal=false;
308    suite.err() << "error for variance:\t" << a.variance() << " " 
309                << b.variance() << std::endl;
310  } 
311  if ( !suite.equal(a.standard_error(), b.standard_error(),tol) ) {
312    equal =false;
313    suite.err() << "error for standard error:\t" << std::endl;
314  }
315  return equal;
316}
317
318bool equal(const AveragerPair& a, const AveragerPair& b, 
319           const unsigned int tol, Suite& suite)
320{
321  bool ok = true;
322  if  ( std::abs(a.covariance()-b.covariance())>tol){
323    suite.add(false);
324    suite.err() << "error covariance: " << a.covariance() << "\t" 
325           << b.covariance() << std::endl;
326  }
327  if ( std::abs(a.correlation()-b.correlation())>tol ) {
328    suite.add(false);
329    suite.err() << "error correlation" << std::endl;
330  } 
331  return ok;
332}
333
334bool equal(const AveragerPair& a, const AveragerPairWeighted& b, 
335           const unsigned int tol, Suite& suite)
336{
337  bool ok = true;
338  if  ( std::abs(a.covariance()-b.covariance())>tol){
339    suite.add(false);
340    suite.err() << "error covariance: " << a.covariance() << "\t" 
341           << b.covariance() << std::endl;
342  }
343  if ( std::abs(a.correlation()-b.correlation())>tol ) {
344    suite.add(false);
345    suite.err() << "error correlation" << std::endl;
346    suite.err() << "unweighted:" << a.correlation() << std::endl;
347    suite.err() << "weighted:" << b.correlation() << std::endl;
348    suite.err() << "difference:" << a.correlation()-b.correlation() 
349                << std::endl;
350  } 
351  if ( !equal(a.x_averager(),b.x_averager(),tol,suite)) {
352    ok =false;
353    suite.err() << "error for x_averager():\t" << std::endl;
354  }
355  return ok;
356}
357bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
358           const unsigned int tol, Suite& suite)
359{
360  bool ok = true;
361  if  ( !suite.equal(a.covariance(),b.covariance(),tol) ){
362    suite.add(false);
363    suite.err() << "error covariance: " << a.covariance() << "\t" 
364           << b.covariance() << std::endl;
365  }
366  if ( !suite.equal(a.correlation(), b.correlation(), tol) ) {
367    suite.add(false);
368    suite.err() << "error correlation" << std::endl;
369    suite.err() << "a:" << a.correlation() << std::endl;
370    suite.err() << "b:" << b.correlation() << std::endl;
371    suite.err() << "difference:" << a.correlation()-b.correlation() 
372                << std::endl;
373    suite.err() << "tol:" << tol << std::endl;
374  } 
375  if ( !equal(a.x_averager(),b.x_averager(),tol,suite)) {
376    ok =false;
377    suite.err() << "error for x_averager():\t" << std::endl;
378  }
379  return ok;
380}
381
382
383
Note: See TracBrowser for help on using the repository browser.