src/utils/Accumulator.hpp

/*
 * Copyright (c) 2012 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 *
 * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
 * research, please cite the appropriate papers when you publish your
 * work.  Good starting points are:
 *                                                                      
 * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
 * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
 * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
 * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
 * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
 */

#ifndef UTILS_ACCUMULATOR_HPP
#define UTILS_ACCUMULATOR_HPP

#include <cmath>
#include <cassert>
#include "math/Vector3.hpp"

namespace OpenMD {


  class BaseAccumulator {
  public:
    virtual void clear() = 0;
    /**
     * get the number of accumulated values
     */
    virtual size_t count()  {
      return Count_;
    }
    virtual ~BaseAccumulator() {};
  protected:
    size_t Count_;

  };   


  /** 
   * Basic Accumulator class for numbers. 
   */  
  class Accumulator : public BaseAccumulator {    

    typedef RealType ElementType;
    typedef RealType ResultType;

  public:
    
    Accumulator() : BaseAccumulator() {
      this->clear();
    }
    
    ~Accumulator() {};

    /**
     * Accumulate another value
     */
    virtual void add(ElementType const& val) {
      Count_++;
      Avg_  += (val       - Avg_ ) / Count_;
      Avg2_ += (val * val - Avg2_) / Count_;
      Val_   = val;
      if (Count_ <= 1) {
        Max_ = val;
        Min_ = val;
      } else {
        Max_ = val > Max_ ? val : Max_;
        Min_ = val < Min_ ? val : Min_;
      }
    }
    
    /**
     * reset the Accumulator to the empty state
     */
    void clear() {
      Count_ = 0;
      Avg_   = 0;
      Avg2_  = 0;
      Val_   = 0;
    }
    

    /**
     * return the most recently added value
     */
    void getLastValue(ElementType &ret)  {
      ret = Val_;
      return;
    }    

    /**
     * compute the Mean
     */
    void getAverage(ResultType &ret)  {
      assert(Count_ != 0);
      ret = Avg_;
      return;
    }

    /**
     * compute the Variance
     */
    void getVariance(ResultType &ret)  {
      assert(Count_ != 0);
      ret = (Avg2_ - Avg_  * Avg_);
      return;
    }
    
    /**
     * compute error of average value
     */
    void getStdDev(ResultType &ret)  {
      assert(Count_ != 0);
      RealType var;
      this->getVariance(var);
      ret = sqrt(var);
      return;
    }

    /**
     * return the largest value
     */
    void getMax(ElementType &ret)  {
      assert(Count_ != 0);
      ret = Max_;
      return;
    }

    /**
     * return the smallest value
     */
    void getMin(ElementType &ret)  {
      assert(Count_ != 0);
      ret = Max_;
      return;
    }

    /**
     * return the 95% confidence interval:
     *
     * That is returns c, such that we have 95% confidence that the
     * true mean is within 2c of the Average (x):
     *
     *   x - c <= true mean <= x + c
     *
     */
    void get95percentConfidenceInterval(ResultType &ret) {
      assert(Count_ != 0);
      RealType sd;
      this->getStdDev(sd);
      ret = 1.960 * sd / sqrt(RealType(Count_));
      return;
    }

  private:
    ElementType Val_;
    ResultType Avg_;
    ResultType Avg2_;
    ElementType Min_;
    ElementType Max_;
  };

  class VectorAccumulator : public BaseAccumulator {
    
    typedef Vector3d ElementType;
    typedef Vector3d ResultType;
    
  public:
    VectorAccumulator() : BaseAccumulator() {
      this->clear();
    }

    /**
     * Accumulate another value
     */
    void add(ElementType const& val) {
      Count_++;
      RealType len(0.0);
      for (unsigned int i =0; i < 3; i++) {
        Avg_[i]  += (val[i]       - Avg_[i] ) / Count_;
        Avg2_[i] += (val[i] * val[i] - Avg2_[i]) / Count_;
        Val_[i]   = val[i];
        len += val[i]*val[i];
      }
      len = sqrt(len);
      AvgLen_  += (len       - AvgLen_ ) / Count_;
      AvgLen2_ += (len * len - AvgLen2_) / Count_;

      if (Count_ <= 1) {
        Max_ = len;
        Min_ = len;
      } else {
        Max_ = len > Max_ ? len : Max_;
        Min_ = len < Min_ ? len : Min_;
      }
    }

