statistics.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
7 // version 2.1 of the License, or (at your option) any later version.
8
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17
18
19 // C++ includes
20 #include <algorithm> // for std::min_element, std::max_element
21 #include <fstream> // std::ofstream
22 #include <numeric> // std::accumulate
23
24 // Local includes
25 #include "libmesh/statistics.h"
27
28 namespace libMesh
29 {
30
31
32
33 // ------------------------------------------------------------
34 // StatisticsVector class member functions
35 template <typename T>
37 {
38  Real normsq = 0.;
39  const dof_id_type n = cast_int<dof_id_type>(this->size());
40  for (dof_id_type i = 0; i != n; ++i)
41  normsq += ((*this)[i] * (*this)[i]);
42
43  return std::sqrt(normsq);
44 }
45
46
47 template <typename T>
49 {
50  LOG_SCOPE ("minimum()", "StatisticsVector");
51
52  const T min = *(std::min_element(this->begin(), this->end()));
53
54  return min;
55 }
56
57
58
59
60 template <typename T>
62 {
63  LOG_SCOPE ("maximum()", "StatisticsVector");
64
65  const T max = *(std::max_element(this->begin(), this->end()));
66
67  return max;
68 }
69
70
71
72
73 template <typename T>
75 {
76  LOG_SCOPE ("mean()", "StatisticsVector");
77
78  const dof_id_type n = cast_int<dof_id_type>(this->size());
79
80  Real the_mean = 0;
81
82  for (dof_id_type i=0; i<n; i++)
83  {
84  the_mean += ( static_cast<Real>((*this)[i]) - the_mean ) /
85  static_cast<Real>(i + 1);
86  }
87
88  return the_mean;
89 }
90
91
92
93
94 template <typename T>
96 {
97  const dof_id_type n = cast_int<dof_id_type>(this->size());
98
99  if (n == 0)
100  return 0.;
101
102  LOG_SCOPE ("median()", "StatisticsVector");
103
104  std::sort(this->begin(), this->end());
105
106  const dof_id_type lhs = (n-1) / 2;
107  const dof_id_type rhs = n / 2;
108
109  Real the_median = 0;
110
111
112  if (lhs == rhs)
113  {
114  the_median = static_cast<Real>((*this)[lhs]);
115  }
116
117  else
118  {
119  the_median = ( static_cast<Real>((*this)[lhs]) +
120  static_cast<Real>((*this)[rhs]) ) / 2.0;
121  }
122
123  return the_median;
124 }
125
126
127
128
129 template <typename T>
131 {
132  StatisticsVector<T> sv = (*this);
133
134  return sv.median();
135 }
136
137
138
139
140 template <typename T>
142 {
143  const dof_id_type n = cast_int<dof_id_type>(this->size());
144
145  LOG_SCOPE ("variance()", "StatisticsVector");
146
147  Real the_variance = 0;
148
149  for (dof_id_type i=0; i<n; i++)
150  {
151  const Real delta = ( static_cast<Real>((*this)[i]) - mean_in );
152  the_variance += (delta * delta - the_variance) /
153  static_cast<Real>(i + 1);
154  }
155
156  if (n > 1)
157  the_variance *= static_cast<Real>(n) / static_cast<Real>(n - 1);
158
159  return the_variance;
160 }
161
162
163 template <typename T>
165 {
166  const dof_id_type n = cast_int<dof_id_type>(this->size());
167  const Real max = this->maximum();
168
169  for (dof_id_type i=0; i<n; i++)
170  (*this)[i] = static_cast<T>((*this)[i] / max);
171 }
172
173
174
175
176
177 template <typename T>
178 void StatisticsVector<T>::histogram(std::vector<dof_id_type> & bin_members,
179  unsigned int n_bins)
180 {
181  // Must have at least 1 bin
182  libmesh_assert (n_bins>0);
183
184  const dof_id_type n = cast_int<dof_id_type>(this->size());
185
186  std::sort(this->begin(), this->end());
187
188  // The StatisticsVector can hold both integer and float types.
189  // We will define all the bins, etc. using Reals.
190  Real min = static_cast<Real>(this->minimum());
191  Real max = static_cast<Real>(this->maximum());
192  Real bin_size = (max - min) / static_cast<Real>(n_bins);
193
194  LOG_SCOPE ("histogram()", "StatisticsVector");
195
196  std::vector<Real> bin_bounds(n_bins+1);
197  for (std::size_t i=0; i<bin_bounds.size(); i++)
198  bin_bounds[i] = min + Real(i) * bin_size;
199
200  // Give the last bin boundary a little wiggle room: we don't want
201  // it to be just barely less than the max, otherwise our bin test below
202  // may fail.
203  bin_bounds.back() += 1.e-6 * bin_size;
204
205  // This vector will store the number of members each bin has.
206  bin_members.resize(n_bins);
207
208  dof_id_type data_index = 0;
209  for (std::size_t j=0; j<bin_members.size(); j++) // bin vector indexing
210  {
211  // libMesh::out << "(debug) Filling bin " << j << std::endl;
212
213  for (dof_id_type i=data_index; i<n; i++) // data vector indexing
214  {
215  //libMesh::out << "(debug) Processing index=" << i << std::endl;
216  Real current_val = static_cast<Real>( (*this)[i] );
217
218  // There may be entries in the vector smaller than the value
219  // reported by this->minimum(). (e.g. inactive elements in an
220  // ErrorVector.) We just skip entries like that.
221  if (current_val < min)
222  {
223  // libMesh::out << "(debug) Skipping entry v[" << i << "]="
224  // << (*this)[i]
225  // << " which is less than the min value: min="
226  // << min << std::endl;
227  continue;
228  }
229
230  if (current_val > bin_bounds[j+1]) // if outside the current bin (bin[j] is bounded
231  // by bin_bounds[j] and bin_bounds[j+1])
232  {
233  // libMesh::out.precision(16);
234  // libMesh::out.setf(std::ios_base::fixed);
235  // libMesh::out << "(debug) (*this)[i]= " << (*this)[i]
236  // << " is greater than bin_bounds[j+1]="
237  // << bin_bounds[j+1] << std::endl;
238  data_index = i; // start searching here for next bin
239  break; // go to next bin
240  }
241
242  // Otherwise, increment current bin's count
243  bin_members[j]++;
244  // libMesh::out << "(debug) Binned index=" << i << std::endl;
245  }
246  }
247
248 #ifdef DEBUG
249  // Check the number of binned entries
250  const dof_id_type n_binned = std::accumulate(bin_members.begin(),
251  bin_members.end(),
252  static_cast<dof_id_type>(0),
253  std::plus<dof_id_type>());
254
255  if (n != n_binned)
256  {
257  libMesh::out << "Warning: The number of binned entries, n_binned="
258  << n_binned
259  << ", did not match the total number of entries, n="
260  << n << "." << std::endl;
261  }
262 #endif
263 }
264
265
266
267
268
269 template <typename T>
271  const std::string & filename,
272  unsigned int n_bins)
273 {
274  // First generate the histogram with the desired number of bins
275  std::vector<dof_id_type> bin_members;
276  this->histogram(bin_members, n_bins);
277
278  // The max, min and bin size are used to generate x-axis values.
279  T min = this->minimum();
280  T max = this->maximum();
281  T bin_size = (max - min) / static_cast<T>(n_bins);
282
283  // On processor 0: Write histogram to file
284  if (my_procid==0)
285  {
286  std::ofstream out_stream (filename.c_str());
287
288  out_stream << "clear all\n";
289  out_stream << "clf\n";
290  //out_stream << "x=linspace(" << min << "," << max << "," << n_bins+1 << ");\n";
291
292  // abscissa values are located at the center of each bin.
293  out_stream << "x=[";
294  for (std::size_t i=0; i<bin_members.size(); ++i)
295  {
296  out_stream << min + (Real(i)+0.5)*bin_size << " ";
297  }
298  out_stream << "];\n";
299
300  out_stream << "y=[";
301  for (std::size_t i=0; i<bin_members.size(); ++i)
302  {
303  out_stream << bin_members[i] << " ";
304  }
305  out_stream << "];\n";
306  out_stream << "bar(x,y);\n";
307  }
308 }
309
310
311
312 template <typename T>
313 void StatisticsVector<T>::histogram(std::vector<dof_id_type> & bin_members,
314  unsigned int n_bins) const
315 {
316  StatisticsVector<T> sv = (*this);
317
318  return sv.histogram(bin_members, n_bins);
319 }
320
321
322
323
324 template <typename T>
325 std::vector<dof_id_type> StatisticsVector<T>::cut_below(Real cut) const
326 {
327  LOG_SCOPE ("cut_below()", "StatisticsVector");
328
329  const dof_id_type n = cast_int<dof_id_type>(this->size());
330
331  std::vector<dof_id_type> cut_indices;
332  cut_indices.reserve(n/2); // Arbitrary
333
334  for (dof_id_type i=0; i<n; i++)
335  {
336  if ((*this)[i] < cut)
337  {
338  cut_indices.push_back(i);
339  }
340  }
341
342  return cut_indices;
343 }
344
345
346
347
348 template <typename T>
349 std::vector<dof_id_type> StatisticsVector<T>::cut_above(Real cut) const
350 {
351  LOG_SCOPE ("cut_above()", "StatisticsVector");
352
353  const dof_id_type n = cast_int<dof_id_type>(this->size());
354
355  std::vector<dof_id_type> cut_indices;
356  cut_indices.reserve(n/2); // Arbitrary
357
358  for (dof_id_type i=0; i<n; i++)
359  if ((*this)[i] > cut)
360  cut_indices.push_back(i);
361
362  return cut_indices;
363 }
364
365
366
367
368 //------------------------------------------------------------
369 // Explicit Instantiations
370 template class StatisticsVector<float>;
371 template class StatisticsVector<double>;
372 #ifdef LIBMESH_DEFAULT_TRIPLE_PRECISION
373 template class StatisticsVector<long double>;
374 #endif
375 template class StatisticsVector<int>;
376 template class StatisticsVector<unsigned int>;
377
378 } // namespace libMesh
virtual T maximum() const
Definition: statistics.C:61
virtual Real mean() const
Definition: statistics.C:74
virtual Real median()
Definition: statistics.C:95
virtual Real l2_norm() const
Definition: statistics.C:36
uint8_t processor_id_type
Definition: id_types.h:99
A std::vector derived class for implementing simple statistical algorithms.
Definition: statistics.h:67
virtual std::vector< dof_id_type > cut_above(Real cut) const
Definition: statistics.C:349
IterBase * end
long double max(long double a, double b)
void plot_histogram(const processor_id_type my_procid, const std::string &filename, unsigned int n_bins)
Definition: statistics.C:270
virtual void histogram(std::vector< dof_id_type > &bin_members, unsigned int n_bins=10)
Definition: statistics.C:178
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual std::vector< dof_id_type > cut_below(Real cut) const
Definition: statistics.C:325
virtual T minimum() const
Definition: statistics.C:48
OStreamProxy out(std::cout)
long double min(long double a, double b)
virtual Real variance() const
Definition: statistics.h:134
uint8_t dof_id_type
Definition: id_types.h:64