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 "AdaptiveMonteCarloUtils.h"
      11             : #include "IndirectSort.h"
      12             : #include "libmesh/int_range.h"
      13             : 
      14             : /* AdaptiveMonteCarloUtils contains functions that are used across the Adaptive Monte
      15             :  Carlo set of algorithms.*/
      16             : 
      17             : namespace AdaptiveMonteCarloUtils
      18             : {
      19             : 
      20             : Real
      21          44 : computeSTD(const std::vector<Real> & data, const unsigned int & start_index)
      22             : {
      23          44 :   if (data.size() < start_index)
      24             :     return 0.0;
      25             :   else
      26             :   {
      27          44 :     const Real mean = computeMean(data, start_index);
      28             :     const Real sq_diff =
      29          44 :         std::accumulate(data.begin() + start_index,
      30             :                         data.end(),
      31             :                         0.0,
      32         946 :                         [&mean](Real x, Real y) { return x + (y - mean) * (y - mean); });
      33          44 :     return std::sqrt(sq_diff / (data.size() - start_index));
      34             :   }
      35             : }
      36             : 
      37             : Real
      38          88 : computeMean(const std::vector<Real> & data, const unsigned int & start_index)
      39             : {
      40          88 :   if (data.size() < start_index)
      41             :     return 0.0;
      42             :   else
      43          88 :     return std::accumulate(data.begin() + start_index, data.end(), 0.0) /
      44          88 :            (data.size() - start_index);
      45             : }
      46             : 
      47             : std::vector<std::vector<Real>>
      48          68 : sortInput(const std::vector<std::vector<Real>> & inputs,
      49             :           const std::vector<Real> & outputs,
      50             :           const unsigned int samplessub,
      51             :           const Real subset_prob)
      52             : {
      53             :   std::vector<size_t> ind;
      54          68 :   Moose::indirectSort(outputs.begin(), outputs.end(), ind);
      55             : 
      56          68 :   std::vector<std::vector<Real>> out(inputs.size(), std::vector<Real>(samplessub * subset_prob));
      57          68 :   const size_t offset = std::ceil(samplessub * (1 - subset_prob));
      58         204 :   for (const auto & j : index_range(out))
      59         408 :     for (const auto & i : index_range(out[j]))
      60         272 :       out[j][i] = inputs[j][ind[i + offset]];
      61             : 
      62          68 :   return out;
      63          68 : }
      64             : 
      65             : std::vector<Real>
      66          34 : sortOutput(const std::vector<Real> & outputs, const unsigned int samplessub, const Real subset_prob)
      67             : {
      68             :   std::vector<size_t> ind;
      69          34 :   Moose::indirectSort(outputs.begin(), outputs.end(), ind);
      70             : 
      71          34 :   std::vector<Real> out(samplessub * subset_prob);
      72          34 :   const size_t offset = std::round(samplessub * (1 - subset_prob));
      73         102 :   for (const auto & i : index_range(out))
      74          68 :     out[i] = outputs[ind[i + offset]];
      75             : 
      76          34 :   return out;
      77          34 : }
      78             : 
      79             : Real
      80          34 : computeMin(const std::vector<Real> & data)
      81             : {
      82          34 :   return *std::min_element(data.begin(), data.end());
      83             : }
      84             : 
      85             : std::vector<Real>
      86           0 : computeVectorABS(const std::vector<Real> & data)
      87             : {
      88           0 :   std::vector<Real> data_abs(data.size());
      89           0 :   for (unsigned int i = 0; i < data.size(); ++i)
      90           0 :     data_abs[i] = std::abs(data[i]);
      91           0 :   return data_abs;
      92             : }
      93             : 
      94             : unsigned int
      95         660 : weightedResample(const std::vector<Real> & weights, Real rnd)
      96             : {
      97             :   unsigned int req_index = 0;
      98        3212 :   for (unsigned int i = 0; i < weights.size(); ++i)
      99             :   {
     100        3212 :     if (rnd < weights[i])
     101             :     {
     102             :       req_index = i;
     103             :       break;
     104             :     }
     105        2552 :     rnd -= weights[i];
     106             :   }
     107         660 :   return req_index;
     108             : }
     109             : } // namespace AdaptiveMonteCarloUtils