libMesh
statistics.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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"
26 #include "libmesh/libmesh_logging.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 + 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 + (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 Real median()
Definition: statistics.C:95
virtual Real l2_norm() const
Definition: statistics.C:36
virtual T maximum() const
Definition: statistics.C:61
void normalize()
Divides all entries by the largest entry and stores the result.
Definition: statistics.C:164
uint8_t processor_id_type
Definition: id_types.h:99
virtual Real mean() const
Definition: statistics.C:74
virtual T minimum() const
Definition: statistics.C:48
The StatisticsVector class is derived from the std::vector<> and therefore has all of its useful feat...
Definition: statistics.h:67
IterBase * end
Also have a polymorphic pointer to the end object, this prevents iterating past the end...
The libMesh namespace provides an interface to certain functionality in the library.
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)
Generates a Matlab/Octave style file which can be used to make a plot of the histogram having the des...
Definition: statistics.C:270
libmesh_assert(j)
virtual std::vector< dof_id_type > cut_below(Real cut) const
Definition: statistics.C:325
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
OStreamProxy out
virtual Real variance() const
Definition: statistics.h:134
long double min(long double a, double b)
uint8_t dof_id_type
Definition: id_types.h:64
virtual std::vector< dof_id_type > cut_above(Real cut) const
Definition: statistics.C:349