Quinoa regression test code coverage report
Current view: top level - DiffEq/OrnsteinUhlenbeck - DiagOrnsteinUhlenbeck.hpp (source / functions) Hit Total Coverage
Commit: Quinoa_v0.3-298-g5666f42 Lines: 16 17 94.1 %
Date: 2021-01-09 20:32:40 Functions: 2 16 12.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 14 24 58.3 %

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/DiffEq/OrnsteinUhlenbeck/DiagOrnsteinUhlenbeck.hpp
       4                 :            :   \copyright 2012-2015 J. Bakosi,
       5                 :            :              2016-2018 Los Alamos National Security, LLC.,
       6                 :            :              2019-2020 Triad National Security, LLC.
       7                 :            :              All rights reserved. See the LICENSE file for details.
       8                 :            :   \brief     System of diagonal Ornstein-Uhlenbeck SDEs
       9                 :            :   \details   This file implements the time integration of a system of stochastic
      10                 :            :     differential equations (SDEs), with linear drift and constant diagonal
      11                 :            :     diffusion, whose invariant is the joint [normal
      12                 :            :     distribution](http://en.wikipedia.org/wiki/Normal_distribution).
      13                 :            : 
      14                 :            :     In a nutshell, the equation integrated governs a set of scalars,
      15                 :            :     \f$Y_\alpha\f$, \f$\alpha\!=\!1,\dots,N\f$, as
      16                 :            : 
      17                 :            :     @m_class{m-show-m}
      18                 :            : 
      19                 :            :     \f[
      20                 :            :        \mathrm{d}Y_\alpha(t) = \theta_\alpha\left(\mu_\alpha - Y_\alpha\right)
      21                 :            :        \mathrm{d}t + \sigma_\alpha\mathrm{d}W_\alpha(t), \qquad \alpha=1,\dots,N
      22                 :            :     \f]
      23                 :            : 
      24                 :            :     @m_class{m-hide-m}
      25                 :            : 
      26                 :            :     \f[ \begin{split}
      27                 :            :        \mathrm{d}Y_\alpha(t) = \theta_\alpha\left(\mu_\alpha - Y_\alpha\right)
      28                 :            :        \mathrm{d}t + \sigma_\alpha\mathrm{d}W_\alpha(t), \\ \alpha=1,\dots,N
      29                 :            :     \end{split} \f]
      30                 :            : 
      31                 :            :     with parameter vectors \f$\theta_\alpha > 0\f$, \f$\mu_\alpha\f$, and
      32                 :            :     \f$\sigma_\alpha > 0\f$. Here \f$\mathrm{d}W_\alpha(t)\f$ is an isotropic
      33                 :            :     vector-valued [Wiener process](http://en.wikipedia.org/wiki/Wiener_process)
      34                 :            :     with independent increments. The invariant distribution is the joint normal
      35                 :            :     distribution. This system of SDEs consists of N independent equations, each
      36                 :            :     being a single-variate [Ornstein-Uhlenbeck
      37                 :            :     process](http://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process).
      38                 :            : 
      39                 :            :     From the Fokker-Planck equation, equivalent to the SDE above, the equations
      40                 :            :     governing the means, \f$ \langle Y_\alpha \rangle\f$, are
      41                 :            :     \f[
      42                 :            :       \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
      43                 :            :       \frac{\partial\irmean{Y_\alpha}}{\partial t} =
      44                 :            :         \theta_\alpha\left(\mu_\alpha - \irmean{Y_\alpha}\right)
      45                 :            :     \f]
      46                 :            :     while the equation governing the covariance matrix, \f$ \langle y_\alpha
      47                 :            :     y_\beta \rangle \equiv \left\langle (Y_\alpha - \langle Y_\alpha \rangle)
      48                 :            :     (Y_\beta - \langle Y_\beta\rangle) \right\rangle \f$, is
      49                 :            : 
      50                 :            :     @m_class{m-show-m}
      51                 :            : 
      52                 :            :     \f[
      53                 :            :       \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
      54                 :            :       \newcommand{\irv}[1]{\langle{#1^2}\rangle}
      55                 :            :       \frac{\partial\irmean{y_\alpha y_\beta}}{\partial t} =
      56                 :            :         -\left(\theta_\alpha+\theta_\beta\right)\irmean{y_\alpha y_\beta} +
      57                 :            :         \left\{ \begin{array}{lr}
      58                 :            :           \sigma_\alpha^2 & \mathrm{for} \quad \alpha = \beta,\\
      59                 :            :           0               & \mathrm{for} \quad \alpha \ne \beta
      60                 :            :         \end{array} \right..
      61                 :            :     \f]
      62                 :            : 
      63                 :            :     @m_class{m-hide-m}
      64                 :            : 
      65                 :            :     \f[ \begin{split}
      66                 :            :       \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
      67                 :            :       \newcommand{\irv}[1]{\langle{#1^2}\rangle}
      68                 :            :       \frac{\partial\irmean{y_\alpha y_\beta}}{\partial t} =
      69                 :            :         -\left(\theta_\alpha+\theta_\beta\right)\irmean{y_\alpha y_\beta} \\ +
      70                 :            :         \left\{ \begin{array}{lr}
      71                 :            :           \sigma_\alpha^2 & \mathrm{for} \quad \alpha = \beta,\\
      72                 :            :           0               & \mathrm{for} \quad \alpha \ne \beta
      73                 :            :         \end{array} \right..
      74                 :            :     \end{split} \f]
      75                 :            : */
      76                 :            : // *****************************************************************************
      77                 :            : #ifndef DiagOrnsteinUhlenbeck_h
      78                 :            : #define DiagOrnsteinUhlenbeck_h
      79                 :            : 
      80                 :            : #include <vector>
      81                 :            : #include <cmath>
      82                 :            : 
      83                 :            : #include "InitPolicy.hpp"
      84                 :            : #include "DiagOrnsteinUhlenbeckCoeffPolicy.hpp"
      85                 :            : #include "RNG.hpp"
      86                 :            : #include "Particles.hpp"
      87                 :            : 
      88                 :            : namespace walker {
      89                 :            : 
      90                 :            : extern ctr::InputDeck g_inputdeck;
      91                 :            : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
      92                 :            : 
      93                 :            : //! \brief Diagonal Ornstein-Uhlenbeck SDE used polymorphically with DiffEq
      94                 :            : //! \details The template arguments specify policies and are used to configure
      95                 :            : //!   the behavior of the class. The policies are:
      96                 :            : //!   - Init - initialization policy, see DiffEq/InitPolicy.h
      97                 :            : //!   - Coefficients - coefficients policy, see
      98                 :            : //!       DiffEq/DiagOrnsteinUhlenbeckCoeffPolicy.h
      99                 :            : template< class Init, class Coefficients >
     100                 :            : class DiagOrnsteinUhlenbeck {
     101                 :            : 
     102                 :            :   private:
     103                 :            :     using ncomp_t = tk::ctr::ncomp_t;
     104                 :            : 
     105                 :            :   public:
     106                 :            :     //! \brief Constructor
     107                 :            :     //! \param[in] c Index specifying which system of diagonal
     108                 :            :     //!   Ornstein-Uhlenbeck SDEs to construct. There can be multiple diag_ou
     109                 :            :     //!   ... end blocks in a control file. This index specifies which diagonal
     110                 :            :     //!   Ornstein-Uhlenbeck SDE system to instantiate. The index corresponds to
     111                 :            :     //!   the order in which the diag_ou ... end blocks are given the control
     112                 :            :     //!   file.
     113                 :         57 :     explicit DiagOrnsteinUhlenbeck( ncomp_t c ) :
     114                 :            :       m_c( c ),
     115                 :            :       m_depvar( g_inputdeck.get< tag::param, tag::diagou, tag::depvar >().at(c) ),
     116                 :            :       m_ncomp( g_inputdeck.get< tag::component >().get< tag::diagou >().at(c) ),
     117                 :         57 :       m_offset( g_inputdeck.get< tag::component >().offset< tag::diagou >(c) ),
     118                 :         57 :       m_rng( g_rng.at( tk::ctr::raw(
     119                 :            :         g_inputdeck.get< tag::param, tag::diagou, tag::rng >().at(c) ) ) ),
     120                 :            :       m_sigmasq(),
     121                 :            :       m_theta(),
     122                 :            :       m_mu(),
     123                 :         57 :       coeff( m_ncomp,
     124                 :            :              g_inputdeck.get< tag::param, tag::diagou, tag::sigmasq >().at(c),
     125                 :            :              g_inputdeck.get< tag::param, tag::diagou, tag::theta >().at(c),
     126                 :            :              g_inputdeck.get< tag::param, tag::diagou, tag::mu >().at(c),
     127 [ -  + ][ -  + ]:        171 :              m_sigmasq, m_theta, m_mu ) {}
         [ -  + ][ -  + ]
                 [ +  - ]
     128                 :            : 
     129                 :            :     //! Initalize SDE, prepare for time integration
     130                 :            :     //! \param[in] stream Thread (or more precisely stream) ID 
     131                 :            :     //! \param[in,out] particles Array of particle properties 
     132                 :            :     void initialize( int stream, tk::Particles& particles ) {
     133                 :            :       //! Set initial conditions using initialization policy
     134                 :            :       Init::template
     135                 :            :         init< tag::diagou >
     136                 :          0 :             ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
     137                 :            :     }
     138                 :            : 
     139                 :            :     //! \brief Advance particles according to the system of diagonal
     140                 :            :     //!   Orsntein-Uhlenbeck SDEs
     141                 :            :     //! \param[in,out] particles Array of particle properties
     142                 :            :     //! \param[in] stream Thread (or more precisely stream) ID
     143                 :            :     //! \param[in] dt Time step size
     144                 :    1650000 :     void advance( tk::Particles& particles,
     145                 :            :                   int stream,
     146                 :            :                   tk::real dt,
     147                 :            :                   tk::real,
     148                 :            :                   const std::map< tk::ctr::Product, tk::real >& )
     149                 :            :     {
     150                 :    1650000 :       const auto npar = particles.nunk();
     151         [ +  + ]:  521650000 :       for (auto p=decltype(npar){0}; p<npar; ++p) {
     152                 :            :         // Generate Gaussian random numbers with zero mean and unit variance
     153                 :  520000000 :         std::vector< tk::real > dW( m_ncomp );
     154 [ +  - ][ +  - ]: 1040000000 :         m_rng.gaussian( stream, m_ncomp, dW.data() );
     155                 :            : 
     156                 :            :         // Advance all m_ncomp scalars
     157         [ +  + ]: 1560000000 :         for (ncomp_t i=0; i<m_ncomp; ++i) {
     158         [ +  - ]: 1040000000 :           tk::real& par = particles( p, i, m_offset );
     159         [ +  - ]: 1040000000 :           tk::real d = m_sigmasq[i] * dt;
     160         [ +  - ]: 1040000000 :           d = (d > 0.0 ? std::sqrt(d) : 0.0);
     161                 : 1040000000 :           par += m_theta[i]*(m_mu[i] - par)*dt + d*dW[i];
     162                 :            :         }
     163                 :            :       }
     164                 :    1650000 :     }
     165                 :            : 
     166                 :            :   private:
     167                 :            :     const ncomp_t m_c;                  //!< Equation system index
     168                 :            :     const char m_depvar;                //!< Dependent variable
     169                 :            :     const ncomp_t m_ncomp;              //!< Number of components
     170                 :            :     const ncomp_t m_offset;             //!< Offset SDE operates from
     171                 :            :     const tk::RNG& m_rng;               //!< Random number generator
     172                 :            : 
     173                 :            :     //! Coefficients
     174                 :            :     std::vector< kw::sde_sigmasq::info::expect::type > m_sigmasq;
     175                 :            :     std::vector< kw::sde_theta::info::expect::type > m_theta;
     176                 :            :     std::vector< kw::sde_mu::info::expect::type > m_mu;
     177                 :            : 
     178                 :            :     //! Coefficients policy
     179                 :            :     Coefficients coeff;
     180                 :            : };
     181                 :            : 
     182                 :            : } // walker::
     183                 :            : 
     184                 :            : #endif // DiagOrnsteinUhlenbeck_h

Generated by: LCOV version 1.12