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 1119 : makeCalculatorEnum()
17 : {
18 : return MultiMooseEnum(
19 2238 : "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 6554462 : Mean<InType, OutType>::initialize()
26 : {
27 6554462 : _count = 0;
28 6554462 : _sum = OutType();
29 6554462 : }
30 :
31 : template <typename InType, typename OutType>
32 : void
33 732239646 : Mean<InType, OutType>::update(const typename InType::value_type & val)
34 : {
35 796143486 : _count++;
36 796143486 : _sum += static_cast<OutType>(val);
37 732239646 : }
38 :
39 : template <typename InType, typename OutType>
40 : void
41 6452338 : Mean<InType, OutType>::finalize(bool is_distributed)
42 : {
43 6452338 : if (is_distributed)
44 : {
45 1687518 : this->_communicator.sum(_count);
46 1687518 : this->_communicator.sum(_sum);
47 : }
48 6452338 : if (_count > 0)
49 6452338 : _sum /= static_cast<OutType>(_count);
50 6452338 : }
51 :
52 : // MEAN ABS ////////////////////////////////////////////////////////////////////////////////////////
53 : template <typename InType, typename OutType>
54 : void
55 63903840 : MeanAbsoluteValue<InType, OutType>::update(const typename InType::value_type & val)
56 : {
57 63903840 : Mean<InType, OutType>::update(std::abs(val));
58 63903840 : }
59 :
60 : // SUM /////////////////////////////////////////////////////////////////////////////////////////////
61 : template <typename InType, typename OutType>
62 : void
63 102124 : Sum<InType, OutType>::finalize(bool is_distributed)
64 : {
65 102124 : if (is_distributed)
66 2071 : this->_communicator.sum(this->_sum);
67 102124 : }
68 :
69 : // STDDEV //////////////////////////////////////////////////////////////////////////////////////////
70 : template <typename InType, typename OutType>
71 : void
72 5669917 : StdDev<InType, OutType>::initialize()
73 : {
74 5669917 : _count = 0;
75 5669917 : _sum = OutType();
76 5669917 : _sum_of_square = OutType();
77 5669917 : }
78 :
79 : template <typename InType, typename OutType>
80 : void
81 720031716 : StdDev<InType, OutType>::update(const typename InType::value_type & val)
82 : {
83 720031716 : _count++;
84 720031716 : _sum += static_cast<OutType>(val);
85 720031716 : _sum_of_square += MathUtils::pow(static_cast<OutType>(val), 2);
86 720031716 : }
87 :
88 : template <typename InType, typename OutType>
89 : void
90 5669917 : StdDev<InType, OutType>::finalize(bool is_distributed)
91 : {
92 5669917 : if (is_distributed)
93 : {
94 1126174 : this->_communicator.sum(_count);
95 1126174 : this->_communicator.sum(_sum);
96 1126174 : this->_communicator.sum(_sum_of_square);
97 : }
98 :
99 5669917 : if (_count <= 1)
100 0 : _sum_of_square = 0;
101 : else
102 5669917 : _sum_of_square = std::sqrt(std::abs(_sum_of_square - _sum * _sum / _count) / (_count - 1));
103 5669917 : }
104 :
105 : // STDERR //////////////////////////////////////////////////////////////////////////////////////////
106 : template <typename InType, typename OutType>
107 : void
108 120141 : StdErr<InType, OutType>::finalize(bool is_distributed)
109 : {
110 120141 : StdDev<InType, OutType>::finalize(is_distributed);
111 120141 : this->_sum_of_square /= std::sqrt(this->_count);
112 120141 : }
113 :
114 : // RATIO ///////////////////////////////////////////////////////////////////////////////////////////
115 : template <typename InType, typename OutType>
116 : void
117 300428 : Ratio<InType, OutType>::initialize()
118 : {
119 300428 : _min = std::numeric_limits<OutType>::max();
120 300428 : _max = std::numeric_limits<OutType>::min();
121 300428 : }
122 :
123 : template <typename InType, typename OutType>
124 : void
125 1502563 : Ratio<InType, OutType>::update(const typename InType::value_type & val)
126 : {
127 1502563 : if (_min > val)
128 572456 : _min = static_cast<OutType>(val);
129 1502563 : if (_max < val)
130 577369 : _max = static_cast<OutType>(val);
131 1502563 : }
132 :
133 : template <typename InType, typename OutType>
134 : void
135 300428 : Ratio<InType, OutType>::finalize(bool is_distributed)
136 : {
137 300428 : if (is_distributed)
138 : {
139 279 : this->_communicator.min(_min);
140 279 : this->_communicator.max(_max);
141 : }
142 300428 : }
143 :
144 : // L2NORM //////////////////////////////////////////////////////////////////////////////////////////
145 : template <typename InType, typename OutType>
146 : void
147 100146 : L2Norm<InType, OutType>::initialize()
148 : {
149 100146 : _l2_norm = OutType();
150 100146 : }
151 :
152 : template <typename InType, typename OutType>
153 : void
154 500871 : L2Norm<InType, OutType>::update(const typename InType::value_type & val)
155 : {
156 500871 : _l2_norm += MathUtils::pow(static_cast<OutType>(val), 2);
157 500871 : }
158 :
159 : template <typename InType, typename OutType>
160 : void
161 100146 : L2Norm<InType, OutType>::finalize(bool is_distributed)
162 : {
163 100146 : if (is_distributed)
164 93 : this->_communicator.sum(_l2_norm);
165 100146 : _l2_norm = std::sqrt(_l2_norm);
166 100146 : }
167 :
168 : // MEDIAN //////////////////////////////////////////////////////////////////////////////////////////
169 : template <typename InType, typename OutType>
170 : void
171 100091 : Median<InType, OutType>::initialize()
172 : {
173 100091 : _storage.clear();
174 100091 : }
175 :
176 : template <typename InType, typename OutType>
177 : void
178 500446 : Median<InType, OutType>::update(const typename InType::value_type & val)
179 : {
180 500446 : _storage.push_back(static_cast<OutType>(val));
181 500446 : }
182 :
183 : template <typename InType, typename OutType>
184 : void
185 100091 : 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 100091 : _median = OutType();
189 100091 : auto count = _storage.size();
190 100091 : if (is_distributed)
191 63 : this->_communicator.sum(count);
192 100091 : if (count == 0)
193 100055 : return;
194 :
195 100091 : if (!is_distributed || this->n_processors() == 1)
196 : {
197 100055 : std::sort(_storage.begin(), _storage.end());
198 100055 : if (count % 2)
199 100028 : _median = _storage[count / 2];
200 : else
201 27 : _median += (_storage[count / 2] + _storage[count / 2 - 1]) / 2;
202 100055 : return;
203 : }
204 :
205 36 : dof_id_type kgt = count % 2 ? (count / 2) : (count / 2 - 1);
206 : dof_id_type klt = kgt;
207 108 : while (true)
208 : {
209 : // Gather all sizes and figure out current number of values
210 108 : std::vector<std::size_t> sz = {_storage.size()};
211 108 : this->_communicator.allgather(sz);
212 108 : dof_id_type n = std::accumulate(sz.begin(), sz.end(), 0);
213 :
214 : // Choose the first value for the first processor with values
215 108 : for (const auto & i : index_range(sz))
216 108 : if (sz[i])
217 : {
218 108 : if (this->processor_id() == i)
219 54 : _median = _storage[0];
220 108 : this->_communicator.broadcast(_median, i);
221 : break;
222 : }
223 :
224 : // Count number of values greater than, less than, and equal to _median
225 108 : std::vector<dof_id_type> m(3, 0);
226 378 : for (const auto & val : _storage)
227 : {
228 270 : if (_median < val)
229 216 : m[0]++;
230 54 : else if (val < _median)
231 0 : m[1]++;
232 : }
233 108 : this->_communicator.sum(m);
234 108 : m[2] = n - m[0] - m[1];
235 :
236 : // Remove greater than equal to
237 108 : 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 108 : else if ((m[1] + m[2]) <= klt)
247 : {
248 72 : _storage.erase(std::remove_if(_storage.begin(),
249 : _storage.end(),
250 198 : [this](const OutType & val) { return val <= _median; }),
251 : _storage.end());
252 72 : klt -= m[1] + m[2];
253 : }
254 : // If the number of points is odd, then we've found it
255 36 : 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 36 : if (m[0] > kgt)
263 : {
264 36 : num2 = std::numeric_limits<OutType>::max();
265 108 : for (const auto & val : _storage)
266 72 : if (_median < val && val < num2)
267 18 : num2 = val;
268 36 : 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 36 : _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 5288 : CalculatorBuilder<InType, OutType>::build(const MooseEnumItem & item,
294 : const libMesh::ParallelObject & other)
295 : {
296 5288 : if (item == "min")
297 93 : return std::make_unique<Min<InType, OutType>>(other, item);
298 :
299 5195 : else if (item == "max")
300 93 : return std::make_unique<Max<InType, OutType>>(other, item);
301 :
302 5102 : else if (item == "sum")
303 1257 : return std::make_unique<Sum<InType, OutType>>(other, item);
304 :
305 3845 : else if (item == "mean" || item == "average") // average is deprecated
306 1856 : return std::make_unique<Mean<InType, OutType>>(other, item);
307 :
308 1989 : else if (item == "stddev")
309 1653 : return std::make_unique<StdDev<InType, OutType>>(other, item);
310 :
311 336 : else if (item == "stderr")
312 90 : return std::make_unique<StdErr<InType, OutType>>(other, item);
313 :
314 246 : else if (item == "norm2")
315 93 : return std::make_unique<L2Norm<InType, OutType>>(other, item);
316 :
317 153 : else if (item == "ratio")
318 81 : return std::make_unique<Ratio<InType, OutType>>(other, item);
319 :
320 72 : else if (item == "median")
321 30 : return std::make_unique<Median<InType, OutType>>(other, item);
322 :
323 42 : else if (item == "meanabs")
324 42 : 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
|
