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 378 : Percentile<std::vector<InType>, std::vector<OutType>>::compute(const std::vector<InType> & data,
18 : const bool is_distributed)
19 : {
20 : // Bootstrap estimates
21 378 : const auto values = this->computeBootstrapEstimates(data, is_distributed);
22 :
23 : // Extract percentiles
24 : std::vector<std::vector<OutType>> output;
25 378 : if (this->processor_id() == 0)
26 1012 : for (const Real & level : this->_levels)
27 : {
28 724 : long unsigned int index = std::lrint(level * (this->_replicates - 1));
29 724 : output.push_back(values[index]);
30 : }
31 :
32 378 : return output;
33 378 : }
34 :
35 : template <typename InType, typename OutType>
36 : std::vector<std::vector<OutType>>
37 16 : BiasCorrectedAccelerated<std::vector<InType>, std::vector<OutType>>::compute(
38 : const std::vector<InType> & data, bool is_distributed)
39 : {
40 : // Bootstrap estimates
41 16 : 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 16 : const std::vector<OutType> value = this->_calc.compute(data, is_distributed);
46 16 : std::vector<Real> bias(value.size());
47 138 : for (const auto & i : index_range(value))
48 : {
49 : Real count = 0.;
50 627940 : for (const auto & val : values)
51 : {
52 627940 : if (val[i] < value[i])
53 627818 : count++;
54 : else
55 : break;
56 : }
57 122 : bias[i] = NormalDistribution::quantile(count / this->_replicates, 0, 1);
58 : }
59 :
60 : // Compute Acceleration, Efron and Tibshirani (2003), Eq. 14.15, p. 186
61 16 : const std::vector<OutType> acc =
62 16 : 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 16 : std::vector<std::vector<OutType>> output(this->_levels.size(),
66 32 : std::vector<OutType>(value.size()));
67 154 : for (const auto & l : index_range(this->_levels))
68 : {
69 138 : const Real z = NormalDistribution::quantile(this->_levels[l], 0, 1);
70 1080 : for (const auto & i : index_range(value))
71 : {
72 942 : const OutType x = bias[i] + (bias[i] + (bias[i] + z) / (1 - acc[i] * (bias[i] + z)));
73 942 : const Real alpha = NormalDistribution::cdf(x, 0, 1);
74 :
75 942 : long unsigned int index = std::lrint(alpha * (this->_replicates - 1));
76 942 : output[l][i] = values[index][i];
77 : }
78 : }
79 16 : return output;
80 16 : }
81 :
82 : template <typename InType, typename OutType>
83 : std::vector<OutType>
84 16 : 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 16 : std::vector<std::size_t> local_sizes = {local_size};
89 16 : if (is_distributed)
90 14 : this->_communicator.allgather(local_sizes);
91 16 : const std::size_t count = std::accumulate(local_sizes.begin(), local_sizes.end(), 0);
92 16 : const processor_id_type rank = is_distributed ? this->processor_id() : 0;
93 :
94 : // Jackknife statistics
95 16 : std::vector<std::vector<OutType>> theta_i(local_size);
96 :
97 : // Compute jackknife estimates, Ch. 11, Eq. 11.2, p. 141
98 40 : for (processor_id_type r = 0; r < local_sizes.size(); ++r)
99 3424 : for (std::size_t i = 0; i < local_sizes[r]; ++i)
100 : {
101 3400 : this->_calc.initializeCalculator();
102 2103400 : for (std::size_t il = 0; il < local_size; ++il)
103 2100000 : if (i != il || r != rank)
104 2097000 : this->_calc.updateCalculator(data[il]);
105 3400 : this->_calc.finalizeCalculator(is_distributed);
106 3400 : if (r == rank)
107 6000 : theta_i[i] = this->_calc.getValue();
108 : }
109 :
110 : // Compute jackknife sum, Ch. 11, Eq. 11.4, p. 141
111 16 : const std::size_t nval = theta_i.size() > 0 ? theta_i[0].size() : 0;
112 16 : std::vector<OutType> theta_dot(nval, 0.);
113 3016 : for (const auto & ti : theta_i)
114 60000 : for (const auto & i : make_range(nval))
115 57000 : theta_dot[i] += ti[i] / count;
116 16 : if (is_distributed)
117 14 : this->_communicator.sum(theta_dot);
118 :
119 : // Acceleration, Ch. 14, Eq. 14.15, p. 185
120 16 : std::vector<OutType> num_den(2 * nval);
121 3016 : for (const auto & jk : theta_i)
122 60000 : for (const auto & i : make_range(nval))
123 : {
124 57000 : num_den[i] += MathUtils::pow(theta_dot[i] - jk[i], 3);
125 57000 : num_den[nval + i] += MathUtils::pow(theta_dot[i] - jk[i], 2);
126 : }
127 16 : if (is_distributed)
128 14 : 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 16 : std::vector<OutType> acc(nval);
133 138 : for (const auto & i : make_range(nval))
134 122 : acc[i] = num_den[i] / (6 * std::pow(num_den[nval + i], 3. / 2.));
135 16 : return acc;
136 16 : }
137 :
138 : template <typename InType, typename OutType>
139 : std::unique_ptr<Calculator<std::vector<InType>, std::vector<OutType>>>
140 188 : CalculatorBuilder<std::vector<InType>, std::vector<OutType>>::build(
141 : const MooseEnumItem & item, const libMesh::ParallelObject & other)
142 : {
143 188 : if (item == "min")
144 8 : return std::make_unique<VectorCalculator<InType, OutType, Min>>(other, item);
145 :
146 180 : else if (item == "max")
147 8 : return std::make_unique<VectorCalculator<InType, OutType, Max>>(other, item);
148 :
149 172 : else if (item == "sum")
150 8 : return std::make_unique<VectorCalculator<InType, OutType, Sum>>(other, item);
151 :
152 164 : else if (item == "mean" || item == "average") // average is deprecated
153 59 : return std::make_unique<VectorCalculator<InType, OutType, Mean>>(other, item);
154 :
155 105 : else if (item == "stddev")
156 59 : return std::make_unique<VectorCalculator<InType, OutType, StdDev>>(other, item);
157 :
158 46 : else if (item == "stderr")
159 8 : return std::make_unique<VectorCalculator<InType, OutType, StdErr>>(other, item);
160 :
161 38 : else if (item == "norm2")
162 8 : return std::make_unique<VectorCalculator<InType, OutType, L2Norm>>(other, item);
163 :
164 30 : else if (item == "ratio")
165 8 : return std::make_unique<VectorCalculator<InType, OutType, Ratio>>(other, item);
166 :
167 22 : else if (item == "median")
168 8 : return std::make_unique<VectorCalculator<InType, OutType, Median>>(other, item);
169 :
170 14 : else if (item == "meanabs")
171 14 : 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 166 : 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 166 : if (item == "percentile")
189 150 : ptr = std::make_unique<Percentile<std::vector<InType>, std::vector<OutType>>>(
190 : other, item, levels, replicates, seed, calc);
191 16 : else if (item == "bca")
192 16 : 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 166 : 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
|
