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

           Branch data     Line data    Source code
       1                 :            : // *****************************************************************************
       2                 :            : /*!
       3                 :            :   \file      src/DiffEq/OrnsteinUhlenbeck/OrnsteinUhlenbeck.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 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 diffusion,
      11                 :            :     whose invariant is the [joint normal
      12                 :            :     distribution](http://en.wikipedia.org/wiki/Multivariate_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 + \sum_{\beta=1}^N \sigma_{\alpha\beta}\mathrm{d}W_\beta(t),
      22                 :            :        \qquad \alpha=1,\dots,N
      23                 :            :     \f]
      24                 :            : 
      25                 :            :     @m_class{m-hide-m}
      26                 :            : 
      27                 :            :     \f[ \begin{split}
      28                 :            :        \mathrm{d}Y_\alpha(t) = \theta_\alpha\left(\mu_\alpha - Y_\alpha\right)
      29                 :            :        \mathrm{d}t + \sum_{\beta=1}^N \sigma_{\alpha\beta}\mathrm{d}W_\beta(t),
      30                 :            :        \\ \alpha=1,\dots,N
      31                 :            :     \end{split} \f]
      32                 :            : 
      33                 :            :     with parameter vectors \f$\theta_\alpha > 0\f$, \f$\mu_\alpha\f$, and
      34                 :            :     symmetric positive semi-definite diffusion matrix
      35                 :            :     \f$\sigma_{\alpha\beta}\f$. Here \f$\mathrm{d}W_\beta(t)\f$ is an isotropic
      36                 :            :     vector-valued [Wiener process](http://en.wikipedia.org/wiki/Wiener_process)
      37                 :            :     with independent increments. The invariant distribution is the joint normal
      38                 :            :     distribution. This system of SDEs consists of N coupled equations, each
      39                 :            :     being a single-variate [Ornstein-Uhlenbeck
      40                 :            :     process](http://en.wikipedia.org/wiki/Ornstein%E2%80%93Uhlenbeck_process).
      41                 :            : 
      42                 :            :     From the Fokker-Planck equation, equivalent to the SDE above, the equations
      43                 :            :     governing the means, \f$ \langle Y_\alpha \rangle\f$, are
      44                 :            :     \f[
      45                 :            :       \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
      46                 :            :       \frac{\partial\irmean{Y_\alpha}}{\partial t} =
      47                 :            :         \theta_\alpha\left(\mu_\alpha - \irmean{Y_\alpha}\right)
      48                 :            :     \f]
      49                 :            :     while the equation governing the covariance matrix, \f$ \langle y_\alpha
      50                 :            :     y_\beta \rangle \equiv \left\langle (Y_\alpha - \langle Y_\alpha \rangle)
      51                 :            :     (Y_\beta - \langle Y_\beta\rangle) \right\rangle \f$, is
      52                 :            : 
      53                 :            :     @m_class{m-show-m}
      54                 :            : 
      55                 :            :     \f[
      56                 :            :       \newcommand{\irmean}[1]{{\langle{#1}\rangle}}
      57                 :            :       \newcommand{\irv}[1]{\langle{#1^2}\rangle}
      58                 :            :       \frac{\partial\irmean{y_\alpha y_\beta}}{\partial t} =
      59                 :            :       -\left(\theta_\alpha+\theta_\beta\right)\irmean{y_\alpha y_\beta}
      60                 :            :       +\sum_{\gamma=1}^N \sigma_{\alpha\gamma} \sigma_{\gamma\beta}.
      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                 :            :       +\sum_{\gamma=1}^N \sigma_{\alpha\gamma} \sigma_{\gamma\beta}.
      71                 :            :     \end{split} \f]
      72                 :            : */
      73                 :            : // *****************************************************************************
      74                 :            : #ifndef OrnsteinUhlenbeck_h
      75                 :            : #define OrnsteinUhlenbeck_h
      76                 :            : 
      77                 :            : #include <vector>
      78                 :            : #include <cmath>
      79                 :            : 
      80                 :            : #include "QuinoaConfig.hpp"
      81                 :            : 
      82                 :            : #ifdef HAS_MKL
      83                 :            :   #include <mkl_lapacke.h>
      84                 :            : #else
      85                 :            :   #include <lapacke.h>
      86                 :            : #endif
      87                 :            : 
      88                 :            : #include "InitPolicy.hpp"
      89                 :            : #include "OrnsteinUhlenbeckCoeffPolicy.hpp"
      90                 :            : #include "RNG.hpp"
      91                 :            : #include "Particles.hpp"
      92                 :            : 
      93                 :            : namespace walker {
      94                 :            : 
      95                 :            : extern ctr::InputDeck g_inputdeck;
      96                 :            : extern std::map< tk::ctr::RawRNGType, tk::RNG > g_rng;
      97                 :            : 
      98                 :            : //! \brief Ornstein-Uhlenbeck SDE used polymorphically with DiffEq
      99                 :            : //! \details The template arguments specify policies and are used to configure
     100                 :            : //!   the behavior of the class. The policies are:
     101                 :            : //!   - Init - initialization policy, see DiffEq/InitPolicy.h
     102                 :            : //!   - Coefficients - coefficients policy, see
     103                 :            : //!       DiffEq/OrnsteinUhlenbeckCoeffPolicy.h
     104                 :            : template< class Init, class Coefficients >
     105                 :            : class OrnsteinUhlenbeck {
     106                 :            : 
     107                 :            :   private:
     108                 :            :     using ncomp_t = tk::ctr::ncomp_t;
     109                 :            : 
     110                 :            :   public:
     111                 :            :     //! \brief Constructor
     112                 :            :     //! \param[in] c Index specifying which system of Ornstein-Uhlenbeck SDEs to
     113                 :            :     //!   construct. There can be multiple ornstein-uhlenbeck ... end blocks in
     114                 :            :     //!   a control file. This index specifies which Ornstein-Uhlenbeck SDE
     115                 :            :     //!   system to instantiate. The index corresponds to the order in which the
     116                 :            :     //!   ornstein-uhlenbeck ... end blocks are given the control file.
     117                 :         46 :     explicit OrnsteinUhlenbeck( ncomp_t c ) :
     118                 :            :       m_c( c ),
     119                 :            :       m_depvar( g_inputdeck.get< tag::param, tag::ou, tag::depvar >().at(c) ),
     120                 :            :       m_ncomp( g_inputdeck.get< tag::component >().get< tag::ou >().at(c) ),
     121                 :         46 :       m_offset( g_inputdeck.get< tag::component >().offset< tag::ou >(c) ),
     122                 :         46 :       m_rng( g_rng.at( tk::ctr::raw(
     123                 :            :         g_inputdeck.get< tag::param, tag::ou, tag::rng >().at(c) ) ) ),
     124                 :            :       m_sigma(),
     125                 :            :       m_theta(),
     126                 :            :       m_mu(),
     127                 :         46 :       coeff( m_ncomp,
     128                 :            :              g_inputdeck.get< tag::param, tag::ou, tag::sigmasq >().at(c),
     129                 :            :              g_inputdeck.get< tag::param, tag::ou, tag::theta >().at(c),
     130                 :            :              g_inputdeck.get< tag::param, tag::ou, tag::mu >().at(c),
     131 [ -  + ][ -  + ]:        138 :              m_sigma, m_theta, m_mu )
         [ -  + ][ -  + ]
                 [ +  - ]
     132                 :            :     {
     133                 :            :       // Compute diffusion matrix using Cholesky-decomposition
     134                 :         46 :       lapack_int n = static_cast< lapack_int >( m_ncomp );
     135                 :            :       #ifndef NDEBUG
     136                 :            :       lapack_int info =
     137                 :            :       #endif
     138         [ +  - ]:         46 :         LAPACKE_dpotrf( LAPACK_ROW_MAJOR, 'U', n, m_sigma.data(), n );
     139                 :            :       Assert( info == 0, "Error in Cholesky-decomposition" );
     140                 :         46 :     }
     141                 :            : 
     142                 :            :     //! Initalize SDE, prepare for time integration
     143                 :            :     //! \param[in] stream Thread (or more precisely stream) ID 
     144                 :            :     //! \param[in,out] particles Array of particle properties 
     145                 :            :     void initialize( int stream, tk::Particles& particles ) {
     146                 :            :       //! Set initial conditions using initialization policy
     147                 :            :       Init::template
     148                 :            :         init< tag::ou >
     149                 :          0 :             ( g_inputdeck, m_rng, stream, particles, m_c, m_ncomp, m_offset );
     150                 :            :     }
     151                 :            : 
     152                 :            :     //! \brief Advance particles according to the system of Orsntein-Uhlenbeck
     153                 :            :     //!   SDEs
     154                 :            :     //! \param[in,out] particles Array of particle properties
     155                 :            :     //! \param[in] stream Thread (or more precisely stream) ID
     156                 :            :     //! \param[in] dt Time step size
     157                 :      59000 :     void advance( tk::Particles& particles,
     158                 :            :                   int stream,
     159                 :            :                   tk::real dt,
     160                 :            :                   tk::real,
     161                 :            :                   const std::map< tk::ctr::Product, tk::real >& )
     162                 :            :     {
     163                 :      59000 :       const auto npar = particles.nunk();
     164         [ +  + ]:  105059000 :       for (auto p=decltype(npar){0}; p<npar; ++p) {
     165                 :            :         // Generate Gaussian random numbers with zero mean and unit variance
     166                 :  105000000 :         std::vector< tk::real > dW( m_ncomp );
     167         [ +  - ]:  105000000 :         m_rng.gaussian( stream, m_ncomp, dW.data() );
     168                 :            : 
     169                 :            :         // Advance all m_ncomp scalars
     170         [ +  + ]:  420000000 :         for (ncomp_t i=0; i<m_ncomp; ++i) {
     171                 :  315000000 :           tk::real& par = particles( p, i, m_offset );
     172                 :  315000000 :           par += m_theta[i]*(m_mu[i] - par)*dt;
     173         [ +  + ]: 1260000000 :           for (ncomp_t j=0; j<m_ncomp; ++j) {
     174                 :  945000000 :             tk::real d = m_sigma[ j*m_ncomp+i ] * sqrt(dt);     // use transpose
     175                 :  945000000 :             par += d*dW[j];
     176                 :            :           }
     177                 :            :         }
     178                 :            :       }
     179                 :      59000 :     }
     180                 :            : 
     181                 :            :   private:
     182                 :            :     const ncomp_t m_c;                  //!< Equation system index
     183                 :            :     const char m_depvar;                //!< Dependent variable
     184                 :            :     const ncomp_t m_ncomp;              //!< Number of components
     185                 :            :     const ncomp_t m_offset;             //!< Offset SDE operates from
     186                 :            :     const tk::RNG& m_rng;               //!< Random number generator
     187                 :            : 
     188                 :            :     //! Coefficients
     189                 :            :     std::vector< kw::sde_sigmasq::info::expect::type > m_sigma;
     190                 :            :     std::vector< kw::sde_theta::info::expect::type > m_theta;
     191                 :            :     std::vector< kw::sde_mu::info::expect::type > m_mu;
     192                 :            : 
     193                 :            :     //! Coefficients policy
     194                 :            :     Coefficients coeff;
     195                 :            : };
     196                 :            : 
     197                 :            : } // walker::
     198                 :            : 
     199                 :            : #endif // OrnsteinUhlenbeck_h

Generated by: LCOV version 1.12