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 "VectorCalculators.h"
11 :
12 : namespace StochasticTools
13 : {
14 :
15 : template <typename InType, typename OutType>
16 : std::vector<std::vector<OutType>>
17 843 : Percentile<std::vector<InType>, std::vector<OutType>>::compute(const std::vector<InType> & data,
18 : const bool is_distributed)
19 : {
20 : // Bootstrap estimates
21 843 : const auto values = this->computeBootstrapEstimates(data, is_distributed);
22 :
23 : // Extract percentiles
24 : std::vector<std::vector<OutType>> output;
25 843 : if (this->processor_id() == 0)
26 2007 : for (const Real & level : this->_levels)
27 : {
28 1434 : long unsigned int index = std::lrint(level * (this->_replicates - 1));
29 1434 : output.push_back(values[index]);
30 : }
31 :
32 843 : return output;
33 843 : }
34 :
35 : template <typename InType, typename OutType>
36 : std::vector<std::vector<OutType>>
37 34 : BiasCorrectedAccelerated<std::vector<InType>, std::vector<OutType>>::compute(
38 : const std::vector<InType> & data, bool is_distributed)
39 : {
40 : // Bootstrap estimates
41 34 : const std::vector<std::vector<OutType>> values =
42 : this->computeBootstrapEstimates(data, is_distributed);
43 :
44 : // Compute bias-correction, Efron and Tibshirani (2003), Eq. 14.14, p. 186
45 34 : const std::vector<OutType> value = this->_calc.compute(data, is_distributed);
46 34 : std::vector<Real> bias(value.size());
47 246 : for (const auto & i : index_range(value))
48 : {
49 : Real count = 0.;
50 1095517 : for (const auto & val : values)
51 : {
52 1095517 : if (val[i] < value[i])
53 1095305 : count++;
54 : else
55 : break;
56 : }
57 212 : bias[i] = NormalDistribution::quantile(count / this->_replicates, 0, 1);
58 : }
59 :
60 : // Compute Acceleration, Efron and Tibshirani (2003), Eq. 14.15, p. 186
61 34 : const std::vector<OutType> acc =
62 34 : data.empty() ? std::vector<OutType>(value.size()) : acceleration(data, is_distributed);
63 :
64 : // Compute intervals, Efron and Tibshirani (2003), Eq. 14.10, p. 185
65 34 : std::vector<std::vector<OutType>> output(this->_levels.size(),
66 68 : std::vector<OutType>(value.size()));
67 334 : for (const auto & l : index_range(this->_levels))
68 : {
69 300 : const Real z = NormalDistribution::quantile(this->_levels[l], 0, 1);
70 2052 : for (const auto & i : index_range(value))
71 : {
72 1752 : const OutType x = bias[i] + (bias[i] + (bias[i] + z) / (1 - acc[i] * (bias[i] + z)));
73 1752 : const Real alpha = NormalDistribution::cdf(x, 0, 1);
74 :
75 1752 : long unsigned int index = std::lrint(alpha * (this->_replicates - 1));
76 1752 : output[l][i] = values[index][i];
77 : }
78 : }
79 34 : return output;
80 34 : }
81 :
82 : template <typename InType, typename OutType>
83 : std::vector<OutType>
84 34 : BiasCorrectedAccelerated<std::vector<InType>, std::vector<OutType>>::acceleration(
85 : const std::vector<InType> & data, const bool is_distributed)
86 : {
87 : const std::size_t local_size = data.size();
88 34 : std::vector<std::size_t> local_sizes = {local_size};
89 34 : if (is_distributed)
90 32 : this->_communicator.allgather(local_sizes);
91 34 : const std::size_t count = std::accumulate(local_sizes.begin(), local_sizes.end(), 0);
92 34 : const processor_id_type rank = is_distributed ? this->processor_id() : 0;
93 :
94 : // Jackknife statistics
95 34 : std::vector<std::vector<OutType>> theta_i(local_size);
96 :
97 : // Compute jackknife estimates, Ch. 11, Eq. 11.2, p. 141
98 92 : for (processor_id_type r = 0; r < local_sizes.size(); ++r)
99 5258 : for (std::size_t i = 0; i < local_sizes[r]; ++i)
100 : {
101 5200 : this->_calc.initializeCalculator();
102 2205200 : for (std::size_t il = 0; il < local_size; ++il)
103 2200000 : if (i != il || r != rank)
104 2196000 : this->_calc.updateCalculator(data[il]);
105 5200 : this->_calc.finalizeCalculator(is_distributed);
106 5200 : if (r == rank)
107 8000 : theta_i[i] = this->_calc.getValue();
108 : }
109 :
110 : // Compute jackknife sum, Ch. 11, Eq. 11.4, p. 141
111 34 : const std::size_t nval = theta_i.size() > 0 ? theta_i[0].size() : 0;
112 34 : std::vector<OutType> theta_dot(nval, 0.);
113 4034 : for (const auto & ti : theta_i)
114 66000 : for (const auto & i : make_range(nval))
115 62000 : theta_dot[i] += ti[i] / count;
116 34 : if (is_distributed)
117 32 : this->_communicator.sum(theta_dot);
118 :
119 : // Acceleration, Ch. 14, Eq. 14.15, p. 185
120 34 : std::vector<OutType> num_den(2 * nval);
121 4034 : for (const auto & jk : theta_i)
122 66000 : for (const auto & i : make_range(nval))
123 : {
124 62000 : num_den[i] += MathUtils::pow(theta_dot[i] - jk[i], 3);
125 62000 : num_den[nval + i] += MathUtils::pow(theta_dot[i] - jk[i], 2);
126 : }
127 34 : if (is_distributed)
128 32 : this->_communicator.sum(num_den);
129 :
130 : mooseAssert(std::find(num_den.begin() + nval, num_den.end(), OutType()) == num_den.end(),
131 : "The acceleration denomenator must not be zero.");
132 34 : std::vector<OutType> acc(nval);
133 246 : for (const auto & i : make_range(nval))
134 212 : acc[i] = num_den[i] / (6 * std::pow(num_den[nval + i], 3. / 2.));
135 34 : return acc;
136 34 : }
137 :
138 : template <typename InType, typename OutType>
139 : std::unique_ptr<Calculator<std::vector<InType>, std::vector<OutType>>>
140 413 : CalculatorBuilder<std::vector<InType>, std::vector<OutType>>::build(
141 : const MooseEnumItem & item, const libMesh::ParallelObject & other)
142 : {
143 413 : if (item == "min")
144 17 : return std::make_unique<VectorCalculator<InType, OutType, Min>>(other, item);
145 :
146 396 : else if (item == "max")
147 17 : return std::make_unique<VectorCalculator<InType, OutType, Max>>(other, item);
148 :
149 379 : else if (item == "sum")
150 17 : return std::make_unique<VectorCalculator<InType, OutType, Sum>>(other, item);
151 :
152 362 : else if (item == "mean" || item == "average") // average is deprecated
153 131 : return std::make_unique<VectorCalculator<InType, OutType, Mean>>(other, item);
154 :
155 231 : else if (item == "stddev")
156 131 : return std::make_unique<VectorCalculator<InType, OutType, StdDev>>(other, item);
157 :
158 100 : else if (item == "stderr")
159 17 : return std::make_unique<VectorCalculator<InType, OutType, StdErr>>(other, item);
160 :
161 83 : else if (item == "norm2")
162 17 : return std::make_unique<VectorCalculator<InType, OutType, L2Norm>>(other, item);
163 :
164 66 : else if (item == "ratio")
165 17 : return std::make_unique<VectorCalculator<InType, OutType, Ratio>>(other, item);
166 :
167 49 : else if (item == "median")
168 17 : return std::make_unique<VectorCalculator<InType, OutType, Median>>(other, item);
169 :
170 32 : else if (item == "meanabs")
171 32 : return std::make_unique<VectorCalculator<InType, OutType, MeanAbsoluteValue>>(other, item);
172 :
173 0 : ::mooseError("Failed to create Statistics::Calculator object for ", item);
174 : return nullptr;
175 : }
176 :
177 : template <typename InType, typename OutType>
178 : std::unique_ptr<BootstrapCalculator<std::vector<InType>, std::vector<OutType>>>
179 373 : BootstrapCalculatorBuilder<std::vector<InType>, std::vector<OutType>>::build(
180 : const MooseEnum & item,
181 : const libMesh::ParallelObject & other,
182 : const std::vector<Real> & levels,
183 : unsigned int replicates,
184 : unsigned int seed,
185 : StochasticTools::Calculator<std::vector<InType>, std::vector<OutType>> & calc)
186 : {
187 : std::unique_ptr<BootstrapCalculator<std::vector<InType>, std::vector<OutType>>> ptr = nullptr;
188 373 : if (item == "percentile")
189 339 : ptr = std::make_unique<Percentile<std::vector<InType>, std::vector<OutType>>>(
190 : other, item, levels, replicates, seed, calc);
191 34 : else if (item == "bca")
192 34 : ptr = std::make_unique<BiasCorrectedAccelerated<std::vector<InType>, std::vector<OutType>>>(
193 : other, item, levels, replicates, seed, calc);
194 : else
195 0 : ::mooseError("Failed to create Statistics::BootstrapCalculator object for ", item);
196 :
197 373 : return ptr;
198 0 : }
199 :
200 : #define createVectorCalculators(InType, OutType) \
201 : template class Percentile<std::vector<InType>, std::vector<OutType>>; \
202 : template class BiasCorrectedAccelerated<std::vector<InType>, std::vector<OutType>>; \
203 : template struct CalculatorBuilder<std::vector<InType>, std::vector<OutType>>; \
204 : template struct BootstrapCalculatorBuilder<std::vector<InType>, std::vector<OutType>>
205 :
206 : createVectorCalculators(std::vector<Real>, Real);
207 :
208 : } // StocasticTools namespace
|
