source: trunk/test/averager_test.cc @ 617

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

Cleaned up stray gslapi stuff.

  • 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 617 2006-08-31 08:58:05Z jari $
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  Averager* a1 = new Averager(1+3+5,1+9+25,3);
54  if (!equal(a,*a1)){
55    ok=false;
56    std::cout << a.sum_x() << std::endl;
57    std::cout << a.sum_xx() << std::endl;
58    std::cout << a.n() << std::endl;
59    std::cout << a.variance() << std::endl;
60    std::cout << a.mean() << std::endl;
61    std::cout << a1->sum_x() << std::endl;
62    std::cout << a1->sum_xx() << std::endl;
63    std::cout << a1->n() << std::endl;
64    std::cout << a1->variance() << std::endl;
65    std::cout << a1->mean() << std::endl;
66    std::cout << a.mean() - a1->mean() << std::endl;
67    std::cout << a.variance() - a1->variance() << std::endl;
68    *error << "error: Averager(const double x,const double xx,const long n)\n";
69  }
70  delete a1;
71
72  a1 = new Averager(a);
73  if (!equal(a,*a1)){
74    ok=false;
75    *error << "error: Copy constructor\n";
76  }
77  delete a1;
78
79  a.add(3,5);
80  if (! a.standard_error()==sqrt(a.variance()/a.n())){
81    ok=false;
82    *error << "error: standard_error\n";
83  }
84     
85
86  if ( fabs(a.variance() - a.std()*a.std())>
87       std::numeric_limits<double>().round_error() ){
88    ok=false;
89    *error << "error: std squared should be variance" << std::endl;
90    *error << "std2: " << a.std()*a.std() << std::endl;
91    *error << "variance: " << a.variance() << std::endl;
92    *error << "difference is: " << a.std()*a.std()-a.variance() << std::endl;
93  }
94 
95  if ( a.variance() != a.variance(a.mean()) ){
96    ok=false;
97    *error << "error: variance incorrect\n" << std::endl;
98    *error << "variance: " << a.variance() << std::endl;
99    *error << "mean: " << a.mean() << std::endl;
100    *error << "variance(mean) " << a.variance(a.mean()) << std::endl;
101  }
102
103  // Testing AveragerWeighted
104  *error << "testing AveragerWeighted" << std::endl;
105  theplu::utility::vector x(3,0);
106  x(0)=0;
107  x(1)=1;
108  x(2)=2;
109  theplu::utility::vector w(3,1);
110  theplu::statistics::AveragerWeighted aw;
111  aw.add_values(x,w);
112  a.reset();
113  a.add_values(x);
114  const double tol=std::numeric_limits<double>().round_error();
115  if (!equal(a,aw,tol,error)){
116    *error << "error: AveragerWeighted with unitary weights should " 
117           << "be equal to Averager" << std::endl;
118    ok=false;
119  }
120
121  AveragerWeighted* aw2 = new AveragerWeighted(aw);
122  if (!equal(aw,*aw2,tol,error)){
123    *error << "error: AveragerWeighted copy constructor " << std::endl;
124    ok=false;
125  }
126   
127  aw2->add(12,0);
128  if (!equal(aw,*aw2,tol,error)){
129    *error << "error: AveragerWeighted adding a data point with weight=0 " 
130           << "should make no change " << std::endl;
131    ok=false;
132  }
133     
134  aw2->reset();
135  w.scale(17);
136  aw2->add_values(x,w);
137  if (!equal(aw,*aw2,tol,error)){
138    *error << "error: AveragerWeighted rescaling weights " 
139           << "should make no change " << std::endl;
140    ok=false;
141  }
142  delete aw2;
143 
144
145  *error << "testing AveragerPair" << std::endl;
146  AveragerPair ap;
147  for (int i=0; i<10; i++)
148    ap.add(static_cast<double>(i),i);
149  if (fabs(ap.correlation()-1)>tol){
150    ok=false;
151    *error << "correlation: " << ap.correlation() << std::endl;
152    *error << "error: correlation between identical vectors should be unity" 
153           << std::endl;
154  }
155  if (ap.x_averager().variance()!=ap.covariance()){
156    ok=false;
157    *error << "error: covariance of identical vectors should equal to variance" 
158           << std::endl;
159  }
160  AveragerPair* ap2 = new AveragerPair(ap);
161  delete ap2;
162
163  *error << "testing AveragerPairWeighted" << std::endl;
164  AveragerPairWeighted apw;
165  x(0)=0; x(1)=1; x(2)=2;
166  theplu::utility::vector y(3,0);
167  x(0)=0; x(1)=0; x(2)=2;
168  apw.add_values(x,y,w,w);
169  ap.reset();
170  ap.add_values(x,y);
171  if (!equal(ap,apw,tol,error)){
172    *error << "error: AveragerPairWeighted with unitary weights should " 
173           << "be equal to AveragerPair" << std::endl;
174    ok=false;
175  }
176
177  AveragerPairWeighted* apw2 = new AveragerPairWeighted(apw);
178  if (!equal(apw,*apw2,tol,error)){
179    *error << "error: AveragerPairWeighted copy constructor " << std::endl;
180    ok=false;
181  }
182   
183  apw2->add(12,23222.03,32.3,0);
184  if (!equal(apw,*apw2,tol,error)){
185    *error << "error: AveragerWeighted adding a data point with weight=0 " 
186           << "should make no change " << std::endl;
187    ok=false;
188  }
189     
190  apw2->reset();
191  w.scale(17);
192  apw2->add_values(x,y,w,w);
193  if (!equal(apw,*apw2,tol,error)){
194    *error << "error: AveragerWeighted rescaling weights " 
195           << "should make no change " << std::endl;
196    ok=false;
197  }
198  delete apw2;
199
200  if (error!=&std::cerr)
201    delete error;
202
203  if (!ok)
204    return -1;
205  return 0;
206}
207
208bool equal(const Averager& a, const Averager& b)
209{
210//  std::cout << (a.n()==b.n()) << std::endl;
211//  std::cout << (a.mean()==b.mean()) << std::endl;
212//  std::cout << (a.variance()==b.variance()) << std::endl;
213  return (a.n()==b.n() && a.mean()==b.mean() && a.variance()==b.variance());
214}
215
216bool equal(const AveragerWeighted& a, const AveragerWeighted& b, 
217           const double tol, std::ostream* error)
218{
219  bool equal = true;
220  if  ( fabs(a.mean()-b.mean())>tol){
221    equal=false;
222    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
223  }
224  if ( fabs(a.variance()-b.variance())>tol ) {
225    equal=false;
226    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
227           << std::endl;
228  } 
229  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
230    equal =false;
231    *error << "error for standard error:\t" << std::endl;
232  }
233  return equal;
234}
235
236bool equal(const Averager& a, const AveragerWeighted& b, const double tol, 
237           std::ostream* error)
238{
239  bool equal = true;
240  if  ( fabs(a.mean()-b.mean())>tol){
241    equal=false;
242    *error << "mean:\t" << a.mean() << "\t" << b.mean() << std::endl;
243  }
244  if ( fabs(a.variance()-b.variance())>tol ) {
245    equal=false;
246    *error << "error for variance:\t" << a.variance() << " " << b.variance() 
247           << std::endl;
248  } 
249  if ( fabs(a.standard_error()-b.standard_error())>tol ) {
250    equal =false;
251    *error << "error for standard error:\t" << std::endl;
252  }
253  return equal;
254}
255
256bool equal(const AveragerPair& a, const AveragerPair& b, 
257           const double tol, std::ostream* error)
258{
259  bool ok = true;
260  if  ( fabs(a.covariance()-b.covariance())>tol){
261    ok=false;
262    *error << "error covariance: " << a.covariance() << "\t" 
263           << b.covariance() << std::endl;
264  }
265  if ( fabs(a.correlation()-b.correlation())>tol ) {
266    ok=false;
267    *error << "error correlation" << std::endl;
268  } 
269  return ok;
270}
271
272bool equal(const AveragerPair& a, const AveragerPairWeighted& b, 
273           const double tol, std::ostream* error)
274{
275  bool ok = true;
276  if  ( fabs(a.covariance()-b.covariance())>tol){
277    ok=false;
278    *error << "error covariance: " << a.covariance() << "\t" 
279           << b.covariance() << std::endl;
280  }
281  if ( fabs(a.correlation()-b.correlation())>tol ) {
282    ok=false;
283    *error << "error correlation" << std::endl;
284  } 
285  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
286    ok =false;
287    *error << "error for x_averager():\t" << std::endl;
288  }
289  return ok;
290}
291bool equal(const AveragerPairWeighted& a, const AveragerPairWeighted& b, 
292           const double tol, std::ostream* error)
293{
294  bool ok = true;
295  if  ( fabs(a.covariance()-b.covariance())>tol){
296    ok=false;
297    *error << "error covariance: " << a.covariance() << "\t" 
298           << b.covariance() << std::endl;
299  }
300  if ( fabs(a.correlation()-b.correlation())>tol ) {
301    ok=false;
302    *error << "error correlation" << std::endl;
303  } 
304  if ( !equal(a.x_averager(),b.x_averager(),tol,error)) {
305    ok =false;
306    *error << "error for x_averager():\t" << std::endl;
307  }
308  return ok;
309}
310
311
312
Note: See TracBrowser for help on using the repository browser.