ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/utils/Accumulator.hpp
Revision: 1979
Committed: Sat Apr 5 20:56:01 2014 UTC (11 years ago) by gezelter
File size: 11216 byte(s)
Log Message:
* Changed the stdDev printouts in RNEMD stats to the 95% confidence intervals
* Added the ability to specify non-bonded constraints in a molecule
* Added the ability to do selection offsets in the pAngle staticProps module

File Contents

# Content
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

Properties

Name Value
svn:eol-style native