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