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