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          32 : computeSTD(const std::vector<Real> & data, const unsigned int & start_index)
      22             : {
      23          32 :   if (data.size() < start_index)
      24             :     return 0.0;
      25             :   else
      26             :   {
      27          32 :     const Real mean = computeMean(data, start_index);
      28             :     const Real sq_diff =
      29          32 :         std::accumulate(data.begin() + start_index,
      30             :                         data.end(),
      31             :                         0.0,
      32         688 :                         [&mean](Real x, Real y) { return x + (y - mean) * (y - mean); });
      33          32 :     return std::sqrt(sq_diff / (data.size() - start_index));
      34             :   }
      35             : }
      36             : 
      37             : Real
      38          64 : computeMean(const std::vector<Real> & data, const unsigned int & start_index)
      39             : {
      40          64 :   if (data.size() < start_index)
      41             :     return 0.0;
      42             :   else
      43          64 :     return std::accumulate(data.begin() + start_index, data.end(), 0.0) /
      44          64 :            (data.size() - start_index);
      45             : }
      46             : 
      47             : std::vector<std::vector<Real>>
      48          30 : 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          30 :   Moose::indirectSort(outputs.begin(), outputs.end(), ind);
      55             : 
      56          30 :   std::vector<std::vector<Real>> out(inputs.size(), std::vector<Real>(samplessub * subset_prob));
      57          30 :   const size_t offset = std::ceil(samplessub * (1 - subset_prob));
      58          90 :   for (const auto & j : index_range(out))
      59         180 :     for (const auto & i : index_range(out[j]))
      60         120 :       out[j][i] = inputs[j][ind[i + offset]];
      61             : 
      62          30 :   return out;
      63          30 : }
      64             : 
      65             : std::vector<Real>
      66          14 : sortOutput(const std::vector<Real> & outputs, const unsigned int samplessub, const Real subset_prob)
      67             : {
      68             :   std::vector<size_t> ind;
      69          14 :   Moose::indirectSort(outputs.begin(), outputs.end(), ind);
      70             : 
      71          14 :   std::vector<Real> out(samplessub * subset_prob);
      72          14 :   const size_t offset = std::round(samplessub * (1 - subset_prob));
      73          42 :   for (const auto & i : index_range(out))
      74          28 :     out[i] = outputs[ind[i + offset]];
      75             : 
      76          14 :   return out;
      77          14 : }
      78             : 
      79             : Real
      80          14 : computeMin(const std::vector<Real> & data)
      81             : {
      82          14 :   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         300 : weightedResample(const std::vector<Real> & weights, Real rnd)
      96             : {
      97             :   unsigned int req_index = 0;
      98        1460 :   for (unsigned int i = 0; i < weights.size(); ++i)
      99             :   {
     100        1460 :     if (rnd < weights[i])
     101             :     {
     102             :       req_index = i;
     103             :       break;
     104             :     }
     105        1160 :     rnd -= weights[i];
     106             :   }
     107         300 :   return req_index;
     108             : }
     109             : } // namespace AdaptiveMonteCarloUtils