source: trunk/test/averager_test.cc @ 1619

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

Addresses #436. GPL license copy reference should also be updated.

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