    /**
     * reset the Accumulator to the empty state
     */
    void clear() {
      Count_ = 0;
      Avg_ = V3Zero;
      Avg2_ = V3Zero;
      Val_ = V3Zero;
      AvgLen_   = 0;
      AvgLen2_  = 0;
    }
    
    /**
     * return the most recently added value
     */
    void getLastValue(ElementType &ret) {
      ret = Val_;
      return;
    }
    
    /**
     * compute the Mean
     */
    void getAverage(ResultType &ret) {
      assert(Count_ != 0);
      ret = Avg_;
      return;
    }
    
    /**
     * compute the Variance
     */
    void getVariance(ResultType &ret) {
      assert(Count_ != 0);
      for (unsigned int i =0; i < 3; i++) {
        ret[i] = (Avg2_[i] - Avg_[i]  * Avg_[i]);
      }
      return;
    }
    
    /**
     * compute error of average value
     */
    void getStdDev(ResultType &ret) {
      assert(Count_ != 0);
      ResultType var;
      this->getVariance(var);
      ret[0] = sqrt(var[0]);
      ret[1] = sqrt(var[1]);
      ret[2] = sqrt(var[2]);
      return;
    }

    /**
     * return the 95% confidence interval:
     *
     * That is returns c, such that we have 95% confidence that the
     * true mean is within 2c of the Average (x):
     *
     *   x - c <= true mean <= x + c
     *
     */
    void get95percentConfidenceInterval(ResultType &ret) {
      assert(Count_ != 0);
      ResultType sd;
      this->getStdDev(sd);
      ret[0] = 1.960 * sd[0] / sqrt(RealType(Count_));
      ret[1] = 1.960 * sd[1] / sqrt(RealType(Count_));
      ret[2] = 1.960 * sd[2] / sqrt(RealType(Count_));
      return;
    }

    /**
     * return the largest length
     */
    void getMaxLength(RealType &ret) {
      assert(Count_ != 0);
      ret = Max_;
      return;
    }

    /**
     * return the smallest length
     */
    void getMinLength(RealType &ret) {
      assert(Count_ != 0);
      ret = Min_;
      return;
    }

    /**
     * return the largest length
     */
    void getAverageLength(RealType &ret) {
      assert(Count_ != 0);
      ret = AvgLen_;
      return;
    }

    /**
     * compute the Variance of the length
     */
    void getLengthVariance(RealType &ret) {
      assert(Count_ != 0);      
      ret= (AvgLen2_ - AvgLen_ * AvgLen_);
      return;
    }
    
    /**
     * compute error of average value
     */
    void getLengthStdDev(RealType &ret) {
      assert(Count_ != 0);
      RealType var;
      this->getLengthVariance(var);
      ret = sqrt(var);
      return;
    }

    /**
     * return the 95% confidence interval:
     *
     * That is returns c, such that we have 95% confidence that the
     * true mean is within 2c of the Average (x):
     *
     *   x - c <= true mean <= x + c
     *
     */
    void getLength95percentConfidenceInterval(ResultType &ret) {
      assert(Count_ != 0);
      RealType sd;
      this->getLengthStdDev(sd);
      ret = 1.960 * sd / sqrt(RealType(Count_));
      return;
    }


  private:
    ResultType Val_;
    ResultType Avg_;
    ResultType Avg2_;
    RealType AvgLen_;
    RealType AvgLen2_;
    RealType Min_;
    RealType Max_;

  };

  class MatrixAccumulator : public BaseAccumulator {
    
    typedef Mat3x3d ElementType;
    typedef Mat3x3d ResultType;
    
  public:
    MatrixAccumulator() : BaseAccumulator() {
      this->clear();
    }

    /**
     * Accumulate another value
     */
    void add(ElementType const& val) {
      Count_++;
      for (unsigned int i = 0; i < 3; i++) {
        for (unsigned int j = 0; j < 3; j++) {          
          Avg_(i,j)  += (val(i,j)       - Avg_(i,j) ) / Count_;
          Avg2_(i,j) += (val(i,j) * val(i,j) - Avg2_(i,j)) / Count_;
          Val_(i,j)   = val(i,j);
        }
      }
    }

    /**
     * reset the Accumulator to the empty state
     */
    void clear() {
      Count_ = 0;
      Avg_ *= 0.0;
      Avg2_ *= 0.0;
      Val_ *= 0.0;
    }
    
    /**
     * return the most recently added value
     */
    void getLastValue(ElementType &ret) {
      ret = Val_;
      return;
    }
    
