Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "Calculators.h"
11 :
12 : namespace StochasticTools
13 : {
14 :
15 : MultiMooseEnum
16 2564 : makeCalculatorEnum()
17 : {
18 : return MultiMooseEnum(
19 5128 : "min=0 max=1 sum=2 mean=3 stddev=4 norm2=5 ratio=6 stderr=7 median=8 meanabs=9");
20 : }
21 :
22 : // MEAN ////////////////////////////////////////////////////////////////////////////////////////////
23 : template <typename InType, typename OutType>
24 : void
25 13060723 : Mean<InType, OutType>::initialize()
26 : {
27 13060723 : _count = 0;
28 13060723 : _sum = OutType();
29 13060723 : }
30 :
31 : template <typename InType, typename OutType>
32 : void
33 897922742 : Mean<InType, OutType>::update(const typename InType::value_type & val)
34 : {
35 1025730422 : _count++;
36 1025730422 : _sum += static_cast<OutType>(val);
37 897922742 : }
38 :
39 : template <typename InType, typename OutType>
40 : void
41 12856695 : Mean<InType, OutType>::finalize(bool is_distributed)
42 : {
43 12856695 : if (is_distributed)
44 : {
45 3856116 : this->_communicator.sum(_count);
46 3856116 : this->_communicator.sum(_sum);
47 : }
48 12856695 : if (_count > 0)
49 12856695 : _sum /= static_cast<OutType>(_count);
50 12856695 : }
51 :
52 : // MEAN ABS ////////////////////////////////////////////////////////////////////////////////////////
53 : template <typename InType, typename OutType>
54 : void
55 127807680 : MeanAbsoluteValue<InType, OutType>::update(const typename InType::value_type & val)
56 : {
57 127807680 : Mean<InType, OutType>::update(std::abs(val));
58 127807680 : }
59 :
60 : // SUM /////////////////////////////////////////////////////////////////////////////////////////////
61 : template <typename InType, typename OutType>
62 : void
63 204028 : Sum<InType, OutType>::finalize(bool is_distributed)
64 : {
65 204028 : if (is_distributed)
66 3940 : this->_communicator.sum(this->_sum);
67 204028 : }
68 :
69 : // STDDEV //////////////////////////////////////////////////////////////////////////////////////////
70 : template <typename InType, typename OutType>
71 : void
72 11131721 : StdDev<InType, OutType>::initialize()
73 : {
74 11131721 : _count = 0;
75 11131721 : _sum = OutType();
76 11131721 : _sum_of_square = OutType();
77 11131721 : }
78 :
79 : template <typename InType, typename OutType>
80 : void
81 873509282 : StdDev<InType, OutType>::update(const typename InType::value_type & val)
82 : {
83 873509282 : _count++;
84 873509282 : _sum += static_cast<OutType>(val);
85 873509282 : _sum_of_square += MathUtils::pow(static_cast<OutType>(val), 2);
86 873509282 : }
87 :
88 : template <typename InType, typename OutType>
89 : void
90 11131721 : StdDev<InType, OutType>::finalize(bool is_distributed)
91 : {
92 11131721 : if (is_distributed)
93 : {
94 2573314 : this->_communicator.sum(_count);
95 2573314 : this->_communicator.sum(_sum);
96 2573314 : this->_communicator.sum(_sum_of_square);
97 : }
98 :
99 11131721 : if (_count <= 1)
100 0 : _sum_of_square = 0;
101 : else
102 11131721 : _sum_of_square = std::sqrt(std::abs(_sum_of_square - _sum * _sum / _count) / (_count - 1));
103 11131721 : }
104 :
105 : // STDERR //////////////////////////////////////////////////////////////////////////////////////////
106 : template <typename InType, typename OutType>
107 : void
108 240282 : StdErr<InType, OutType>::finalize(bool is_distributed)
109 : {
110 240282 : StdDev<InType, OutType>::finalize(is_distributed);
111 240282 : this->_sum_of_square /= std::sqrt(this->_count);
112 240282 : }
113 :
114 : // RATIO ///////////////////////////////////////////////////////////////////////////////////////////
115 : template <typename InType, typename OutType>
116 : void
117 600856 : Ratio<InType, OutType>::initialize()
118 : {
119 600856 : _min = std::numeric_limits<OutType>::max();
120 600856 : _max = std::numeric_limits<OutType>::min();
121 600856 : }
122 :
123 : template <typename InType, typename OutType>
124 : void
125 3004748 : Ratio<InType, OutType>::update(const typename InType::value_type & val)
126 : {
127 3004748 : if (_min > val)
128 1144834 : _min = static_cast<OutType>(val);
129 3004748 : if (_max < val)
130 1154594 : _max = static_cast<OutType>(val);
131 3004748 : }
132 :
133 : template <typename InType, typename OutType>
134 : void
135 600856 : Ratio<InType, OutType>::finalize(bool is_distributed)
136 : {
137 600856 : if (is_distributed)
138 : {
139 612 : this->_communicator.min(_min);
140 612 : this->_communicator.max(_max);
141 : }
142 600856 : }
143 :
144 : // L2NORM //////////////////////////////////////////////////////////////////////////////////////////
145 : template <typename InType, typename OutType>
146 : void
147 200292 : L2Norm<InType, OutType>::initialize()
148 : {
149 200292 : _l2_norm = OutType();
150 200292 : }
151 :
152 : template <typename InType, typename OutType>
153 : void
154 1001616 : L2Norm<InType, OutType>::update(const typename InType::value_type & val)
155 : {
156 1001616 : _l2_norm += MathUtils::pow(static_cast<OutType>(val), 2);
157 1001616 : }
158 :
159 : template <typename InType, typename OutType>
160 : void
161 200292 : L2Norm<InType, OutType>::finalize(bool is_distributed)
162 : {
163 200292 : if (is_distributed)
164 204 : this->_communicator.sum(_l2_norm);
165 200292 : _l2_norm = std::sqrt(_l2_norm);
166 200292 : }
167 :
168 : // MEDIAN //////////////////////////////////////////////////////////////////////////////////////////
169 : template <typename InType, typename OutType>
170 : void
171 200182 : Median<InType, OutType>::initialize()
172 : {
173 : _storage.clear();
174 200182 : }
175 :
176 : template <typename InType, typename OutType>
177 : void
178 1000766 : Median<InType, OutType>::update(const typename InType::value_type & val)
179 : {
180 1000766 : _storage.push_back(static_cast<OutType>(val));
181 1000766 : }
182 :
183 : template <typename InType, typename OutType>
184 : void
185 200182 : Median<InType, OutType>::finalize(bool is_distributed)
186 : {
187 : // Make sure we aren't doing anything silly like taking the median of an empty vector
188 200182 : _median = OutType();
189 200182 : auto count = _storage.size();
190 200182 : if (is_distributed)
191 144 : this->_communicator.sum(count);
192 200182 : if (count == 0)
193 200074 : return;
194 :
195 200182 : if (!is_distributed || this->n_processors() == 1)
196 : {
197 200074 : std::sort(_storage.begin(), _storage.end());
198 200074 : if (count % 2)
199 200038 : _median = _storage[count / 2];
200 : else
201 36 : _median += (_storage[count / 2] + _storage[count / 2 - 1]) / 2;
202 200074 : return;
203 : }
204 :
205 108 : dof_id_type kgt = count % 2 ? (count / 2) : (count / 2 - 1);
206 : dof_id_type klt = kgt;
207 : while (true)
208 : {
209 : // Gather all sizes and figure out current number of values
210 324 : std::vector<std::size_t> sz = {_storage.size()};
211 324 : this->_communicator.allgather(sz);
212 324 : dof_id_type n = std::accumulate(sz.begin(), sz.end(), 0);
213 :
214 : // Choose the first value for the first processor with values
215 324 : for (const auto & i : index_range(sz))
216 324 : if (sz[i])
217 : {
218 324 : if (this->processor_id() == i)
219 162 : _median = _storage[0];
220 324 : this->_communicator.broadcast(_median, i);
221 : break;
222 : }
223 :
224 : // Count number of values greater than, less than, and equal to _median
225 324 : std::vector<dof_id_type> m(3, 0);
226 1134 : for (const auto & val : _storage)
227 : {
228 810 : if (_median < val)
229 648 : m[0]++;
230 162 : else if (val < _median)
231 0 : m[1]++;
232 : }
233 324 : this->_communicator.sum(m);
234 324 : m[2] = n - m[0] - m[1];
235 :
236 : // Remove greater than equal to
237 324 : if ((m[0] + m[2]) <= kgt)
238 : {
239 0 : _storage.erase(std::remove_if(_storage.begin(),
240 : _storage.end(),
241 0 : [this](const OutType & val) { return val >= _median; }),
242 : _storage.end());
243 0 : kgt -= m[0] + m[2];
244 : }
245 : // Remove less than equal to
246 324 : else if ((m[1] + m[2]) <= klt)
247 : {
248 216 : _storage.erase(std::remove_if(_storage.begin(),
249 : _storage.end(),
250 594 : [this](const OutType & val) { return val <= _median; }),
251 : _storage.end());
252 216 : klt -= m[1] + m[2];
253 : }
254 : // If the number of points is odd, then we've found it
255 108 : else if (count % 2)
256 : break;
257 : // Get average of the two middle numbers
258 : else
259 : {
260 : OutType num2;
261 : // Find the next greater than
262 108 : if (m[0] > kgt)
263 : {
264 108 : num2 = std::numeric_limits<OutType>::max();
265 324 : for (const auto & val : _storage)
266 216 : if (_median < val && val < num2)
267 54 : num2 = val;
268 108 : this->_communicator.min(num2);
269 : }
270 : // Find the next less than
271 0 : else if (m[1] > klt)
272 : {
273 0 : num2 = std::numeric_limits<OutType>::min();
274 0 : for (const auto & val : _storage)
275 0 : if (val < _median && num2 < val)
276 0 : num2 += val;
277 0 : this->_communicator.max(num2);
278 : }
279 : // Otherwise we know the other number is equal
280 : else
281 0 : num2 = _median;
282 :
283 108 : _median = (_median + num2) / 2;
284 : break;
285 : }
286 : }
287 : }
288 :
289 : // CalculatorBuilder
290 : // //////////////////////////////////////////////////////////////////////////////////
291 : template <typename InType, typename OutType>
292 : std::unique_ptr<Calculator<InType, OutType>>
293 11642 : CalculatorBuilder<InType, OutType>::build(const MooseEnumItem & item,
294 : const libMesh::ParallelObject & other)
295 : {
296 11642 : if (item == "min")
297 200 : return std::make_unique<Min<InType, OutType>>(other, item);
298 :
299 11442 : else if (item == "max")
300 200 : return std::make_unique<Max<InType, OutType>>(other, item);
301 :
302 11242 : else if (item == "sum")
303 2728 : return std::make_unique<Sum<InType, OutType>>(other, item);
304 :
305 8514 : else if (item == "mean" || item == "average") // average is deprecated
306 4108 : return std::make_unique<Mean<InType, OutType>>(other, item);
307 :
308 4406 : else if (item == "stddev")
309 3674 : return std::make_unique<StdDev<InType, OutType>>(other, item);
310 :
311 732 : else if (item == "stderr")
312 196 : return std::make_unique<StdErr<InType, OutType>>(other, item);
313 :
314 536 : else if (item == "norm2")
315 200 : return std::make_unique<L2Norm<InType, OutType>>(other, item);
316 :
317 336 : else if (item == "ratio")
318 174 : return std::make_unique<Ratio<InType, OutType>>(other, item);
319 :
320 162 : else if (item == "median")
321 66 : return std::make_unique<Median<InType, OutType>>(other, item);
322 :
323 96 : else if (item == "meanabs")
324 96 : return std::make_unique<MeanAbsoluteValue<InType, OutType>>(other, item);
325 :
326 0 : ::mooseError("Failed to create Statistics::Calculator object for ", item);
327 : return nullptr;
328 : }
329 :
330 : #define createCalculators(InType, OutType) \
331 : template class Mean<InType, OutType>; \
332 : template class Max<InType, OutType>; \
333 : template class Min<InType, OutType>; \
334 : template class Sum<InType, OutType>; \
335 : template class StdDev<InType, OutType>; \
336 : template class StdErr<InType, OutType>; \
337 : template class Ratio<InType, OutType>; \
338 : template class L2Norm<InType, OutType>; \
339 : template class Median<InType, OutType>; \
340 : template struct CalculatorBuilder<InType, OutType>
341 :
342 : createCalculators(std::vector<Real>, Real);
343 : createCalculators(std::vector<int>, Real);
344 :
345 : } // StocasticTools namespace
|
