DGtal  0.9.2
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  myUnbiasedVar[ k ] = mySamples[ k ] * myVar[ k ]
254  / ( mySamples[ k ] - 1 );
255  }
256  myIsTerminate = true;
257 }
258 
259 
260 void
261 DGtal::MultiStatistics::init( unsigned int size, bool storeSamples )
262 {
263  erase();
264  if ( size == 0 ) return;
265  myNb = size;
266  mySamples = new unsigned int[ size ];
267  myExp = new double[ size ];
268  myExp2 = new double[ size ];
269  myVar = new double[ size ];
270  myUnbiasedVar = new double[ size ];
271  myMax = new double[ size ];
272  myMin = new double[ size ];
273  myIndiceMax = new unsigned int[ size ];
274  myIndiceMin = new unsigned int[ size ];
275  myStoreSamples = storeSamples;
276  if ( myStoreSamples )
277  myValues = new std::vector<double>[ size ];
278  clear();
279  myIsTerminate = false;
280 }
281 
282 
283 
284 
285 void
286 DGtal::MultiStatistics::clear()
287 {
288  if ( myNb == 0 ) return;
289  for ( unsigned int i = 0; i < myNb; ++ i )
290  {
291  mySamples[ i ] = 0;
292  myExp[ i ] = 0.0;
293  myExp2[ i ] = 0.0;
294  myVar[ i ] = 0.0;
295  myUnbiasedVar[ i ] = 0.0;
296  myMax[ i ] = 0.0;
297  myMin[ i ] = 0.0;
298  myIndiceMin[ i ] = 0;
299  myIndiceMax[ i ] = 0;
300 
301  if ( myStoreSamples ) {
302  myValues[ i ].clear();
303  myValues[ i ].reserve( 128 );
304  }
305  }
306  myIsTerminate = false;
307 }
308 
309 
310 
311 void
312 DGtal::MultiStatistics::erase()
313 {
314  if ( mySamples != 0 ) delete[] mySamples;
315  if ( myExp != 0 ) delete[] myExp;
316  if ( myExp2 != 0 ) delete[] myExp2;
317  if ( myVar != 0 ) delete[] myVar;
318  if ( myUnbiasedVar != 0 ) delete[] myUnbiasedVar;
319  if ( myMax != 0 ) delete[] myMax;
320  if ( myMin != 0 ) delete[] myMin;
321  if ( myIndiceMax != 0 ) delete[] myIndiceMax;
322  if ( myIndiceMin != 0 ) delete[] myIndiceMin;
323  if ( myValues != 0 ) delete[] myValues;
324 
325  mySamples = 0;
326  myExp = 0;
327  myExp2 = 0;
328  myVar = 0;
329  myUnbiasedVar = 0;
330  myNb = 0;
331  myMin = 0;
332  myMax = 0;
333  myIndiceMin = 0;
334  myIndiceMax = 0;
335  myValues = 0;
336 
337  myStoreSamples = false;
338 }
339 
340 
341 
342 double
343 DGtal::MultiStatistics::covariance( const unsigned int x, const unsigned int y,
344  const unsigned int s, unsigned int e ) const
345 {
346  ASSERT( ( x < myNb ) && ( y < myNb ) && ( x != y )
347  && myStoreSamples
348  && ( samples( x ) == samples( y ) ) );
349  unsigned int size = samples( x );
350  double coVariance = 0.0;
351  double mx = 0.0;
352  double my = 0.0;
353 
354  if ( e == 0 ) e = size;
355  ASSERT( e > s );
356  for( unsigned int k = s; k != e; ++k )
357  {
358  coVariance += value( x, k ) * value( y, k );
359  mx += value( x, k );
360  my += value( y, k );
361  }
362  mx /= ( e - s );
363  my /= ( e - s );
364  coVariance = coVariance / ( e - s );
365  coVariance = coVariance - mx * my;
366  return coVariance;
367 }
368 
369 
370 
371 std::pair<double,double>
372 DGtal::MultiStatistics::linearRegression( const unsigned int x, const unsigned int y ) const
373 {
374  ASSERT( ( x < myNb ) && ( y < myNb ) && ( x != y )
375  && myStoreSamples
376  && ( samples( x ) == samples( y ) ) );
377 
378  double cov = covariance( x, y );
379  double slope = cov / unbiasedVariance( x );
380  double b = mean( y ) - slope * mean( x );
381  return std::make_pair( slope, b );
382 }
383 
384 
385 
386 double
387 DGtal::MultiStatistics::median( unsigned int k )
388 {
389  ASSERT( myStoreSamples );
390 
391  nth_element( myValues[ k ].begin(), myValues[ k ].begin()+(myValues[k].size()/2),
392  myValues[ k ].end());
393  return *(myValues[ k ].begin()+(myValues[k].size()/2));
394 }
395 
396 
397 
398 
399 ///////////////////////////////////////////////////////////////////////////////
400 // Interface - public :
401 
402 
403 
404 void
405 DGtal::MultiStatistics::selfDisplay( std::ostream& out ) const
406 {
407  out << "[Statistics] nb=" << nb() << std::endl;
408 }
409 
410 
411 bool
412 DGtal::MultiStatistics::isValid() const
413 {
414  return true;
415 }
416 
417 
418 
419 ///////////////////////////////////////////////////////////////////////////////
420 // Implementation of inline functions and external operators //
421 
422 
423 
424 
425 
426 
427 
428 /**
429  * Overloads 'operator<<' for displaying objects of class 'MultiStatistics'.
430  * @param out the output stream where the object is written.
431  * @param object the object of class 'MultiStatistics' to write.
432  * @return the output stream after the writing.
433  */
434 inline
435 std::ostream&
436 DGtal::operator<< ( std::ostream & out,
437  const MultiStatistics & object )
438 {
439  object.selfDisplay ( out );
440  return out;
441 }
442 
443 // //
444 ///////////////////////////////////////////////////////////////////////////////
445 
446