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