1 |
/* |
2 |
* Copyright (c) 2012 The University of Notre Dame. All Rights Reserved. |
3 |
* |
4 |
* The University of Notre Dame grants you ("Licensee") a |
5 |
* non-exclusive, royalty free, license to use, modify and |
6 |
* redistribute this software in source and binary code form, provided |
7 |
* that the following conditions are met: |
8 |
* |
9 |
* 1. Redistributions of source code must retain the above copyright |
10 |
* notice, this list of conditions and the following disclaimer. |
11 |
* |
12 |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
* notice, this list of conditions and the following disclaimer in the |
14 |
* documentation and/or other materials provided with the |
15 |
* distribution. |
16 |
* |
17 |
* This software is provided "AS IS," without a warranty of any |
18 |
* kind. All express or implied conditions, representations and |
19 |
* warranties, including any implied warranty of merchantability, |
20 |
* fitness for a particular purpose or non-infringement, are hereby |
21 |
* excluded. The University of Notre Dame and its licensors shall not |
22 |
* be liable for any damages suffered by licensee as a result of |
23 |
* using, modifying or distributing the software or its |
24 |
* derivatives. In no event will the University of Notre Dame or its |
25 |
* licensors be liable for any lost revenue, profit or data, or for |
26 |
* direct, indirect, special, consequential, incidental or punitive |
27 |
* damages, however caused and regardless of the theory of liability, |
28 |
* arising out of the use of or inability to use software, even if the |
29 |
* University of Notre Dame has been advised of the possibility of |
30 |
* such damages. |
31 |
* |
32 |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
* research, please cite the appropriate papers when you publish your |
34 |
* work. Good starting points are: |
35 |
* |
36 |
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 |
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). |
39 |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
*/ |
42 |
|
43 |
#ifndef UTILS_ACCUMULATOR_HPP |
44 |
#define UTILS_ACCUMULATOR_HPP |
45 |
|
46 |
#include <cmath> |
47 |
#include <cassert> |
48 |
#include "math/Vector3.hpp" |
49 |
|
50 |
namespace OpenMD { |
51 |
|
52 |
|
53 |
class BaseAccumulator { |
54 |
public: |
55 |
virtual void clear() = 0; |
56 |
/** |
57 |
* get the number of accumulated values |
58 |
*/ |
59 |
virtual size_t count() { |
60 |
return Count_; |
61 |
} |
62 |
protected: |
63 |
size_t Count_; |
64 |
|
65 |
}; |
66 |
|
67 |
|
68 |
|
69 |
/** |
70 |
* Basic Accumulator class for numbers. |
71 |
*/ |
72 |
class Accumulator : public BaseAccumulator { |
73 |
|
74 |
typedef RealType ElementType; |
75 |
typedef RealType ResultType; |
76 |
|
77 |
public: |
78 |
|
79 |
Accumulator() : BaseAccumulator() { |
80 |
this->clear(); |
81 |
} |
82 |
|
83 |
/** |
84 |
* Accumulate another value |
85 |
*/ |
86 |
virtual void add(ElementType const& val) { |
87 |
Count_++; |
88 |
Avg_ += (val - Avg_ ) / Count_; |
89 |
Avg2_ += (val * val - Avg2_) / Count_; |
90 |
Val_ = val; |
91 |
if (Count_ <= 1) { |
92 |
Max_ = val; |
93 |
Min_ = val; |
94 |
} else { |
95 |
Max_ = val > Max_ ? val : Max_; |
96 |
Min_ = val < Min_ ? val : Min_; |
97 |
} |
98 |
} |
99 |
|
100 |
/** |
101 |
* reset the Accumulator to the empty state |
102 |
*/ |
103 |
void clear() { |
104 |
Count_ = 0; |
105 |
Avg_ = 0; |
106 |
Avg2_ = 0; |
107 |
Val_ = 0; |
108 |
} |
109 |
|
110 |
|
111 |
/** |
112 |
* return the most recently added value |
113 |
*/ |
114 |
void getLastValue(ElementType &ret) { |
115 |
ret = Val_; |
116 |
return; |
117 |
} |
118 |
|
119 |
/** |
120 |
* compute the Mean |
121 |
*/ |
122 |
void getAverage(ResultType &ret) { |
123 |
assert(Count_ != 0); |
124 |
ret = Avg_; |
125 |
return; |
126 |
} |
127 |
|
128 |
/** |
129 |
* compute the Variance |
130 |
*/ |
131 |
void getVariance(ResultType &ret) { |
132 |
assert(Count_ != 0); |
133 |
ret = (Avg2_ - Avg_ * Avg_); |
134 |
return; |
135 |
} |
136 |
|
137 |
/** |
138 |
* compute error of average value |
139 |
*/ |
140 |
void getStdDev(ResultType &ret) { |
141 |
assert(Count_ != 0); |
142 |
RealType var; |
143 |
this->getVariance(var); |
144 |
ret = sqrt(var); |
145 |
return; |
146 |
} |
147 |
|
148 |
/** |
149 |
* return the largest value |
150 |
*/ |
151 |
void getMax(ElementType &ret) { |
152 |
assert(Count_ != 0); |
153 |
ret = Max_; |
154 |
return; |
155 |
} |
156 |
|
157 |
/** |
158 |
* return the smallest value |
159 |
*/ |
160 |
void getMin(ElementType &ret) { |
161 |
assert(Count_ != 0); |
162 |
ret = Max_; |
163 |
return; |
164 |
} |
165 |
|
166 |
/** |
167 |
* return the 95% confidence interval: |
168 |
* |
169 |
* That is returns c, such that we have 95% confidence that the |
170 |
* true mean is within 2c of the Average (x): |
171 |
* |
172 |
* x - c <= true mean <= x + c |
173 |
* |
174 |
*/ |
175 |
void get95percentConfidenceInterval(ResultType &ret) { |
176 |
assert(Count_ != 0); |
177 |
RealType sd; |
178 |
this->getStdDev(sd); |
179 |
ret = 1.960 * sd / sqrt(Count_); |
180 |
return; |
181 |
} |
182 |
|
183 |
private: |
184 |
ElementType Val_; |
185 |
ResultType Avg_; |
186 |
ResultType Avg2_; |
187 |
ElementType Min_; |
188 |
ElementType Max_; |
189 |
}; |
190 |
|
191 |
class VectorAccumulator : public BaseAccumulator { |
192 |
|
193 |
typedef Vector3d ElementType; |
194 |
typedef Vector3d ResultType; |
195 |
|
196 |
public: |
197 |
VectorAccumulator() : BaseAccumulator() { |
198 |
this->clear(); |
199 |
} |
200 |
|
201 |
/** |
202 |
* Accumulate another value |
203 |
*/ |
204 |
void add(ElementType const& val) { |
205 |
Count_++; |
206 |
RealType len(0.0); |
207 |
for (unsigned int i =0; i < 3; i++) { |
208 |
Avg_[i] += (val[i] - Avg_[i] ) / Count_; |
209 |
Avg2_[i] += (val[i] * val[i] - Avg2_[i]) / Count_; |
210 |
Val_[i] = val[i]; |
211 |
len += val[i]*val[i]; |
212 |
} |
213 |
len = sqrt(len); |
214 |
AvgLen_ += (len - AvgLen_ ) / Count_; |
215 |
AvgLen2_ += (len * len - AvgLen2_) / Count_; |
216 |
|
217 |
if (Count_ <= 1) { |
218 |
Max_ = len; |
219 |
Min_ = len; |
220 |
} else { |
221 |
Max_ = len > Max_ ? len : Max_; |
222 |
Min_ = len < Min_ ? len : Min_; |
223 |
} |
224 |
} |
225 |
|
226 |
/** |
227 |
* reset the Accumulator to the empty state |
228 |
*/ |
229 |
void clear() { |
230 |
Count_ = 0; |
231 |
Avg_ = V3Zero; |
232 |
Avg2_ = V3Zero; |
233 |
Val_ = V3Zero; |
234 |
AvgLen_ = 0; |
235 |
AvgLen2_ = 0; |
236 |
} |
237 |
|
238 |
/** |
239 |
* return the most recently added value |
240 |
*/ |
241 |
void getLastValue(ElementType &ret) { |
242 |
ret = Val_; |
243 |
return; |
244 |
} |
245 |
|
246 |
/** |
247 |
* compute the Mean |
248 |
*/ |
249 |
void getAverage(ResultType &ret) { |
250 |
assert(Count_ != 0); |
251 |
ret = Avg_; |
252 |
return; |
253 |
} |
254 |
|
255 |
/** |
256 |
* compute the Variance |
257 |
*/ |
258 |
void getVariance(ResultType &ret) { |
259 |
assert(Count_ != 0); |
260 |
for (unsigned int i =0; i < 3; i++) { |
261 |
ret[i] = (Avg2_[i] - Avg_[i] * Avg_[i]); |
262 |
} |
263 |
return; |
264 |
} |
265 |
|
266 |
/** |
267 |
* compute error of average value |
268 |
*/ |
269 |
void getStdDev(ResultType &ret) { |
270 |
assert(Count_ != 0); |
271 |
ResultType var; |
272 |
this->getVariance(var); |
273 |
ret[0] = sqrt(var[0]); |
274 |
ret[1] = sqrt(var[1]); |
275 |
ret[2] = sqrt(var[2]); |
276 |
return; |
277 |
} |
278 |
|
279 |
/** |
280 |
* return the 95% confidence interval: |
281 |
* |
282 |
* That is returns c, such that we have 95% confidence that the |
283 |
* true mean is within 2c of the Average (x): |
284 |
* |
285 |
* x - c <= true mean <= x + c |
286 |
* |
287 |
*/ |
288 |
void get95percentConfidenceInterval(ResultType &ret) { |
289 |
assert(Count_ != 0); |
290 |
ResultType sd; |
291 |
this->getStdDev(sd); |
292 |
ret[0] = 1.960 * sd[0] / sqrt(Count_); |
293 |
ret[1] = 1.960 * sd[1] / sqrt(Count_); |
294 |
ret[2] = 1.960 * sd[2] / sqrt(Count_); |
295 |
return; |
296 |
} |
297 |
|
298 |
/** |
299 |
* return the largest length |
300 |
*/ |
301 |
void getMaxLength(RealType &ret) { |
302 |
assert(Count_ != 0); |
303 |
ret = Max_; |
304 |
return; |
305 |
} |
306 |
|
307 |
/** |
308 |
* return the smallest length |
309 |
*/ |
310 |
void getMinLength(RealType &ret) { |
311 |
assert(Count_ != 0); |
312 |
ret = Min_; |
313 |
return; |
314 |
} |
315 |
|
316 |
/** |
317 |
* return the largest length |
318 |
*/ |
319 |
void getAverageLength(RealType &ret) { |
320 |
assert(Count_ != 0); |
321 |
ret = AvgLen_; |
322 |
return; |
323 |
} |
324 |
|
325 |
/** |
326 |
* compute the Variance of the length |
327 |
*/ |
328 |
void getLengthVariance(RealType &ret) { |
329 |
assert(Count_ != 0); |
330 |
ret= (AvgLen2_ - AvgLen_ * AvgLen_); |
331 |
return; |
332 |
} |
333 |
|
334 |
/** |
335 |
* compute error of average value |
336 |
*/ |
337 |
void getLengthStdDev(RealType &ret) { |
338 |
assert(Count_ != 0); |
339 |
RealType var; |
340 |
this->getLengthVariance(var); |
341 |
ret = sqrt(var); |
342 |
return; |
343 |
} |
344 |
|
345 |
/** |
346 |
* return the 95% confidence interval: |
347 |
* |
348 |
* That is returns c, such that we have 95% confidence that the |
349 |
* true mean is within 2c of the Average (x): |
350 |
* |
351 |
* x - c <= true mean <= x + c |
352 |
* |
353 |
*/ |
354 |
void getLength95percentConfidenceInterval(ResultType &ret) { |
355 |
assert(Count_ != 0); |
356 |
RealType sd; |
357 |
this->getLengthStdDev(sd); |
358 |
ret = 1.960 * sd / sqrt(Count_); |
359 |
return; |
360 |
} |
361 |
|
362 |
|
363 |
private: |
364 |
ResultType Val_; |
365 |
ResultType Avg_; |
366 |
ResultType Avg2_; |
367 |
RealType AvgLen_; |
368 |
RealType AvgLen2_; |
369 |
RealType Min_; |
370 |
RealType Max_; |
371 |
|
372 |
}; |
373 |
|
374 |
class MatrixAccumulator : public BaseAccumulator { |
375 |
|
376 |
typedef Mat3x3d ElementType; |
377 |
typedef Mat3x3d ResultType; |
378 |
|
379 |
public: |
380 |
MatrixAccumulator() : BaseAccumulator() { |
381 |
this->clear(); |
382 |
} |
383 |
|
384 |
/** |
385 |
* Accumulate another value |
386 |
*/ |
387 |
void add(ElementType const& val) { |
388 |
Count_++; |
389 |
for (unsigned int i = 0; i < 3; i++) { |
390 |
for (unsigned int j = 0; j < 3; j++) { |
391 |
Avg_(i,j) += (val(i,j) - Avg_(i,j) ) / Count_; |
392 |
Avg2_(i,j) += (val(i,j) * val(i,j) - Avg2_(i,j)) / Count_; |
393 |
Val_(i,j) = val(i,j); |
394 |
} |
395 |
} |
396 |
} |
397 |
|
398 |
/** |
399 |
* reset the Accumulator to the empty state |
400 |
*/ |
401 |
void clear() { |
402 |
Count_ = 0; |
403 |
Avg_ *= 0.0; |
404 |
Avg2_ *= 0.0; |
405 |
Val_ *= 0.0; |
406 |
} |
407 |
|
408 |
/** |
409 |
* return the most recently added value |
410 |
*/ |
411 |
void getLastValue(ElementType &ret) { |
412 |
ret = Val_; |
413 |
return; |
414 |
} |
415 |
|
416 |
/** |
417 |
* compute the Mean |
418 |
*/ |
419 |
void getAverage(ResultType &ret) { |
420 |
assert(Count_ != 0); |
421 |
ret = Avg_; |
422 |
return; |
423 |
} |
424 |
|
425 |
/** |
426 |
* compute the Variance |
427 |
*/ |
428 |
void getVariance(ResultType &ret) { |
429 |
assert(Count_ != 0); |
430 |
for (unsigned int i = 0; i < 3; i++) { |
431 |
for (unsigned int j = 0; j < 3; j++) { |
432 |
ret(i,j) = (Avg2_(i,j) - Avg_(i,j) * Avg_(i,j)); |
433 |
} |
434 |
} |
435 |
return; |
436 |
} |
437 |
|
438 |
/** |
439 |
* compute error of average value |
440 |
*/ |
441 |
void getStdDev(ResultType &ret) { |
442 |
assert(Count_ != 0); |
443 |
Mat3x3d var; |
444 |
this->getVariance(var); |
445 |
for (unsigned int i = 0; i < 3; i++) { |
446 |
for (unsigned int j = 0; j < 3; j++) { |
447 |
ret(i,j) = sqrt(var(i,j)); |
448 |
} |
449 |
} |
450 |
return; |
451 |
} |
452 |
|
453 |
/** |
454 |
* return the 95% confidence interval: |
455 |
* |
456 |
* That is returns c, such that we have 95% confidence that the |
457 |
* true mean is within 2c of the Average (x): |
458 |
* |
459 |
* x - c <= true mean <= x + c |
460 |
* |
461 |
*/ |
462 |
void get95percentConfidenceInterval(ResultType &ret) { |
463 |
assert(Count_ != 0); |
464 |
Mat3x3d sd; |
465 |
this->getStdDev(sd); |
466 |
for (unsigned int i = 0; i < 3; i++) { |
467 |
for (unsigned int j = 0; j < 3; j++) { |
468 |
ret(i,j) = 1.960 * sd(i,j) / sqrt(Count_); |
469 |
} |
470 |
} |
471 |
return; |
472 |
} |
473 |
|
474 |
private: |
475 |
ElementType Val_; |
476 |
ResultType Avg_; |
477 |
ResultType Avg2_; |
478 |
}; |
479 |
|
480 |
|
481 |
} |
482 |
|
483 |
#endif |