    /**
     * compute the Mean
     */
    void getAverage(ResultType &ret) {
      assert(Count_ != 0);
      ret = Avg_;
      return;
    }

    /**
     * compute the Variance
     */
    void getVariance(ResultType &ret) {
      assert(Count_ != 0);
      for (unsigned int i = 0; i < 3; i++) {
        for (unsigned int j = 0; j < 3; j++) {          
          ret(i,j) = (Avg2_(i,j) - Avg_(i,j)  * Avg_(i,j));
        }
      }
      return;
    }
    
    /**
     * compute error of average value
     */
    void getStdDev(ResultType &ret) {
      assert(Count_ != 0);
      Mat3x3d var;
      this->getVariance(var);
      for (unsigned int i = 0; i < 3; i++) {
        for (unsigned int j = 0; j < 3; j++) {
          ret(i,j) = sqrt(var(i,j));  
        }
      }
      return;
    }
      
    /**
     * return the 95% confidence interval:
     *
     * That is returns c, such that we have 95% confidence that the
     * true mean is within 2c of the Average (x):
     *
     *   x - c <= true mean <= x + c
     *
     */
    void get95percentConfidenceInterval(ResultType &ret) {
      assert(Count_ != 0);
      Mat3x3d sd;
      this->getStdDev(sd);
      for (unsigned int i = 0; i < 3; i++) {
        for (unsigned int j = 0; j < 3; j++) {
          ret(i,j) = 1.960 * sd(i,j) / sqrt(RealType(Count_));
        }
      }
      return;
    }

  private:
    ElementType Val_;
    ResultType Avg_;
    ResultType Avg2_;
  };


} 

#endif
Revision:	2071
Committed:	Sat Mar 7 21:41:51 2015 UTC (10 years, 2 months ago) by gezelter
File size:	11339 byte(s)
Log Message:	Reducing the number of warnings when using g++ to compile.
#	User	Rev	Content
1	gezelter	1765	/*
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	gezelter	1879	* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	gezelter	1765	* [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	gezelter	1791
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	gezelter	2071	virtual ~BaseAccumulator() {};
63	gezelter	1791	protected:
64			size_t Count_;
65
66			};
67
68
69
70	gezelter	1765	/**
71			* Basic Accumulator class for numbers.
72	gezelter	1791	*/
73			class Accumulator : public BaseAccumulator {
74	gezelter	1765
75			typedef RealType ElementType;
76			typedef RealType ResultType;
77
78			public:
79
80	gezelter	1791	Accumulator() : BaseAccumulator() {
81	gezelter	1765	this->clear();
82			}
83	gezelter	2071
84			~Accumulator() {};
85	gezelter	1765
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	gezelter	1979	/**
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	gezelter	1980	ret = 1.960 * sd / sqrt(RealType(Count_));
183	gezelter	1979	return;
184			}
185
186	gezelter	1765	private:
187			ElementType Val_;
188			ResultType Avg_;
189			ResultType Avg2_;
190			ElementType Min_;
191			ElementType Max_;
192			};
193
194	gezelter	1791	class VectorAccumulator : public BaseAccumulator {
195	gezelter	1765
196			typedef Vector3d ElementType;
197			typedef Vector3d ResultType;
198
199			public:
200	gezelter	1791	VectorAccumulator() : BaseAccumulator() {
201	gezelter	1765	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	gezelter	1979	* 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	gezelter	1980	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	gezelter	1979	return;
299			}
300
301			/**
302	gezelter	1765	* 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	gezelter	1879
348	gezelter	1979	/**
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	gezelter	1980	ret = 1.960 * sd / sqrt(RealType(Count_));
362	gezelter	1979	return;
363			}
364
365
366	gezelter	1765	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	gezelter	1791	class MatrixAccumulator : public BaseAccumulator {
378	gezelter	1765
379			typedef Mat3x3d ElementType;
380			typedef Mat3x3d ResultType;
381
382			public:
383	gezelter	1791	MatrixAccumulator() : BaseAccumulator() {
384	gezelter	1765	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	gezelter	1979
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	gezelter	1980	ret(i,j) = 1.960 * sd(i,j) / sqrt(RealType(Count_));
472	gezelter	1979	}
473			}
474			return;
475			}
476
477	gezelter	1765	private:
478			ElementType Val_;
479			ResultType Avg_;
480			ResultType Avg2_;
481			};
482
483
484			}
485
486			#endif