statistics.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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  for (dof_id_type i = 0; i != this->size(); ++i)
40  normsq += ((*this)[i] * (*this)[i]);
41 
42  return std::sqrt(normsq);
43 }
44 
45 
46 template <typename T>
48 {
49  START_LOG ("minimum()", "StatisticsVector");
50 
51  const T min = *(std::min_element(this->begin(), this->end()));
52 
53  STOP_LOG ("minimum()", "StatisticsVector");
54 
55  return min;
56 }
57 
58 
59 
60 
61 template <typename T>
63 {
64  START_LOG ("maximum()", "StatisticsVector");
65 
66  const T max = *(std::max_element(this->begin(), this->end()));
67 
68  STOP_LOG ("maximum()", "StatisticsVector");
69 
70  return max;
71 }
72 
73 
74 
75 
76 template <typename T>
78 {
79  START_LOG ("mean()", "StatisticsVector");
80 
81  const dof_id_type n = this->size();
82 
83  Real the_mean = 0;
84 
85  for (dof_id_type i=0; i<n; i++)
86  {
87  the_mean += ( static_cast<Real>((*this)[i]) - the_mean ) /
88  static_cast<Real>(i + 1);
89  }
90 
91  STOP_LOG ("mean()", "StatisticsVector");
92 
93  return the_mean;
94 }
95 
96 
97 
98 
99 template <typename T>
101 {
102  const dof_id_type n = this->size();
103 
104  if (n == 0)
105  return 0.;
106 
107  START_LOG ("median()", "StatisticsVector");
108 
109  std::sort(this->begin(), this->end());
110 
111  const dof_id_type lhs = (n-1) / 2;
112  const dof_id_type rhs = n / 2;
113 
114  Real the_median = 0;
115 
116 
117  if (lhs == rhs)
118  {
119  the_median = static_cast<Real>((*this)[lhs]);
120  }
121 
122  else
123  {
124  the_median = ( static_cast<Real>((*this)[lhs]) +
125  static_cast<Real>((*this)[rhs]) ) / 2.0;
126  }
127 
128  STOP_LOG ("median()", "StatisticsVector");
129 
130  return the_median;
131 }
132 
133 
134 
135 
136 template <typename T>
138 {
139  StatisticsVector<T> sv = (*this);
140 
141  return sv.median();
142 }
143 
144 
145 
146 
147 template <typename T>
149 {
150  const dof_id_type n = this->size();
151 
152  START_LOG ("variance()", "StatisticsVector");
153 
154  Real the_variance = 0;
155 
156  for (dof_id_type i=0; i<n; i++)
157  {
158  const Real delta = ( static_cast<Real>((*this)[i]) - mean_in );
159  the_variance += (delta * delta - the_variance) /
160  static_cast<Real>(i + 1);
161  }
162 
163  if (n > 1)
164  the_variance *= static_cast<Real>(n) / static_cast<Real>(n - 1);
165 
166  STOP_LOG ("variance()", "StatisticsVector");
167 
168  return the_variance;
169 }
170 
171 
172 template <typename T>
174 {
175  const dof_id_type n = this->size();
176  const Real max = this->maximum();
177 
178  for (dof_id_type i=0; i<n; i++)
179  {
180  (*this)[i] = static_cast<T>((*this)[i] / max);
181  }
182 }
183 
184 
185 
186 
187 
188 template <typename T>
189 void StatisticsVector<T>::histogram(std::vector<dof_id_type>& bin_members,
190  unsigned int n_bins)
191 {
192  // Must have at least 1 bin
193  libmesh_assert (n_bins>0);
194 
195  const dof_id_type n = this->size();
196 
197  std::sort(this->begin(), this->end());
198 
199  // The StatisticsVector can hold both integer and float types.
200  // We will define all the bins, etc. using Reals.
201  Real min = static_cast<Real>(this->minimum());
202  Real max = static_cast<Real>(this->maximum());
203  Real bin_size = (max - min) / static_cast<Real>(n_bins);
204 
205  START_LOG ("histogram()", "StatisticsVector");
206 
207  std::vector<Real> bin_bounds(n_bins+1);
208  for (unsigned int i=0; i<bin_bounds.size(); i++)
209  bin_bounds[i] = min + i * bin_size;
210 
211  // Give the last bin boundary a little wiggle room: we don't want
212  // it to be just barely less than the max, otherwise our bin test below
213  // may fail.
214  bin_bounds.back() += 1.e-6 * bin_size;
215 
216  // This vector will store the number of members each bin has.
217  bin_members.resize(n_bins);
218 
219  dof_id_type data_index = 0;
220  for (unsigned int j=0; j<bin_members.size(); j++) // bin vector indexing
221  {
222  // libMesh::out << "(debug) Filling bin " << j << std::endl;
223 
224  for (dof_id_type i=data_index; i<n; i++) // data vector indexing
225  {
226  // libMesh::out << "(debug) Processing index=" << i << std::endl;
227  Real current_val = static_cast<Real>( (*this)[i] );
228 
229  // There may be entries in the vector smaller than the value
230  // reported by this->minimum(). (e.g. inactive elements in an
231  // ErrorVector.) We just skip entries like that.
232  if ( current_val < min )
233  {
234  // libMesh::out << "(debug) Skipping entry v[" << i << "]="
235  // << (*this)[i]
236  // << " which is less than the min value: min="
237  // << min << std::endl;
238  continue;
239  }
240 
241  if ( current_val > bin_bounds[j+1] ) // if outside the current bin (bin[j] is bounded
242  // by bin_bounds[j] and bin_bounds[j+1])
243  {
244  // libMesh::out.precision(16);
245  // libMesh::out.setf(std::ios_base::fixed);
246  // libMesh::out << "(debug) (*this)[i]= " << (*this)[i]
247  // << " is greater than bin_bounds[j+1]="
248  // << bin_bounds[j+1] << std::endl;
249  data_index = i; // start searching here for next bin
250  break; // go to next bin
251  }
252 
253  // Otherwise, increment current bin's count
254  bin_members[j]++;
255  // libMesh::out << "(debug) Binned index=" << i << std::endl;
256  }
257  }
258 
259 #ifdef DEBUG
260  // Check the number of binned entries
261  const dof_id_type n_binned = std::accumulate(bin_members.begin(),
262  bin_members.end(),
263  static_cast<dof_id_type>(0),
264  std::plus<dof_id_type>());
265 
266  if (n != n_binned)
267  {
268  libMesh::out << "Warning: The number of binned entries, n_binned="
269  << n_binned
270  << ", did not match the total number of entries, n="
271  << n << "." << std::endl;
272  //libmesh_error();
273  }
274 #endif
275 
276 
277  STOP_LOG ("histogram()", "StatisticsVector");
278 }
279 
280 
281 
282 
283 
284 template <typename T>
286  const std::string& filename,
287  unsigned int n_bins)
288 {
289  // First generate the histogram with the desired number of bins
290  std::vector<dof_id_type> bin_members;
291  this->histogram(bin_members, n_bins);
292 
293  // The max, min and bin size are used to generate x-axis values.
294  T min = this->minimum();
295  T max = this->maximum();
296  T bin_size = (max - min) / static_cast<T>(n_bins);
297 
298  // On processor 0: Write histogram to file
299  if (my_procid==0)
300  {
301  std::ofstream out_stream (filename.c_str());
302 
303  out_stream << "clear all\n";
304  out_stream << "clf\n";
305  //out_stream << "x=linspace(" << min << "," << max << "," << n_bins+1 << ");\n";
306 
307  // abscissa values are located at the center of each bin.
308  out_stream << "x=[";
309  for (unsigned int i=0; i<bin_members.size(); ++i)
310  {
311  out_stream << min + (i+0.5)*bin_size << " ";
312  }
313  out_stream << "];\n";
314 
315  out_stream << "y=[";
316  for (unsigned int i=0; i<bin_members.size(); ++i)
317  {
318  out_stream << bin_members[i] << " ";
319  }
320  out_stream << "];\n";
321  out_stream << "bar(x,y);\n";
322  }
323 }
324 
325 
326 
327 template <typename T>
328 void StatisticsVector<T>::histogram(std::vector<dof_id_type>& bin_members,
329  unsigned int n_bins) const
330 {
331  StatisticsVector<T> sv = (*this);
332 
333  return sv.histogram(bin_members, n_bins);
334 }
335 
336 
337 
338 
339 template <typename T>
340 std::vector<dof_id_type> StatisticsVector<T>::cut_below(Real cut) const
341 {
342  START_LOG ("cut_below()", "StatisticsVector");
343 
344  const dof_id_type n = this->size();
345 
346  std::vector<dof_id_type> cut_indices;
347  cut_indices.reserve(n/2); // Arbitrary
348 
349  for (dof_id_type i=0; i<n; i++)
350  {
351  if ((*this)[i] < cut)
352  {
353  cut_indices.push_back(i);
354  }
355  }
356 
357  STOP_LOG ("cut_below()", "StatisticsVector");
358 
359  return cut_indices;
360 }
361 
362 
363 
364 
365 template <typename T>
366 std::vector<dof_id_type> StatisticsVector<T>::cut_above(Real cut) const
367 {
368  START_LOG ("cut_above()", "StatisticsVector");
369 
370  const dof_id_type n = this->size();
371 
372  std::vector<dof_id_type> cut_indices;
373  cut_indices.reserve(n/2); // Arbitrary
374 
375  for (dof_id_type i=0; i<n; i++)
376  {
377  if ((*this)[i] > cut)
378  {
379  cut_indices.push_back(i);
380  }
381  }
382 
383  STOP_LOG ("cut_above()", "StatisticsVector");
384 
385  return cut_indices;
386 }
387 
388 
389 
390 
391 //------------------------------------------------------------
392 // Explicit Instantions
393 template class StatisticsVector<float>;
394 template class StatisticsVector<double>;
395 #ifdef LIBMESH_DEFAULT_TRIPLE_PRECISION
396 template class StatisticsVector<long double>;
397 #endif
398 template class StatisticsVector<int>;
399 template class StatisticsVector<unsigned int>;
400 
401 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo