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

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

Addresses #361. Averager accepts negative n for removal of data from the Averager.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.1 KB
Line 
1// $Id: averager_test.cc 1290 2008-05-09 11:44:25Z 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  suite.out() << "testing AveragerPairWeighted" << std::endl;
212  AveragerPairWeighted apw;
213  x(0)=0; x(1)=1; x(2)=2;
214  theplu::yat::utility::Vector y(3,0);
215  y(0)=0; y(1)=0; y(2)=2;
216  add(apw, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
217  ap.reset();
218  add(ap, x.begin(), x.end(), y.begin());
219  if (!equal(ap,apw,tol,suite)){
220    suite.err() << "error: AveragerPairWeighted with unitary weights should " 
221           << "be equal to AveragerPair" << std::endl;
222    suite.add(false);
223  }
224
225  AveragerPairWeighted* apw2 = new AveragerPairWeighted(apw);
226  if (!equal(apw,*apw2,tol,suite)){
227    suite.err() << "error: AveragerPairWeighted copy constructor " << std::endl;
228    suite.add(false);
229  }
230   
231  apw2->add(12,23222.03,32.3,0);
232  if (!equal(apw,*apw2,tol,suite)){
233    suite.err() << "error: AveragerWeighted adding a data point with weight=0 " 
234           << "should make no change " << std::endl;
235    suite.add(false);
236  }
237     
238  apw2->reset();
239  w*=17;
240  add(*apw2, x.begin(), x.end(), y.begin(), w.begin(), w.begin());
241  if (!equal(apw,*apw2,tol,suite)){
242    suite.err() << "error: AveragerWeighted rescaling weights " 
243           << "should make no change " << std::endl;
244    suite.add(false);
245  }
246  delete apw2;
247
248  return suite.return_value();
249}
250
251bool equal(const Averager& a, const Averager& b, unsigned int tol, 
252           Suite& suite)
253{
254//  std::cout << (a.n()==b.n()) << std::endl;
255//  std::cout << (a.mean()==b.mean()) << std::endl;
256//  std::cout << (std::abs(a.variance()-b.variance()<1e-15)) << std::endl;
257  return (a.n()==b.n() && suite.equal(a.mean(),b.mean(),tol) &&
258          suite.equal(a.variance(),b.variance(),tol));
259}
260
261bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
262           const unsigned int tol, Suite& suite)
263{
264  bool equal = true;
265  if  ( !suite.equal(a.mean(),b.mean(),tol)){
266    equal=false;
267    suite.err() << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
268    suite.err() << "difference:\t" << a.mean()-b.mean() << std::endl;
269  }
270  if ( !suite.equal(a.variance(),b.variance(),tol) ) {
271    equal=false;
272    suite.err() << "error for variance:\t" << a.variance() << " " 
273                << b.variance() << std::endl;
274  } 
275  if ( !suite.equal(a.standard_error(),b.standard_error(),tol) ) {
276    equal =false;
277    suite.err() << "error for standard error:\t" << std::endl;
278  }
279  return equal;
280}
281
282bool equal(const Averager& a, const AveragerWeighted& b, const unsigned int tol, 
283           Suite& suite)
284{
285  bool equal = true;
286  if  ( !suite.equal(a.mean(),b.mean(),tol)){
287    equal=false;
288    suite.err() << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
289  }
290  if ( !suite.equal(a.variance(),b.variance(),tol) ) {
291    equal=false;
292    suite.err() << "error for variance:\t" << a.variance() << " " 
293                << b.variance() << std::endl;
294  } 
295  if ( !suite.equal(a.standard_error(), b.standard_error(),tol) ) {
296    equal =false;
297    suite.err() << "error for standard error:\t" << std::endl;
298  }
299  return equal;
300}
301
302bool equal(const AveragerPair& a, const AveragerPair& b, 
303           const unsigned int tol, Suite& suite)
304{
305  bool ok = true;
306  if  ( std::abs(a.covariance()-b.covariance())>tol){
307    suite.add(false);
308    suite.err() << "error covariance: " << a.covariance() << "\t" 
309           << b.covariance() << std::endl;
310  }
311  if ( std::abs(a.correlation()-b.correlation())>tol ) {
312    suite.add(false);
313    suite.err() << "error correlation" << std::endl;
314  } 
315  return ok;
316}
317
318bool equal(const AveragerPair& a, const AveragerPairWeighted& 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    suite.err() << "unweighted:" << a.correlation() << std::endl;
331    suite.err() << "weighted:" << b.correlation() << std::endl;
332    suite.err() << "difference:" << a.correlation()-b.correlation() 
333                << std::endl;
334  } 
335  if ( !equal(a.x_averager(),b.x_averager(),tol,suite)) {
336    ok =false;
337    suite.err() << "error for x_averager():\t" << std::endl;
338  }
339  return ok;
340}
341bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
342           const unsigned int tol, Suite& suite)
343{
344  bool ok = true;
345  if  ( !suite.equal(a.covariance(),b.covariance(),tol) ){
346    suite.add(false);
347    suite.err() << "error covariance: " << a.covariance() << "\t" 
348           << b.covariance() << std::endl;
349  }
350  if ( !suite.equal(a.correlation(), b.correlation(), tol) ) {
351    suite.add(false);
352    suite.err() << "error correlation" << std::endl;
353    suite.err() << "a:" << a.correlation() << std::endl;
354    suite.err() << "b:" << b.correlation() << std::endl;
355    suite.err() << "difference:" << a.correlation()-b.correlation() 
356                << std::endl;
357    suite.err() << "tol:" << tol << std::endl;
358  } 
359  if ( !equal(a.x_averager(),b.x_averager(),tol,suite)) {
360    ok =false;
361    suite.err() << "error for x_averager():\t" << std::endl;
362  }
363  return ok;
364}
365
366
367
Note: See TracBrowser for help on using the repository browser.