source: trunk/test/averager_test.cc @ 620

Last change on this file since 620 was 620, checked in by Markus Ringnér, 16 years ago

Fixed problems with Averager that sometimes appeared as failures in averager_test when optimization is turned on

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 8.4 KB
Line 
1// $Id: averager_test.cc 620 2006-09-04 15:16:57Z markus $
2
3#include <c++_tools/statistics/Averager.h>
4#include <c++_tools/statistics/AveragerPair.h>
5#include <c++_tools/statistics/AveragerPairWeighted.h>
6#include <c++_tools/statistics/AveragerWeighted.h>
7#include <c++_tools/utility/vector.h>
8
9#include <fstream>
10#include <limits>
11#include <iostream>
12
13using namespace theplu::statistics;
14
15//Forward declarations
16bool equal(const Averager&, const Averager&);
17bool equal(const AveragerWeighted&, const AveragerWeighted&, const double, 
18           std::ostream* error);
19bool equal(const Averager&, const AveragerWeighted&, const double, 
20           std::ostream* error);
21bool equal(const AveragerPair&, const AveragerPair&, 
22           const double, std::ostream* error);
23bool equal(const AveragerPair&, const AveragerPairWeighted&, 
24           const double, std::ostream* error);
25bool equal(const AveragerPairWeighted&, const AveragerPairWeighted&, 
26           const double, std::ostream* error);
27
28
29int main(const int argc,const char* argv[])
30{ 
31
32  std::ostream* error;
33  if (argc>1 && argv[1]==std::string("-v"))
34    error = &std::cerr;
35  else {
36    error = new std::ofstream("/dev/null");
37    if (argc>1)
38      std::cout << "averager_test -v : for printing extra information\n";
39  }
40  bool ok = true;
41
42  // Testing Averager
43  *error << "testing Averager" << std::endl;
44  Averager a;
45  a.add(1);
46  a.add(3);
47  a.add(5);
48  if (a.n()!=3 || a.mean()!=3 || a.sum_xx()!=35){
49    ok=false;
50    *error << "error: add\n";
51  }
52
53
54  Averager* a1 = new Averager(1.0+3+5,1.0+9+25,3);
55  if (!equal(a,*a1)){
56    ok=false;
57    std::cout.precision(25);
58    std::cout << equal(a,*a1) << std::endl;
59    std::cout << a.sum_x() << '\t' << a1->sum_x() << std::endl;
60    std::cout << a.sum_xx() << '\t' << a1->sum_xx() << std::endl;
61    std::cout << a.n() << '\t' << a1->n() << std::endl;
62    std::cout << a.variance() << '\t' << a1->variance() << std::endl;
63    std::cout << a.mean() << '\t' << a1->mean() << std::endl;
64    std::cout << a.mean() - a1->mean() << std::endl;
65    std::cout << a.variance() - a1->variance() << std::endl;
66    *error << "error: Averager(double x, double xx, u_long n)\n";
67  }
68  delete a1;
69
70  a1 = new Averager(a);
71  if (!equal(a,*a1)){
72    ok=false;
73    *error << "error: Copy constructor\n";
74  }
75  delete a1;
76
77  a.add(3,5);
78  if (! a.standard_error()==sqrt(a.variance()/a.n())){
79    ok=false;
80    *error << "error: standard_error\n";
81  }
82     
83
84  if ( fabs(a.variance() - a.std()*a.std())>
85       std::numeric_limits<double>().round_error() ){
86    ok=false;
87    *error << "error: std squared should be variance" << std::endl;
88    *error << "std2: " << a.std()*a.std() << std::endl;
89    *error << "variance: " << a.variance() << std::endl;
90    *error << "difference is: " << a.std()*a.std()-a.variance() << std::endl;
91  }
92 
93  if ( a.variance() != a.variance(a.mean()) ){
94    ok=false;
95    *error << "error: variance incorrect\n" << std::endl;
96    *error << "variance: " << a.variance() << std::endl;
97    *error << "mean: " << a.mean() << std::endl;
98    *error << "variance(mean) " << a.variance(a.mean()) << std::endl;
99  }
100
101  // Testing AveragerWeighted
102  *error << "testing AveragerWeighted" << std::endl;
103  theplu::utility::vector x(3,0);
104  x(0)=0;
105  x(1)=1;
106  x(2)=2;
107  theplu::utility::vector w(3,1);
108  theplu::statistics::AveragerWeighted aw;
109  aw.add_values(x,w);
110  a.reset();
111  a.add_values(x);
112  const double tol=std::numeric_limits<double>().round_error();
113  if (!equal(a,aw,tol,error)){
114    *error << "error: AveragerWeighted with unitary weights should " 
115           << "be equal to Averager" << std::endl;
116    ok=false;
117  }
118
119  AveragerWeighted* aw2 = new AveragerWeighted(aw);
120  if (!equal(aw,*aw2,tol,error)){
121    *error << "error: AveragerWeighted copy constructor " << std::endl;
122    ok=false;
123  }
124   
125  aw2->add(12,0);
126  if (!equal(aw,*aw2,tol,error)){
127    *error << "error: AveragerWeighted adding a data point with weight=0 " 
128           << "should make no change " << std::endl;
129    ok=false;
130  }
131     
132  aw2->reset();
133  w.scale(17);
134  aw2->add_values(x,w);
135  if (!equal(aw,*aw2,tol,error)){
136    *error << "error: AveragerWeighted rescaling weights " 
137           << "should make no change " << std::endl;
138    ok=false;
139  }
140  delete aw2;
141 
142
143  *error << "testing AveragerPair" << std::endl;
144  AveragerPair ap;
145  for (int i=0; i<10; i++)
146    ap.add(static_cast<double>(i),i);
147  if (fabs(ap.correlation()-1)>tol){
148    ok=false;
149    *error << "correlation: " << ap.correlation() << std::endl;
150    *error << "error: correlation between identical vectors should be unity" 
151           << std::endl;
152  }
153  if (ap.x_averager().variance()!=ap.covariance()){
154    ok=false;
155    *error << "error: covariance of identical vectors should equal to variance" 
156           << std::endl;
157  }
158  AveragerPair* ap2 = new AveragerPair(ap);
159  delete ap2;
160
161  *error << "testing AveragerPairWeighted" << std::endl;
162  AveragerPairWeighted apw;
163  x(0)=0; x(1)=1; x(2)=2;
164  theplu::utility::vector y(3,0);
165  x(0)=0; x(1)=0; x(2)=2;
166  apw.add_values(x,y,w,w);
167  ap.reset();
168  ap.add_values(x,y);
169  if (!equal(ap,apw,tol,error)){
170    *error << "error: AveragerPairWeighted with unitary weights should " 
171           << "be equal to AveragerPair" << std::endl;
172    ok=false;
173  }
174
175  AveragerPairWeighted* apw2 = new AveragerPairWeighted(apw);
176  if (!equal(apw,*apw2,tol,error)){
177    *error << "error: AveragerPairWeighted copy constructor " << std::endl;
178    ok=false;
179  }
180   
181  apw2->add(12,23222.03,32.3,0);
182  if (!equal(apw,*apw2,tol,error)){
183    *error << "error: AveragerWeighted adding a data point with weight=0 " 
184           << "should make no change " << std::endl;
185    ok=false;
186  }
187     
188  apw2->reset();
189  w.scale(17);
190  apw2->add_values(x,y,w,w);
191  if (!equal(apw,*apw2,tol,error)){
192    *error << "error: AveragerWeighted rescaling weights " 
193           << "should make no change " << std::endl;
194    ok=false;
195  }
196  delete apw2;
197
198  if (error!=&std::cerr)
199    delete error;
200
201  if (!ok)
202    return -1;
203  return 0;
204}
205
206bool equal(const Averager& a, const Averager& b)
207{
208//  std::cout << (a.n()==b.n()) << std::endl;
209//  std::cout << (a.mean()==b.mean()) << std::endl;
210//  std::cout << (a.variance()==b.variance()) << std::endl;
211  return (a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance());
212}
213
214bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
215           const double tol, std::ostream* error)
216{
217  bool equal = true;
218  if  ( fabs(a.mean()-b.mean())>tol){
219    equal=false;
220    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
221  }
222  if ( fabs(a.variance()-b.variance())>tol ) {
223    equal=false;
224    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
225           << std::endl;
226  } 
227  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
228    equal =false;
229    *error << "error for standard error:\t" << std::endl;
230  }
231  return equal;
232}
233
234bool equal(const Averager& a, const AveragerWeighted& b, const double tol, 
235           std::ostream* error)
236{
237  bool equal = true;
238  if  ( fabs(a.mean()-b.mean())>tol){
239    equal=false;
240    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
241  }
242  if ( fabs(a.variance()-b.variance())>tol ) {
243    equal=false;
244    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
245           << std::endl;
246  } 
247  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
248    equal =false;
249    *error << "error for standard error:\t" << std::endl;
250  }
251  return equal;
252}
253
254bool equal(const AveragerPair& a, const AveragerPair& b, 
255           const double tol, std::ostream* error)
256{
257  bool ok = true;
258  if  ( fabs(a.covariance()-b.covariance())>tol){
259    ok=false;
260    *error << "error covariance: " << a.covariance() << "\t" 
261           << b.covariance() << std::endl;
262  }
263  if ( fabs(a.correlation()-b.correlation())>tol ) {
264    ok=false;
265    *error << "error correlation" << std::endl;
266  } 
267  return ok;
268}
269
270bool equal(const AveragerPair& a, const AveragerPairWeighted& b, 
271           const double tol, std::ostream* error)
272{
273  bool ok = true;
274  if  ( fabs(a.covariance()-b.covariance())>tol){
275    ok=false;
276    *error << "error covariance: " << a.covariance() << "\t" 
277           << b.covariance() << std::endl;
278  }
279  if ( fabs(a.correlation()-b.correlation())>tol ) {
280    ok=false;
281    *error << "error correlation" << std::endl;
282  } 
283  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
284    ok =false;
285    *error << "error for x_averager():\t" << std::endl;
286  }
287  return ok;
288}
289bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
290           const double tol, std::ostream* error)
291{
292  bool ok = true;
293  if  ( fabs(a.covariance()-b.covariance())>tol){
294    ok=false;
295    *error << "error covariance: " << a.covariance() << "\t" 
296           << b.covariance() << std::endl;
297  }
298  if ( fabs(a.correlation()-b.correlation())>tol ) {
299    ok=false;
300    *error << "error correlation" << std::endl;
301  } 
302  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
303    ok =false;
304    *error << "error for x_averager():\t" << std::endl;
305  }
306  return ok;
307}
308
309
310
Note: See TracBrowser for help on using the repository browser.