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
6 // License as published by the Free Software Foundation; either
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