OpenStructure
stat_accumulator.hh
Go to the documentation of this file.
1 //------------------------------------------------------------------------------
2 // This file is part of the OpenStructure project <www.openstructure.org>
3 //
4 // Copyright (C) 2008-2020 by the OpenStructure authors
5 //
6 // This library is free software; you can redistribute it and/or modify it under
7 // the terms of the GNU Lesser General Public License as published by the Free
8 // Software Foundation; either version 3.0 of the License, or (at your option)
9 // any later version.
10 // This library is distributed in the hope that it will be useful, but WITHOUT
11 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU Lesser General Public License
16 // along with this library; if not, write to the Free Software Foundation, Inc.,
17 // 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
18 //------------------------------------------------------------------------------
19 
20 #ifndef OST_STAT_ACCUMULATOR_HH
21 #define OST_STAT_ACCUMULATOR_HH
22 
23 #include <limits>
24 #include <boost/math/special_functions/binomial.hpp>
25 #include <ost/base.hh>
26 #include <ost/message.hh>
28 
29 namespace ost { namespace img { namespace alg {
30 
31 namespace {
32 
33 template<class D>
34 struct MinMax{
35  static bool less_cmp_(const D& v1, const D& v2){return v1<v2;}
36 };
37 template<>
38 struct MinMax<Complex>{
39  static bool less_cmp_(const Complex& v1, const Complex& v2){return abs(v1)<abs(v2);}
40 };
41 
42 } //anon ns
43 
44 template<unsigned int MAX_MOMENT=4,typename DATATYPE=Real>
46 {
47 public:
49  sum_(0.0),
50  sum2_(0.0),
51  max_(-std::numeric_limits<DATATYPE>::max()),
52  min_(std::numeric_limits<DATATYPE>::max()),
53  m_(),
54  w_(0.0),
55  n_(0)
56  {
57  for(unsigned int i=0;i<MAX_MOMENT;++i){
58  m_[i]=0.0;
59  }
60 
61  }
62 
63  StatAccumulator(DATATYPE val, Real w=1.0):
64  sum_(val),
65  sum2_(val*val),
66  max_(val),
67  min_(val),
68  m_(),
69  w_(w),
70  n_(1)
71  {
72  if(MAX_MOMENT>0){
73  m_[0]=val;
74  for(unsigned int i=1;i<MAX_MOMENT;++i){
75  m_[i]=0.0;
76  }
77  }
78  }
79 
80  void operator()(DATATYPE val, Real w=1.0)
81  {
82  *this+=StatAccumulator(val,w);
83  }
84 
86  {
88  acc+=*this;
89  return acc;
90  }
91 
93  {
94  if(0.0>=w_){
95  *this=acc;
96  return *this;
97  }
98  if(0.0>=acc.w_){
99  return *this;
100  }
101  n_+=acc.n_;
102  Real wn=w_+acc.w_;
103 
104  sum_+=acc.sum_;
105  sum2_+=acc.sum2_;
106  max_=std::max<DATATYPE>(max_,acc.max_,MinMax<DATATYPE>::less_cmp_);
107  min_=std::min<DATATYPE>(min_,acc.min_,MinMax<DATATYPE>::less_cmp_);
108  if(MAX_MOMENT>0){
109  DATATYPE delta=acc.m_[0]-m_[0];
110  DATATYPE delta_w=delta/wn;
111  if(MAX_MOMENT>2){
112  DATATYPE w1w2_delta_w=w_*acc.w_*delta_w;
113  DATATYPE w1w2_delta_wp=w1w2_delta_w*w1w2_delta_w;
114  Real iw2=1.0/acc.w_;
115  Real iw2pm1=iw2;
116  Real miw1=-1.0/w_;
117  Real miw1pm1=miw1;
118  DATATYPE mn[MAX_MOMENT]; // only MAX_MOMENT-2 values needed but overdimensioned to allow compilation for MAX_MOMENT==0 (even though it gets kicked out in the end by the dead code elimination)
119  for(unsigned int p=3;p<=MAX_MOMENT;++p){
120  w1w2_delta_wp*=w1w2_delta_w;
121  iw2pm1*=iw2;
122  miw1pm1*=miw1;
123  DATATYPE delta_wk=1.0;
124  DATATYPE s=0.0;
125  Real mw2k=1.0;
126  Real w1k=1.0;
127  for(unsigned int k=1;k<=p-2;++k){
128  w1k*=w_;
129  mw2k*=-acc.w_;
130  delta_wk*=delta_w;
131  s+=boost::math::binomial_coefficient<Real>(p, k)*(mw2k*m_[p-k-1]+w1k*acc.m_[p-k-1])*delta_wk;
132  }
133  mn[p-3]=acc.m_[p-1]+s+w1w2_delta_wp*(iw2pm1-miw1pm1);
134  }
135  for(unsigned int p=3;p<=MAX_MOMENT;++p){
136  m_[p-1]+=mn[p-3];
137  }
138  }
139  if(MAX_MOMENT>1){
140  m_[1]+=acc.m_[1]+delta_w*delta*acc.w_*w_;
141  }
142  m_[0]+=delta_w*acc.w_;
143  w_=wn;
144  }
145  return *this;
146  }
147 
149  {
150  return n_;
151  }
152 
153  DATATYPE GetMaximum() const
154  {
155  return max_;
156  }
157 
158  DATATYPE GetMinimum() const
159  {
160  return min_;
161  }
162 
163  Real GetWeight() const
164  {
165  return w_;
166  }
167 
168  DATATYPE GetMoment(unsigned int i) const
169  {
170  if(1>i){
171  throw Error("Invalid moment.");
172  }
173  if(MAX_MOMENT<i){
174  throw Error("Moment was not calculated.");
175  }
176  return m_[i-1];
177  }
178 
179  DATATYPE GetMean() const
180  {
181  if(MAX_MOMENT<1){
182  throw Error("Mean was not calculated.");
183  }
184  return m_[0];
185  }
186 
187  DATATYPE GetSum() const
188  {
189  return sum_;
190  }
191 
192  DATATYPE GetVariance() const
193  {
194  if(MAX_MOMENT<2){
195  throw Error("Variance was not calculated.");
196  }
197  if(n_==0.0){
198  return 0.0;
199  }
200  return m_[1]/(w_);
201  }
202 
203  DATATYPE GetStandardDeviation() const
204  {
205  return sqrt(GetVariance());
206  }
207 
208  DATATYPE GetRootMeanSquare() const
209  {
210  if(n_==0.0){
211  return 0.0;
212  }
213  return sqrt(sum2_/(w_));
214  }
215 
216  DATATYPE GetSkewness() const
217  {
218  if(MAX_MOMENT<3){
219  throw Error("Skewness was not calculated.");
220  }
221  if(m_[1]<1e-20){
222  return 0.0;
223  }else{
224  return m_[2]/sqrt(m_[1]*m_[1]*m_[1]);
225  }
226  }
227 
228  DATATYPE GetKurtosis() const
229  {
230  if(MAX_MOMENT<4){
231  throw Error("Kurtosis was not calculated.");
232  }
233  if(m_[1]<1e-20){
234  return 0.0;
235  }else{
236  return w_*m_[3] / (m_[1]*m_[1]);
237  }
238  }
239 
240 private:
241  DATATYPE sum_;
242  DATATYPE sum2_;
243  DATATYPE max_;
244  DATATYPE min_;
245  DATATYPE m_[MAX_MOMENT];
246  Real w_;
247  unsigned int n_;
248 };
249 
250 }}} //ns
251 #endif // OST_STAT_ACCUMULATOR_HH
DATATYPE GetMoment(unsigned int i) const
StatAccumulator(DATATYPE val, Real w=1.0)
StatAccumulator< MAX_MOMENT, DATATYPE > & operator+=(const StatAccumulator< MAX_MOMENT, DATATYPE > &acc)
StatAccumulator< MAX_MOMENT, DATATYPE > operator+(const StatAccumulator< MAX_MOMENT, DATATYPE > &acc2) const
void operator()(DATATYPE val, Real w=1.0)
float Real
Definition: base.hh:44
std::complex< Real > Complex
Definition: base.hh:51
Definition: base.dox:1