DGtal  1.5.beta
MultiStatistics.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file MultiStatistics.ih
19  * @author Bertrand Kerautret (\c kerautre@loria.fr )
20  * LORIA (CNRS, UMR 7503), University of Nancy, France
21  *
22  * @author Jacques-Olivier Lachaud
23  * @date 2015/11/09
24  *
25  * Implementation of inline methods defined in MultiStatistics.h
26  *
27  * This file is part of the DGtal library.
28  */
29 
30 ///////////////////////////////////////////////////////////////////////////////
31 // IMPLEMENTATION of inline methods.
32 ///////////////////////////////////////////////////////////////////////////////
33 
34 //////////////////////////////////////////////////////////////////////////////
35 #include <cstdlib>
36 #include <iostream>
37 //////////////////////////////////////////////////////////////////////////////
38 
39 
40 
41 
42 ///////////////////////////////////////////////////////////////////////////////
43 // Implementation of inline
44 
45 
46 inline
47 DGtal::MultiStatistics::MultiStatistics( const unsigned int size,
48  const bool storeSamples )
49 {
50  mySamples = 0;
51  myExp = 0;
52  myExp2 = 0;
53  myVar = 0;
54  myUnbiasedVar = 0;
55  myMax = 0;
56  myMin = 0;
57  myValues = 0;
58  myIndiceMax = 0;
59  myIndiceMin = 0;
60  myIsTerminate = false;
61  init( size, storeSamples );
62 }
63 
64 
65 
66 
67 
68 
69 
70 inline
71 DGtal::MultiStatistics::~MultiStatistics()
72 {
73  erase();
74 }
75 
76 
77 
78 
79 void
80 DGtal::MultiStatistics::read( std::istream & in, MultiStatistics & samples,
81  const std::vector<unsigned int> & indices )
82 {
83  std::string str;
84  getline( in, str );
85  while ( in.good() )
86  {
87  if ( ( str != "" )
88  && ( str[ 0 ] != '#' ) )
89  {
90  std::istringstream in_str( str );
91  unsigned int idx = 1;
92  double val;
93  while ( in_str.good() )
94  {
95  in_str >> val;
96  for ( unsigned int j = 0; j < indices.size(); ++j )
97  if ( indices[ j ] == idx )
98  {
99  samples.addValue( j, val );
100  // cout << "Adding " << val << " to " << j << endl;
101  }
102  ++idx;
103  }
104  }
105  getline( in, str );
106  }
107 }
108 
109 
110 
111 inline
112 unsigned int
113 DGtal::MultiStatistics::nb() const
114 {
115  return myNb;
116 }
117 
118 
119 inline
120 unsigned int
121 DGtal::MultiStatistics::samples( const unsigned int k ) const
122 {
123  return mySamples[ k ];
124 }
125 
126 
127 inline
128 double
129 DGtal::MultiStatistics::mean( const unsigned int k ) const
130 {
131  ASSERT(myIsTerminate);
132  return myExp[ k ];
133 }
134 
135 
136 
137 inline
138 double
139 DGtal::MultiStatistics::variance( const unsigned int k ) const
140 {
141  ASSERT(myIsTerminate);
142  return myVar[ k ];
143 }
144 
145 
146 
147 inline
148 double
149 DGtal::MultiStatistics::unbiasedVariance( const unsigned int k ) const
150 {
151  ASSERT(myIsTerminate);
152  return myUnbiasedVar[ k ];
153 }
154 
155 
156 inline
157 double
158 DGtal::MultiStatistics::max( const unsigned int k ) const
159 {
160  return myMax[ k ];
161 }
162 
163 
164 inline
165 unsigned int
166 DGtal::MultiStatistics::maxIndice( const unsigned int k ) const
167 {
168  return myIndiceMax[ k ];
169 }
170 
171 
172 
173 inline
174 double
175 DGtal::MultiStatistics::min( const unsigned int k ) const
176 {
177  return myMin[ k ];
178 }
179 
180 inline
181 unsigned int
182 DGtal::MultiStatistics::minIndice( const unsigned int k ) const
183 {
184  return myIndiceMin[ k ];
185 }
186 
187 
188 
189 inline
190 double
191 DGtal::MultiStatistics::value( const unsigned int k, const unsigned int i ) const
192 {
193  if ( myStoreSamples ) {
194  ASSERT( ( k < myNb )
195  && ( myValues != 0 )
196  && ( i < mySamples[ k ] ) );
197  return myValues[ k ][ i ];
198  }
199  return 0.0;
200 }
201 
202 
203 inline
204 void
205 DGtal::MultiStatistics::addValue( unsigned int k, double v )
206 {
207  ASSERT( k < myNb );
208  mySamples[ k ] += 1;
209  myExp[ k ] += v;
210  myExp2[ k ] += v*v;
211  if ( mySamples[ k ] == 1 )
212  {
213  myMax[ k ] = v;
214  myMin[ k ] = v;
215  myIndiceMax[k]=0;
216  myIndiceMin[k]=0;
217 
218  }
219  else if ( v > myMax[ k ] ){
220  myMax[ k ] = v;
221  myIndiceMax[k]=mySamples[k]-1;
222  }
223  else if ( v < myMin[ k ] ){
224  myMin[ k ] = v;
225  myIndiceMin[k]=mySamples[k]-1;
226  }
227  if ( myStoreSamples )
228  myValues[ k ].push_back( v );
229  myIsTerminate = false;
230 }
231 
232 
233 template <class Iter>
234 inline
235 void
236 DGtal::MultiStatistics::addValues( const unsigned int k, Iter b, Iter e )
237 {
238  for ( ; b != e; ++b )
239  addValue( k, *b );
240  myIsTerminate = false;
241 }
242 
243 
244 
245 void
246 DGtal::MultiStatistics::terminate()
247 {
248  for ( unsigned int k = 0; k < myNb; ++k )
249  {
250  myExp[ k ] /= mySamples[ k ];
251  myExp2[ k ] /= mySamples[ k ];
252  myVar[ k ] = myExp2[ k ] - myExp[ k ] * myExp[ k ];
253  if ( mySamples[ k ] > 1 )
254  myUnbiasedVar[ k ] = mySamples[ k ] * myVar[ k ] / ( mySamples[ k ] - 1 );
255  else
256  myUnbiasedVar[ k ] = 0;
257 
258  }
259  myIsTerminate = true;
260 }
261 
262 
263 void
264 DGtal::MultiStatistics::init( unsigned int size, bool storeSamples )
265 {
266  erase();
267  if ( size == 0 ) return;
268  myNb = size;
269  mySamples = new unsigned int[ size ];
270  myExp = new double[ size ];
271  myExp2 = new double[ size ];
272  myVar = new double[ size ];
273  myUnbiasedVar = new double[ size ];
274  myMax = new double[ size ];
275  myMin = new double[ size ];
276  myIndiceMax = new unsigned int[ size ];
277  myIndiceMin = new unsigned int[ size ];
278  myStoreSamples = storeSamples;
279  if ( myStoreSamples )
280  myValues = new std::vector<double>[ size ];
281  clear();
282  myIsTerminate = false;
283 }
284 
285 
286 
287 
288 void
289 DGtal::MultiStatistics::clear()
290 {
291  if ( myNb == 0 ) return;
292  for ( unsigned int i = 0; i < myNb; ++ i )
293  {
294  mySamples[ i ] = 0;
295  myExp[ i ] = 0.0;
296  myExp2[ i ] = 0.0;
297  myVar[ i ] = 0.0;
298  myUnbiasedVar[ i ] = 0.0;
299  myMax[ i ] = 0.0;
300  myMin[ i ] = 0.0;
301  myIndiceMin[ i ] = 0;
302  myIndiceMax[ i ] = 0;
303 
304  if ( myStoreSamples ) {
305  myValues[ i ].clear();
306  myValues[ i ].reserve( 128 );
307  }
308  }
309  myIsTerminate = false;
310 }
311 
312 
313 
314 void
315 DGtal::MultiStatistics::erase()
316 {
317  if ( mySamples != 0 ) delete[] mySamples;
318  if ( myExp != 0 ) delete[] myExp;
319  if ( myExp2 != 0 ) delete[] myExp2;
320  if ( myVar != 0 ) delete[] myVar;
321  if ( myUnbiasedVar != 0 ) delete[] myUnbiasedVar;
322  if ( myMax != 0 ) delete[] myMax;
323  if ( myMin != 0 ) delete[] myMin;
324  if ( myIndiceMax != 0 ) delete[] myIndiceMax;
325  if ( myIndiceMin != 0 ) delete[] myIndiceMin;
326  if ( myValues != 0 ) delete[] myValues;
327 
328  mySamples = 0;
329  myExp = 0;
330  myExp2 = 0;
331  myVar = 0;
332  myUnbiasedVar = 0;
333  myNb = 0;
334  myMin = 0;
335  myMax = 0;
336  myIndiceMin = 0;
337  myIndiceMax = 0;
338  myValues = 0;
339 
340  myStoreSamples = false;
341 }
342 
343 
344 
345 double
346 DGtal::MultiStatistics::covariance( const unsigned int x, const unsigned int y,
347  const unsigned int s, unsigned int e ) const
348 {
349  ASSERT( ( x < myNb ) && ( y < myNb ) && ( x != y )
350  && myStoreSamples
351  && ( samples( x ) == samples( y ) ) );
352  unsigned int size = samples( x );
353  double coVariance = 0.0;
354  double mx = 0.0;
355  double my = 0.0;
356 
357  if ( e == 0 ) e = size;
358  ASSERT( e > s );
359  for( unsigned int k = s; k != e; ++k )
360  {
361  coVariance += value( x, k ) * value( y, k );
362  mx += value( x, k );
363  my += value( y, k );
364  }
365  mx /= ( e - s );
366  my /= ( e - s );
367  coVariance = coVariance / ( e - s );
368  coVariance = coVariance - mx * my;
369  return coVariance;
370 }
371 
372 
373 
374 std::pair<double,double>
375 DGtal::MultiStatistics::linearRegression( const unsigned int x, const unsigned int y ) const
376 {
377  ASSERT( ( x < myNb ) && ( y < myNb ) && ( x != y )
378  && myStoreSamples
379  && ( samples( x ) == samples( y ) ) );
380 
381  double cov = covariance( x, y );
382  double slope = cov / unbiasedVariance( x );
383  double b = mean( y ) - slope * mean( x );
384  return std::make_pair( slope, b );
385 }
386 
387 
388 
389 double
390 DGtal::MultiStatistics::median( unsigned int k )
391 {
392  ASSERT( myStoreSamples );
393 
394  nth_element( myValues[ k ].begin(), myValues[ k ].begin()+(myValues[k].size()/2),
395  myValues[ k ].end());
396  return *(myValues[ k ].begin()+(myValues[k].size()/2));
397 }
398 
399 
400 
401 
402 ///////////////////////////////////////////////////////////////////////////////
403 // Interface - public :
404 
405 
406 
407 void
408 DGtal::MultiStatistics::selfDisplay( std::ostream& out ) const
409 {
410  out << "[Statistics] nb=" << nb() << std::endl;
411 }
412 
413 
414 bool
415 DGtal::MultiStatistics::isValid() const
416 {
417  return true;
418 }
419 
420 
421 
422 ///////////////////////////////////////////////////////////////////////////////
423 // Implementation of inline functions and external operators //
424 
425 
426 
427 
428 
429 
430 
431 /**
432  * Overloads 'operator<<' for displaying objects of class 'MultiStatistics'.
433  * @param out the output stream where the object is written.
434  * @param object the object of class 'MultiStatistics' to write.
435  * @return the output stream after the writing.
436  */
437 inline
438 std::ostream&
439 DGtal::operator<< ( std::ostream & out,
440  const MultiStatistics & object )
441 {
442  object.selfDisplay ( out );
443  return out;
444 }
445 
446 // //
447 ///////////////////////////////////////////////////////////////////////////////
448 
449