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 "KernelDensity1D.h"
      11             : #include "math.h"
      12             : #include "libmesh/utility.h"
      13             : #include "DelimitedFileReader.h"
      14             : #include "SystemBase.h"
      15             : #include "Assembly.h"
      16             : #include "Normal.h"
      17             : #include "Uniform.h"
      18             : 
      19             : registerMooseObject("StochasticToolsApp", KernelDensity1D);
      20             : 
      21             : InputParameters
      22         160 : KernelDensity1D::validParams()
      23             : {
      24         160 :   InputParameters params = Distribution::validParams();
      25         320 :   MooseEnum bandwidth_rule("silverman standarddeviation userdefined");
      26         320 :   MooseEnum kernel_function("gaussian uniform");
      27         160 :   params.addClassDescription("KernelDensity1D distribution");
      28         320 :   params.addRequiredParam<MooseEnum>(
      29             :       "bandwidth_rule", bandwidth_rule, "Bandwidth rule for evaluating the bandwith.");
      30         320 :   params.addRequiredParam<MooseEnum>(
      31             :       "kernel_function",
      32             :       kernel_function,
      33             :       "Kernel function determines the shape of the underlying kernel for the kernel density.");
      34         480 :   params.addRangeCheckedParam<Real>("bandwidth",
      35         320 :                                     1.0,
      36             :                                     "bandwidth > 0",
      37             :                                     "Bandwidth controls the smoothness of the kernel density.");
      38         320 :   params.addParam<std::vector<Real>>("data", "The data vector.");
      39         320 :   params.addParam<FileName>("file_name", "Name of the CSV file.");
      40         320 :   params.addParam<std::string>(
      41             :       "file_column_name", "Name of column in csv file to use, by default first column is used.");
      42         160 :   return params;
      43         160 : }
      44             : 
      45          80 : KernelDensity1D::KernelDensity1D(const InputParameters & parameters)
      46             :   : Distribution(parameters),
      47          80 :     _bandwidth_rule(getParam<MooseEnum>("bandwidth_rule")),
      48         160 :     _kernel_function(getParam<MooseEnum>("kernel_function")),
      49         240 :     _bandwidth(getParam<Real>("bandwidth"))
      50             : {
      51         256 :   if (isParamValid("data") && isParamValid("file_name"))
      52           0 :     paramError("data", "data and file_name both cannot be set at the same time.");
      53         160 :   else if (isParamValid("file_name"))
      54             :   {
      55          64 :     MooseUtils::DelimitedFileReader reader(getParam<FileName>("file_name"));
      56          32 :     reader.read();
      57          64 :     if (isParamValid("file_column_name"))
      58           0 :       _data = reader.getData(getParam<std::string>("file_column_name"));
      59             :     else
      60          32 :       _data = reader.getData(0);
      61          32 :   }
      62          96 :   else if (isParamValid("data"))
      63         144 :     _data = getParam<std::vector<Real>>("data");
      64             :   else
      65           0 :     mooseError("Either 'data' or 'file_name' parameters must be specified to represent input data");
      66             :   Real mu = 0;
      67             :   Real sd = 0;
      68        3760 :   for (unsigned i = 0; i < _data.size(); ++i)
      69             :   {
      70        3680 :     mu += _data[i] / _data.size();
      71             :   }
      72        3760 :   for (unsigned i = 0; i < _data.size(); ++i)
      73             :   {
      74        3680 :     sd += Utility::pow<2>((_data[i] - mu)) / _data.size();
      75             :   }
      76          80 :   sd = std::pow(sd, 0.5);
      77          80 :   if (_bandwidth_rule == "silverman")
      78             :   {
      79          48 :     _bandwidth = 1.06 * sd * std::pow(_data.size(), -0.2);
      80             :   }
      81          32 :   else if (_bandwidth_rule == "standarddeviation")
      82             :   {
      83          16 :     _bandwidth = sd;
      84             :   }
      85          80 : }
      86             : 
      87             : Real
      88          80 : KernelDensity1D::pdf(const Real & x,
      89             :                      const Real & bandwidth,
      90             :                      const std::vector<Real> & data,
      91             :                      const MooseEnum & kernel_function)
      92             : {
      93             :   Real value = 0;
      94          80 :   if (kernel_function == "gaussian")
      95             :   {
      96        2144 :     for (unsigned i = 0; i < data.size(); ++i)
      97             :     {
      98        2080 :       value += 1 / (data.size() * bandwidth) * Normal::pdf(((x - data[i]) / bandwidth), 0.0, 1.0);
      99             :     }
     100             :   }
     101          16 :   else if (kernel_function == "uniform")
     102             :   {
     103        1616 :     for (unsigned i = 0; i < data.size(); ++i)
     104             :     {
     105        1600 :       value += 1 / (data.size() * bandwidth) * Uniform::pdf(((x - data[i]) / bandwidth), -1.0, 1.0);
     106             :     }
     107             :   }
     108             :   else
     109           0 :     ::mooseError("Invalid kernel function type ", std::string(kernel_function));
     110          80 :   return value;
     111             : }
     112             : 
     113             : Real
     114          80 : KernelDensity1D::cdf(const Real & x,
     115             :                      const Real & bandwidth,
     116             :                      const std::vector<Real> & data,
     117             :                      const MooseEnum & kernel_function)
     118             : {
     119             :   Real value = 0;
     120          80 :   if (kernel_function == "gaussian")
     121             :   {
     122        2144 :     for (unsigned i = 0; i < data.size(); ++i)
     123             :     {
     124        2080 :       value += 1.0 / (data.size()) * Normal::cdf(x, data[i], bandwidth);
     125             :     }
     126             :   }
     127          16 :   else if (kernel_function == "uniform")
     128             :   {
     129        1616 :     for (unsigned i = 0; i < data.size(); ++i)
     130             :     {
     131        1600 :       value += 1.0 / (data.size()) * Uniform::cdf((x - data[i]) / bandwidth, -1.0, 1.0);
     132             :     }
     133             :   }
     134             :   else
     135           0 :     ::mooseError("Invalid kernel function type ", std::string(kernel_function));
     136          80 :   return value;
     137             : }
     138             : 
     139             : Real
     140          80 : KernelDensity1D::quantile(const Real & p,
     141             :                           const Real & bandwidth,
     142             :                           const std::vector<Real> & data,
     143             :                           const MooseEnum & kernel_function)
     144             : {
     145             :   Real value = 0;
     146          80 :   if (kernel_function == "gaussian")
     147             :   {
     148          64 :     int index = std::round(p * (data.size() - 1));
     149          64 :     value = (Normal::quantile(p, data[index], bandwidth));
     150             :   }
     151          16 :   else if (kernel_function == "uniform")
     152             :   {
     153          16 :     int index = std::round(p * (data.size() - 1));
     154          16 :     value = (Uniform::quantile(p, -1.0, 1.0) * bandwidth + data[index]);
     155             :   }
     156             :   else
     157           0 :     ::mooseError("Invalid kernel function type ", std::string(kernel_function));
     158          80 :   return value;
     159             : }
     160             : 
     161             : Real
     162          80 : KernelDensity1D::pdf(const Real & x) const
     163             : {
     164          80 :   return pdf(x, _bandwidth, _data, _kernel_function);
     165             : }
     166             : 
     167             : Real
     168          80 : KernelDensity1D::cdf(const Real & x) const
     169             : {
     170          80 :   return cdf(x, _bandwidth, _data, _kernel_function);
     171             : }
     172             : 
     173             : Real
     174          80 : KernelDensity1D::quantile(const Real & p) const
     175             : {
     176          80 :   return quantile(p, _bandwidth, _data, _kernel_function);
     177             : }