| 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 | virtual ~BaseAccumulator() {}; | 
| 63 | protected: | 
| 64 | size_t Count_; | 
| 65 |  | 
| 66 | }; | 
| 67 |  | 
| 68 |  | 
| 69 |  | 
| 70 | /** | 
| 71 | * Basic Accumulator class for numbers. | 
| 72 | */ | 
| 73 | class Accumulator : public BaseAccumulator { | 
| 74 |  | 
| 75 | typedef RealType ElementType; | 
| 76 | typedef RealType ResultType; | 
| 77 |  | 
| 78 | public: | 
| 79 |  | 
| 80 | Accumulator() : BaseAccumulator() { | 
| 81 | this->clear(); | 
| 82 | } | 
| 83 |  | 
| 84 | ~Accumulator() {}; | 
| 85 |  | 
| 86 | /** | 
| 87 | * Accumulate another value | 
| 88 | */ | 
| 89 | virtual void add(ElementType const& val) { | 
| 90 | Count_++; | 
| 91 | Avg_  += (val       - Avg_ ) / Count_; | 
| 92 | Avg2_ += (val * val - Avg2_) / Count_; | 
| 93 | Val_   = val; | 
| 94 | if (Count_ <= 1) { | 
| 95 | Max_ = val; | 
| 96 | Min_ = val; | 
| 97 | } else { | 
| 98 | Max_ = val > Max_ ? val : Max_; | 
| 99 | Min_ = val < Min_ ? val : Min_; | 
| 100 | } | 
| 101 | } | 
| 102 |  | 
| 103 | /** | 
| 104 | * reset the Accumulator to the empty state | 
| 105 | */ | 
| 106 | void clear() { | 
| 107 | Count_ = 0; | 
| 108 | Avg_   = 0; | 
| 109 | Avg2_  = 0; | 
| 110 | Val_   = 0; | 
| 111 | } | 
| 112 |  | 
| 113 |  | 
| 114 | /** | 
| 115 | * return the most recently added value | 
| 116 | */ | 
| 117 | void getLastValue(ElementType &ret)  { | 
| 118 | ret = Val_; | 
| 119 | return; | 
| 120 | } | 
| 121 |  | 
| 122 | /** | 
| 123 | * compute the Mean | 
| 124 | */ | 
| 125 | void getAverage(ResultType &ret)  { | 
| 126 | assert(Count_ != 0); | 
| 127 | ret = Avg_; | 
| 128 | return; | 
| 129 | } | 
| 130 |  | 
| 131 | /** | 
| 132 | * compute the Variance | 
| 133 | */ | 
| 134 | void getVariance(ResultType &ret)  { | 
| 135 | assert(Count_ != 0); | 
| 136 | ret = (Avg2_ - Avg_  * Avg_); | 
| 137 | return; | 
| 138 | } | 
| 139 |  | 
| 140 | /** | 
| 141 | * compute error of average value | 
| 142 | */ | 
| 143 | void getStdDev(ResultType &ret)  { | 
| 144 | assert(Count_ != 0); | 
| 145 | RealType var; | 
| 146 | this->getVariance(var); | 
| 147 | ret = sqrt(var); | 
| 148 | return; | 
| 149 | } | 
| 150 |  | 
| 151 | /** | 
| 152 | * return the largest value | 
| 153 | */ | 
| 154 | void getMax(ElementType &ret)  { | 
| 155 | assert(Count_ != 0); | 
| 156 | ret = Max_; | 
| 157 | return; | 
| 158 | } | 
| 159 |  | 
| 160 | /** | 
| 161 | * return the smallest value | 
| 162 | */ | 
| 163 | void getMin(ElementType &ret)  { | 
| 164 | assert(Count_ != 0); | 
| 165 | ret = Max_; | 
| 166 | return; | 
| 167 | } | 
| 168 |  | 
| 169 | /** | 
| 170 | * return the 95% confidence interval: | 
| 171 | * | 
| 172 | * That is returns c, such that we have 95% confidence that the | 
| 173 | * true mean is within 2c of the Average (x): | 
| 174 | * | 
| 175 | *   x - c <= true mean <= x + c | 
| 176 | * | 
| 177 | */ | 
| 178 | void get95percentConfidenceInterval(ResultType &ret) { | 
| 179 | assert(Count_ != 0); | 
| 180 | RealType sd; | 
| 181 | this->getStdDev(sd); | 
| 182 | ret = 1.960 * sd / sqrt(RealType(Count_)); | 
| 183 | return; | 
| 184 | } | 
| 185 |  | 
| 186 | private: | 
| 187 | ElementType Val_; | 
| 188 | ResultType Avg_; | 
| 189 | ResultType Avg2_; | 
| 190 | ElementType Min_; | 
| 191 | ElementType Max_; | 
| 192 | }; | 
| 193 |  | 
| 194 | class VectorAccumulator : public BaseAccumulator { | 
| 195 |  | 
| 196 | typedef Vector3d ElementType; | 
| 197 | typedef Vector3d ResultType; | 
| 198 |  | 
| 199 | public: | 
| 200 | VectorAccumulator() : BaseAccumulator() { | 
| 201 | this->clear(); | 
| 202 | } | 
| 203 |  | 
| 204 | /** | 
| 205 | * Accumulate another value | 
| 206 | */ | 
| 207 | void add(ElementType const& val) { | 
| 208 | Count_++; | 
| 209 | RealType len(0.0); | 
| 210 | for (unsigned int i =0; i < 3; i++) { | 
| 211 | Avg_[i]  += (val[i]       - Avg_[i] ) / Count_; | 
| 212 | Avg2_[i] += (val[i] * val[i] - Avg2_[i]) / Count_; | 
| 213 | Val_[i]   = val[i]; | 
| 214 | len += val[i]*val[i]; | 
| 215 | } | 
| 216 | len = sqrt(len); | 
| 217 | AvgLen_  += (len       - AvgLen_ ) / Count_; | 
| 218 | AvgLen2_ += (len * len - AvgLen2_) / Count_; | 
| 219 |  | 
| 220 | if (Count_ <= 1) { | 
| 221 | Max_ = len; | 
| 222 | Min_ = len; | 
| 223 | } else { | 
| 224 | Max_ = len > Max_ ? len : Max_; | 
| 225 | Min_ = len < Min_ ? len : Min_; | 
| 226 | } | 
| 227 | } | 
| 228 |  | 
| 229 | /** | 
| 230 | * reset the Accumulator to the empty state | 
| 231 | */ | 
| 232 | void clear() { | 
| 233 | Count_ = 0; | 
| 234 | Avg_ = V3Zero; | 
| 235 | Avg2_ = V3Zero; | 
| 236 | Val_ = V3Zero; | 
| 237 | AvgLen_   = 0; | 
| 238 | AvgLen2_  = 0; | 
| 239 | } | 
| 240 |  | 
| 241 | /** | 
| 242 | * return the most recently added value | 
| 243 | */ | 
| 244 | void getLastValue(ElementType &ret) { | 
| 245 | ret = Val_; | 
| 246 | return; | 
| 247 | } | 
| 248 |  | 
| 249 | /** | 
| 250 | * compute the Mean | 
| 251 | */ | 
| 252 | void getAverage(ResultType &ret) { | 
| 253 | assert(Count_ != 0); | 
| 254 | ret = Avg_; | 
| 255 | return; | 
| 256 | } | 
| 257 |  | 
| 258 | /** | 
| 259 | * compute the Variance | 
| 260 | */ | 
| 261 | void getVariance(ResultType &ret) { | 
| 262 | assert(Count_ != 0); | 
| 263 | for (unsigned int i =0; i < 3; i++) { | 
| 264 | ret[i] = (Avg2_[i] - Avg_[i]  * Avg_[i]); | 
| 265 | } | 
| 266 | return; | 
| 267 | } | 
| 268 |  | 
| 269 | /** | 
| 270 | * compute error of average value | 
| 271 | */ | 
| 272 | void getStdDev(ResultType &ret) { | 
| 273 | assert(Count_ != 0); | 
| 274 | ResultType var; | 
| 275 | this->getVariance(var); | 
| 276 | ret[0] = sqrt(var[0]); | 
| 277 | ret[1] = sqrt(var[1]); | 
| 278 | ret[2] = sqrt(var[2]); | 
| 279 | return; | 
| 280 | } | 
| 281 |  | 
| 282 | /** | 
| 283 | * return the 95% confidence interval: | 
| 284 | * | 
| 285 | * That is returns c, such that we have 95% confidence that the | 
| 286 | * true mean is within 2c of the Average (x): | 
| 287 | * | 
| 288 | *   x - c <= true mean <= x + c | 
| 289 | * | 
| 290 | */ | 
| 291 | void get95percentConfidenceInterval(ResultType &ret) { | 
| 292 | assert(Count_ != 0); | 
| 293 | ResultType sd; | 
| 294 | this->getStdDev(sd); | 
| 295 | ret[0] = 1.960 * sd[0] / sqrt(RealType(Count_)); | 
| 296 | ret[1] = 1.960 * sd[1] / sqrt(RealType(Count_)); | 
| 297 | ret[2] = 1.960 * sd[2] / sqrt(RealType(Count_)); | 
| 298 | return; | 
| 299 | } | 
| 300 |  | 
| 301 | /** | 
| 302 | * return the largest length | 
| 303 | */ | 
| 304 | void getMaxLength(RealType &ret) { | 
| 305 | assert(Count_ != 0); | 
| 306 | ret = Max_; | 
| 307 | return; | 
| 308 | } | 
| 309 |  | 
| 310 | /** | 
| 311 | * return the smallest length | 
| 312 | */ | 
| 313 | void getMinLength(RealType &ret) { | 
| 314 | assert(Count_ != 0); | 
| 315 | ret = Min_; | 
| 316 | return; | 
| 317 | } | 
| 318 |  | 
| 319 | /** | 
| 320 | * return the largest length | 
| 321 | */ | 
| 322 | void getAverageLength(RealType &ret) { | 
| 323 | assert(Count_ != 0); | 
| 324 | ret = AvgLen_; | 
| 325 | return; | 
| 326 | } | 
| 327 |  | 
| 328 | /** | 
| 329 | * compute the Variance of the length | 
| 330 | */ | 
| 331 | void getLengthVariance(RealType &ret) { | 
| 332 | assert(Count_ != 0); | 
| 333 | ret= (AvgLen2_ - AvgLen_ * AvgLen_); | 
| 334 | return; | 
| 335 | } | 
| 336 |  | 
| 337 | /** | 
| 338 | * compute error of average value | 
| 339 | */ | 
| 340 | void getLengthStdDev(RealType &ret) { | 
| 341 | assert(Count_ != 0); | 
| 342 | RealType var; | 
| 343 | this->getLengthVariance(var); | 
| 344 | ret = sqrt(var); | 
| 345 | return; | 
| 346 | } | 
| 347 |  | 
| 348 | /** | 
| 349 | * return the 95% confidence interval: | 
| 350 | * | 
| 351 | * That is returns c, such that we have 95% confidence that the | 
| 352 | * true mean is within 2c of the Average (x): | 
| 353 | * | 
| 354 | *   x - c <= true mean <= x + c | 
| 355 | * | 
| 356 | */ | 
| 357 | void getLength95percentConfidenceInterval(ResultType &ret) { | 
| 358 | assert(Count_ != 0); | 
| 359 | RealType sd; | 
| 360 | this->getLengthStdDev(sd); | 
| 361 | ret = 1.960 * sd / sqrt(RealType(Count_)); | 
| 362 | return; | 
| 363 | } | 
| 364 |  | 
| 365 |  | 
| 366 | private: | 
| 367 | ResultType Val_; | 
| 368 | ResultType Avg_; | 
| 369 | ResultType Avg2_; | 
| 370 | RealType AvgLen_; | 
| 371 | RealType AvgLen2_; | 
| 372 | RealType Min_; | 
| 373 | RealType Max_; | 
| 374 |  | 
| 375 | }; | 
| 376 |  | 
| 377 | class MatrixAccumulator : public BaseAccumulator { | 
| 378 |  | 
| 379 | typedef Mat3x3d ElementType; | 
| 380 | typedef Mat3x3d ResultType; | 
| 381 |  | 
| 382 | public: | 
| 383 | MatrixAccumulator() : BaseAccumulator() { | 
| 384 | this->clear(); | 
| 385 | } | 
| 386 |  | 
| 387 | /** | 
| 388 | * Accumulate another value | 
| 389 | */ | 
| 390 | void add(ElementType const& val) { | 
| 391 | Count_++; | 
| 392 | for (unsigned int i = 0; i < 3; i++) { | 
| 393 | for (unsigned int j = 0; j < 3; j++) { | 
| 394 | Avg_(i,j)  += (val(i,j)       - Avg_(i,j) ) / Count_; | 
| 395 | Avg2_(i,j) += (val(i,j) * val(i,j) - Avg2_(i,j)) / Count_; | 
| 396 | Val_(i,j)   = val(i,j); | 
| 397 | } | 
| 398 | } | 
| 399 | } | 
| 400 |  | 
| 401 | /** | 
| 402 | * reset the Accumulator to the empty state | 
| 403 | */ | 
| 404 | void clear() { | 
| 405 | Count_ = 0; | 
| 406 | Avg_ *= 0.0; | 
| 407 | Avg2_ *= 0.0; | 
| 408 | Val_ *= 0.0; | 
| 409 | } | 
| 410 |  | 
| 411 | /** | 
| 412 | * return the most recently added value | 
| 413 | */ | 
| 414 | void getLastValue(ElementType &ret) { | 
| 415 | ret = Val_; | 
| 416 | return; | 
| 417 | } | 
| 418 |  | 
| 419 | /** | 
| 420 | * compute the Mean | 
| 421 | */ | 
| 422 | void getAverage(ResultType &ret) { | 
| 423 | assert(Count_ != 0); | 
| 424 | ret = Avg_; | 
| 425 | return; | 
| 426 | } | 
| 427 |  | 
| 428 | /** | 
| 429 | * compute the Variance | 
| 430 | */ | 
| 431 | void getVariance(ResultType &ret) { | 
| 432 | assert(Count_ != 0); | 
| 433 | for (unsigned int i = 0; i < 3; i++) { | 
| 434 | for (unsigned int j = 0; j < 3; j++) { | 
| 435 | ret(i,j) = (Avg2_(i,j) - Avg_(i,j)  * Avg_(i,j)); | 
| 436 | } | 
| 437 | } | 
| 438 | return; | 
| 439 | } | 
| 440 |  | 
| 441 | /** | 
| 442 | * compute error of average value | 
| 443 | */ | 
| 444 | void getStdDev(ResultType &ret) { | 
| 445 | assert(Count_ != 0); | 
| 446 | Mat3x3d var; | 
| 447 | this->getVariance(var); | 
| 448 | for (unsigned int i = 0; i < 3; i++) { | 
| 449 | for (unsigned int j = 0; j < 3; j++) { | 
| 450 | ret(i,j) = sqrt(var(i,j)); | 
| 451 | } | 
| 452 | } | 
| 453 | return; | 
| 454 | } | 
| 455 |  | 
| 456 | /** | 
| 457 | * return the 95% confidence interval: | 
| 458 | * | 
| 459 | * That is returns c, such that we have 95% confidence that the | 
| 460 | * true mean is within 2c of the Average (x): | 
| 461 | * | 
| 462 | *   x - c <= true mean <= x + c | 
| 463 | * | 
| 464 | */ | 
| 465 | void get95percentConfidenceInterval(ResultType &ret) { | 
| 466 | assert(Count_ != 0); | 
| 467 | Mat3x3d sd; | 
| 468 | this->getStdDev(sd); | 
| 469 | for (unsigned int i = 0; i < 3; i++) { | 
| 470 | for (unsigned int j = 0; j < 3; j++) { | 
| 471 | ret(i,j) = 1.960 * sd(i,j) / sqrt(RealType(Count_)); | 
| 472 | } | 
| 473 | } | 
| 474 | return; | 
| 475 | } | 
| 476 |  | 
| 477 | private: | 
| 478 | ElementType Val_; | 
| 479 | ResultType Avg_; | 
| 480 | ResultType Avg2_; | 
| 481 | }; | 
| 482 |  | 
| 483 |  | 
| 484 | } | 
| 485 |  | 
| 486 | #endif |