Line data Source code
1 : // System includes
2 : #include <glob.h>
3 :
4 : // STL includes
5 : #include <cmath>
6 : #include <numeric>
7 :
8 : // MOOSE contrib
9 : #include "pcrecpp.h"
10 : #include "tinydir.h"
11 : #include "MooseUtils.h"
12 : #include "MooseRandom.h"
13 : #include "Lognormal.h"
14 : #include "Normal.h"
15 :
16 : // MASTODON includes
17 : #include "MastodonUtils.h"
18 :
19 : std::vector<std::vector<Real>>
20 354 : MastodonUtils::responseSpectrum(const Real & freq_start,
21 : const Real & freq_end,
22 : const unsigned int & freq_num,
23 : const std::vector<Real> & history_acc,
24 : const Real & xi,
25 : const Real & reg_dt)
26 : {
27 : std::vector<Real> freq_vec, per_vec, aspec_vec, vspec_vec, dspec_vec;
28 : Real logdf, om_n, om_d, dt2, dis1, vel1, acc1, dis2, vel2, acc2, pdmax, kd;
29 125208 : for (std::size_t n = 0; n < freq_num; ++n)
30 : {
31 : // Building the frequency vector and the period vector.
32 : // Frequencies are distributed uniformly in the log scale.
33 : // Periods are calculated as inverse of frequency
34 124854 : logdf = (std::log10(freq_end) - std::log10(freq_start)) / (freq_num - 1);
35 124854 : freq_vec.push_back(pow(10.0, std::log10(freq_start) + n * logdf));
36 124854 : per_vec.push_back(1.0 / freq_vec[n]);
37 124854 : om_n = 2.0 * 3.141593 * freq_vec[n]; // om_n = 2*pi*f
38 124854 : om_d = om_n * xi;
39 : dis1 = 0.0;
40 : vel1 = 0.0;
41 124854 : dt2 = reg_dt * reg_dt;
42 124854 : pdmax = 0.0;
43 124854 : acc1 = -1.0 * history_acc[0] - 2.0 * om_d * vel1 - om_n * om_n * dis1;
44 124854 : kd = 1.0 + om_d * reg_dt + dt2 * om_n * om_n / 4.0;
45 8010876 : for (std::size_t j = 0; j < history_acc.size(); ++j)
46 : {
47 7886022 : dis2 = ((1.0 + om_d * reg_dt) * dis1 + (reg_dt + 1.0 / 2.0 * om_d * dt2) * vel1 +
48 7886022 : dt2 / 4.0 * acc1 - dt2 / 4.0 * history_acc[j]) /
49 : kd;
50 7886022 : acc2 = 4.0 / dt2 * (dis2 - dis1) - 4.0 / reg_dt * vel1 - acc1;
51 7886022 : vel2 = vel1 + reg_dt / 2.0 * (acc1 + acc2);
52 7886022 : if (std::abs(dis2) > pdmax)
53 3448719 : pdmax = std::abs(dis2);
54 : dis1 = dis2;
55 : vel1 = vel2;
56 : acc1 = acc2;
57 : }
58 124854 : dspec_vec.push_back(pdmax);
59 124854 : vspec_vec.push_back(pdmax * om_n);
60 124854 : aspec_vec.push_back(pdmax * om_n * om_n);
61 : }
62 2478 : return {freq_vec, per_vec, dspec_vec, vspec_vec, aspec_vec};
63 708 : }
64 :
65 : std::vector<std::vector<Real>>
66 354 : MastodonUtils::regularize(const std::vector<Real> & history_acc,
67 : const std::vector<Real> & history_time,
68 : const Real & reg_dt)
69 : {
70 : std::vector<Real> reg_acc;
71 : std::vector<Real> reg_tme;
72 354 : Real cur_tme = history_time[0];
73 : Real cur_acc;
74 28134 : for (std::size_t i = 0; i < history_time.size() - 1; ++i)
75 : {
76 55782 : while (cur_tme >= history_time[i] && cur_tme <= history_time[i + 1])
77 : {
78 28002 : cur_acc = history_acc[i] + (cur_tme - history_time[i]) /
79 28002 : (history_time[i + 1] - history_time[i]) *
80 28002 : (history_acc[i + 1] - history_acc[i]);
81 28002 : reg_acc.push_back(cur_acc);
82 28002 : reg_tme.push_back(cur_tme);
83 28002 : cur_tme += reg_dt;
84 : }
85 : }
86 1416 : return {reg_tme, reg_acc};
87 708 : }
88 :
89 : bool
90 274 : MastodonUtils::checkEqualSize(const std::vector<std::vector<Real>> & vectors)
91 : {
92 553 : for (const auto & v : vectors)
93 : {
94 280 : if (v.size() != vectors[0].size())
95 : return false;
96 : }
97 : return true;
98 : }
99 :
100 : bool
101 17 : MastodonUtils::checkEqualSize(const std::vector<std::vector<Real> *> & vectors)
102 : {
103 68 : for (std::size_t i = 0; i < vectors.size(); i++)
104 : {
105 51 : if ((*vectors[i]).size() != (*vectors[0]).size())
106 : return false;
107 : }
108 : return true;
109 : }
110 :
111 : bool
112 249 : MastodonUtils::checkEqual(const std::vector<Real> & vector1,
113 : const std::vector<Real> & vector2,
114 : const Real percent_error)
115 : {
116 249 : if (vector1.size() != vector2.size())
117 : return false;
118 952 : for (std::size_t i = 0; i < vector1.size(); ++i)
119 : {
120 714 : if (!MooseUtils::absoluteFuzzyEqual(
121 714 : vector1[i], vector2[i], std::abs(vector1[i] * percent_error / 100)))
122 : return false;
123 : }
124 : return true;
125 : }
126 :
127 : bool
128 811779 : MastodonUtils::isNegativeOrZero(const std::vector<Real> & vector)
129 : {
130 4418932 : for (const auto & element : vector)
131 3607154 : if (element <= 0)
132 : return true;
133 : return false;
134 : }
135 :
136 : Real
137 1553 : MastodonUtils::mean(const std::vector<Real> & vector)
138 : {
139 : Real sum = std::accumulate(vector.begin(), vector.end(), 0.0);
140 1553 : return sum / vector.size();
141 : }
142 :
143 : std::vector<Real>
144 17 : MastodonUtils::mean(const std::vector<std::vector<Real> *> & history_acc)
145 : {
146 : // bool check = ;
147 17 : if (MastodonUtils::checkEqualSize(history_acc) < 1)
148 0 : mooseError("Input vectors are all not of equal size.");
149 : else
150 : {
151 : std::vector<Real> mean_acc;
152 17 : mean_acc.resize((*history_acc[0]).size());
153 : Real tmp_var = 0;
154 118 : for (std::size_t i = 0; i < (*history_acc[0]).size(); i++)
155 : {
156 404 : for (std::size_t j = 0; j < history_acc.size(); j++)
157 : {
158 303 : tmp_var = tmp_var + (*history_acc[j])[i];
159 : }
160 101 : mean_acc[i] = tmp_var / history_acc.size();
161 : tmp_var = 0;
162 : }
163 17 : return mean_acc;
164 17 : }
165 : }
166 :
167 : Real
168 46 : MastodonUtils::median(const std::vector<Real> & vector, const std::string & interpolation)
169 : {
170 46 : std::vector<Real> sorted_vector = vector;
171 46 : std::sort(sorted_vector.begin(), sorted_vector.end());
172 46 : if (sorted_vector.size() % 2 != 0)
173 25 : return sorted_vector[(sorted_vector.size() - 1) / 2];
174 21 : else if (interpolation == "linear")
175 19 : return (sorted_vector[sorted_vector.size() / 2] + sorted_vector[sorted_vector.size() / 2 - 1]) /
176 19 : 2.0;
177 2 : else if (interpolation == "lower")
178 1 : return sorted_vector[sorted_vector.size() / 2 - 1];
179 1 : else if (interpolation == "higher")
180 1 : return sorted_vector[sorted_vector.size() / 2];
181 : else
182 0 : mooseError("Invalid interpolation type in median calculation.");
183 46 : }
184 :
185 : Real
186 6 : MastodonUtils::percentile(const std::vector<Real> & vector,
187 : const Real & percent,
188 : const std::string & interpolation)
189 : {
190 6 : std::vector<Real> sorted_vector = vector;
191 6 : std::sort(sorted_vector.begin(), sorted_vector.end());
192 6 : if (percent < 0.0 || percent > 100.0)
193 0 : mooseError("Percent should be between 0 and 100.\n");
194 : std::size_t low_index;
195 6 : if (floor(percent / 100 * sorted_vector.size()) == 0)
196 : low_index = 0;
197 : else
198 5 : low_index = floor(percent / 100 * sorted_vector.size()) - 1;
199 6 : if (interpolation == "lower")
200 1 : return sorted_vector[low_index];
201 5 : else if (interpolation == "higher")
202 1 : return sorted_vector[low_index + 1];
203 4 : else if (interpolation == "linear")
204 : {
205 4 : Real index_remainder = fmod(percent / 100.0 * sorted_vector.size(), 1.0);
206 : Real percentile_value =
207 4 : sorted_vector[low_index] +
208 4 : index_remainder * (sorted_vector[low_index + 1] - sorted_vector[low_index]);
209 4 : return percentile_value;
210 : }
211 : else
212 0 : mooseError("Invalid interpolation type in percentile calculation.");
213 6 : }
214 :
215 : Real
216 44 : MastodonUtils::standardDeviation(const std::vector<Real> & vector)
217 : {
218 : Real sum = 0;
219 820 : for (std::size_t i = 0; i < vector.size(); i++)
220 : {
221 776 : sum += (vector[i] - MastodonUtils::mean(vector)) * (vector[i] - MastodonUtils::mean(vector));
222 : }
223 44 : Real sigma = sqrt(sum / (vector.size() - 1));
224 44 : return sigma;
225 : }
226 :
227 : Real
228 43 : MastodonUtils::lognormalStandardDeviation(const std::vector<Real> & vector)
229 : {
230 43 : if (MastodonUtils::isNegativeOrZero(vector))
231 0 : mooseError("One or more elements in the sample for calculating beta are non positive.\n");
232 : std::vector<Real> logvector;
233 809 : for (const auto & element : vector)
234 766 : logvector.push_back(log(element));
235 86 : return MastodonUtils::standardDeviation(logvector);
236 43 : }
237 :
238 : std::string
239 218 : MastodonUtils::zeropad(const unsigned int n, const unsigned int n_tot)
240 : {
241 218 : std::size_t zeropadsize = (std::to_string(n_tot)).length() - (std::to_string(n)).length();
242 218 : std::string pad = "";
243 403 : for (std::size_t i = 0; i < zeropadsize; i++)
244 : {
245 : pad += "0";
246 : }
247 436 : return pad + std::to_string(n);
248 : }
249 :
250 : std::vector<std::string>
251 22 : MastodonUtils::glob(const std::string & pattern)
252 : {
253 : glob_t glob_result;
254 22 : glob(pattern.c_str(), GLOB_TILDE, NULL, &glob_result);
255 : std::vector<string> output;
256 198 : for (unsigned int i = 0; i < glob_result.gl_pathc; ++i)
257 352 : output.push_back(string(glob_result.gl_pathv[i]));
258 22 : globfree(&glob_result);
259 22 : return output;
260 0 : }
261 :
262 : std::vector<Real>
263 1027 : MastodonUtils::adjust(const std::vector<Real> & input, const Real & scale, const Real & offset)
264 : {
265 1027 : std::vector<Real> output(input.size());
266 11156 : for (std::size_t i = 0; i < input.size(); ++i)
267 10129 : output[i] = scale * input[i] + offset;
268 1027 : return output;
269 : }
270 :
271 : std::vector<Real>
272 16 : MastodonUtils::log10(const std::vector<Real> & input)
273 : {
274 16 : std::vector<Real> output(input.size());
275 88 : for (std::size_t i = 0; i < input.size(); ++i)
276 : {
277 74 : if (input[i] <= 0)
278 2 : ::mooseError("Cannot take the log of ", input[i], ".");
279 72 : output[i] = std::log10(input[i]);
280 : }
281 14 : return output;
282 2 : }
283 :
284 : Real
285 43 : MastodonUtils::greaterProbability(Real demand_median,
286 : Real demand_scale,
287 : Real capacity_median,
288 : Real capacity_scale)
289 : {
290 : // Need P(D-C > 0) = P(lnD-lnC > 0)
291 : // lnD and lnC are Normally distributed
292 : // Evaluating the mean and std dev of lnD-lnC as
293 : // mean = mean(lnD) - mean(lnC) = ln(thetaD) - ln(thetaC)
294 : // std dev = srss(betaD, betaC)
295 : // now, greater probability = 1 - CDF(0.0)
296 43 : Real greater_prob_location = log(demand_median) - log(capacity_median);
297 : Real greater_prob_scale =
298 43 : std::sqrt(demand_scale * demand_scale + capacity_scale * capacity_scale);
299 43 : return (1.0 - Normal::cdf(0.0, greater_prob_location, greater_prob_scale));
300 : }
301 :
302 : Real
303 405597 : MastodonUtils::calcLogLikelihood(const std::vector<Real> & im,
304 : const std::vector<Real> & pf,
305 : const Real & loc,
306 : const Real & sca,
307 : const unsigned int & n)
308 : {
309 : // error check
310 405597 : if (im.size() != pf.size())
311 0 : mooseError("While calculating loglikelihood, intensity measure and failure probability vectors "
312 : "should be of the same size.");
313 405597 : if (MastodonUtils::isNegativeOrZero(im) || MastodonUtils::isNegativeOrZero(pf))
314 0 : mooseError("While calculating loglikelihood, intensity measure or failure probability has a "
315 : "non-positive value.");
316 2207988 : for (std::size_t i = 0; i < pf.size(); ++i)
317 : {
318 1802391 : if (pf[i] > 1.0)
319 0 : mooseError("While calculating loglikelihood, a value greater than 1 is found in the failure "
320 : "probability vector.");
321 : }
322 405597 : if (sca <= 0)
323 0 : mooseError("While calculating loglikelihood, scale parameter should be positive.");
324 405597 : if (loc <= 0)
325 0 : mooseError("While calculating loglikelihood, location parameter should be positive.");
326 :
327 : // Calculating the likelihood at each IM value
328 : Real loglikelihood = 0;
329 2207988 : for (std::size_t i = 0; i < im.size(); ++i)
330 : {
331 : // Binomial pdf in the below calculation is made in the log scale due to
332 : // numerical errors (+inf) in calculating nCr for large trial sizes (n).
333 : // Calculating log10(nCr)
334 1802391 : unsigned int r = floor(n * pf[i]);
335 : Real log10_nCr = 0.0;
336 259605552 : for (std::size_t k = 1; k <= r; k++)
337 257803161 : log10_nCr += std::log10(n - r + k) - std::log10(k);
338 :
339 1802391 : Real p = Lognormal::cdf(im[i], log(loc), sca);
340 1802391 : if (p == 0)
341 : p = std::numeric_limits<Real>::min();
342 :
343 : // calculating sum of loglikelihoods. Each loglikelihood is the log(binomial pdf),
344 : // which ends up to be the summation below.
345 1802391 : loglikelihood += log10_nCr + r * std::log10(p) + (n - r) * std::log10(1.0 - p);
346 : }
347 405597 : return loglikelihood;
348 : }
349 :
350 : std::vector<Real>
351 11 : MastodonUtils::maximizeLogLikelihood(const std::vector<Real> & im,
352 : const std::vector<Real> & pf,
353 : const std::vector<Real> & loc_space,
354 : const std::vector<Real> & sca_space,
355 : const unsigned int & n,
356 : const bool & brute_force,
357 : const Real tolerance,
358 : const Real gamma,
359 : const int num_rnd,
360 : const int seed)
361 : {
362 11 : std::vector<Real> params_return = {0, 0};
363 11 : if (brute_force)
364 : {
365 4 : Real loglikelihoodmax = MastodonUtils::calcLogLikelihood(im, pf, loc_space[0], sca_space[0], n);
366 1547 : for (Real loc = loc_space[0]; loc < loc_space[1]; loc += 0.01)
367 304846 : for (Real sca = sca_space[0]; sca < sca_space[1]; sca += 0.01)
368 303303 : if (MastodonUtils::calcLogLikelihood(im, pf, loc, sca, n) >= loglikelihoodmax)
369 : {
370 776 : loglikelihoodmax = MastodonUtils::calcLogLikelihood(im, pf, loc, sca, n);
371 776 : params_return = {loc, sca};
372 : }
373 : }
374 : else
375 : // Using Randomized Gradient Descent to maximize likelihood (as defined above) or minimize
376 : // -loglikelihood
377 : {
378 : Real loc_rand;
379 : Real sca_rand;
380 7 : std::vector<Real> params_now = {0, 0}; // Initializing new parameter vector here.
381 7 : std::vector<Real> params_before = {0, 0}; // Initializing old parameter vector here.
382 7 : std::vector<Real> gradient_now = {-1, -2}; // Initializing the gradient vector here.
383 : // This variable will get updated within each iteration of the Gradient Descent algorithm.
384 : Real dparam =
385 : 0.01; // Initilizing an arbitrarily small deviation to the random seed parameter vector.
386 : // Note that Gradient Descent algorithm requires two likelihood values from two seeds.
387 : Real likelihood_now; // Initializing a variable.
388 : Real likelihood_before; // Initializing a variable.
389 : Real likelihood_base = std::numeric_limits<Real>::max();
390 : // This variable will get updated if a parameter vector has greater likelihood.
391 7 : MooseRandom::seed(seed); // Setting up the random number generator.
392 18341 : for (int index = 0; index < num_rnd; index++)
393 : {
394 18334 : loc_rand = loc_space[0] + (loc_space[1] - loc_space[0]) * MooseRandom::rand();
395 18334 : sca_rand = sca_space[0] + (sca_space[1] - sca_space[0]) * MooseRandom::rand();
396 18334 : likelihood_now = -MastodonUtils::calcLogLikelihood(im, pf, loc_rand, sca_rand, n);
397 18334 : likelihood_before =
398 18334 : -MastodonUtils::calcLogLikelihood(im, pf, loc_rand + dparam, sca_rand + dparam, n);
399 18334 : params_now = {loc_rand, sca_rand};
400 18334 : params_before = {loc_rand + dparam, sca_rand + dparam};
401 18334 : if (likelihood_now > likelihood_before) // Swap params and likelihoods before and now
402 : {
403 10406 : likelihood_now =
404 10406 : -MastodonUtils::calcLogLikelihood(im, pf, loc_rand + dparam, sca_rand + dparam, n);
405 10406 : likelihood_before = -MastodonUtils::calcLogLikelihood(im, pf, loc_rand, sca_rand, n);
406 10406 : params_now = {loc_rand + dparam, sca_rand + dparam};
407 10406 : params_before = {loc_rand, sca_rand};
408 : }
409 62358 : while (std::abs(likelihood_now - likelihood_before) > tolerance)
410 : {
411 56258 : gradient_now[0] = (likelihood_now - likelihood_before) / (params_now[0] - params_before[0]);
412 56258 : gradient_now[1] = (likelihood_now - likelihood_before) / (params_now[1] - params_before[1]);
413 : likelihood_before = likelihood_now;
414 56258 : params_now[0] = (params_now[0] - gamma * gradient_now[0]);
415 56258 : params_now[1] = (params_now[1] - gamma * gradient_now[1]);
416 52110 : if ((params_now[0] < loc_space[0]) || (params_now[0] > loc_space[1]) ||
417 101739 : (params_now[1] < sca_space[0]) || (params_now[1] > sca_space[1]))
418 : {
419 12234 : index = index - 1;
420 12234 : break;
421 : }
422 : else
423 : {
424 44024 : likelihood_now =
425 44024 : -MastodonUtils::calcLogLikelihood(im, pf, params_now[0], params_now[1], n);
426 44024 : if (likelihood_now < likelihood_base)
427 : {
428 : likelihood_base = likelihood_now;
429 128 : params_return = params_now;
430 : }
431 : }
432 : }
433 : }
434 7 : }
435 11 : return params_return;
436 0 : }
|
