---
title: "Flexible and consistent simulation of a matrix of Monte Carlo variates"
author: "Riccardo Porreca, Roland Schmid"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Matrix MC simulation}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---


```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE)
options(digits = 5) # for kable
linking_ok <- rTRNG::check_rTRNG_linking()
```


```{r code, include=FALSE, cache=FALSE}
source("utils/read_chunk_wrap.R", echo = FALSE, print.eval = FALSE)
read_chunk_wrap("code/mcMat.R")
read_chunk_wrap("code/mcMat.cpp")
read_chunk_wrap("code/mcMatParallel.cpp")
if (linking_ok) {
  Rcpp::sourceCpp("code/mcMat.cpp", verbose = FALSE, embeddedR = FALSE)
  Rcpp::sourceCpp("code/mcMatParallel.cpp", verbose = FALSE, embeddedR = FALSE)
}
```

Consider the Monte Carlo simulation of a matrix of i.i.d. normal random
variables. We will show how **rTRNG** can be used to perform a consistent
(*fair-playing*) simulation of a subset of the variables and simulations.

![](figures/mcMat.png){ width=85% }

## Consistent simulation in R

We rely on the TRNG engines exposed to R as reference classes by **rTRNG**.
```{r}
library(rTRNG)
```
The `mcMatR` function below performs the full sequential Monte Carlo simulation of `nrow`
normal i.i.d. samples of `ncol` variables using the `yarn2` generator.
```{r mcMatR}
```
A second function `mcSubMatR` relies on `jump` and `split` operations to perform
only a chunk [`startRow`, `endRow`] of simulations for a subset `subCols` of the
variables.
```{r mcSubMatR}
```
The parallel nature of the `yarn2` generator ensures the sub-simulation obtained
via `mcSubMatR` is consistent with the full sequential simulation.
```{r subMatExampleR}
```
```{r, echo=FALSE}
knitr::kable(cbind.data.frame(M = M, S = S), row.names = TRUE)
```


## Consistent simulation with Rcpp

We now use **Rcpp** to define functions `mcMatRcpp` and `mcSubMatRcpp` for the
full sequential simulation and the sub-simulation, respectively. The
`Rcpp::depends` attribute makes sure the TRNG library and headers shipped with
**rTRNG** are available to the C++ code. Moreover, `Rcpp::plugins(cpp11)`
enforces the C++11 standard required by TRNG >= 4.22.

```{Rcpp depends-h-ns, eval=FALSE}
```
```{Rcpp mcMatRcpp, eval=FALSE}
```
```{Rcpp mcSubMatRcpp, eval=FALSE}
```
As seen above for the R case, consistency of the simulation obtained via
`mcSubMatRcpp` with the full sequential simulation is guaranteed.
```{r subMatExampleRcpp, eval=linking_ok}
```
```{r, echo=FALSE, eval=linking_ok}
knitr::kable(cbind.data.frame(M = M, S = S), row.names = TRUE)
```


## Consistent parallel simulation with RcppParallel

The same technique used for generating a sub-set of the simulations can be
exploited for performing a parallel simulation in C++. We can embed the body of
`mcSubMatRcpp` above into an `RcppParallel::Worker` for performing chunks of
Monte Carlo simulations in parallel, for any subset `subCols` of the variables.
```{Rcpp mcMatRcppParallel, eval=FALSE}
```
The parallel nature of the `yarn2` generator ensures the parallel simulation is
playing fair, i.e. is consistent with the sequential simulation.
```{r fullMatExampleRcppParallel, eval=linking_ok}
```
Similarly, we can achieve a consistent parallel simulation of a subset of the
variables only.
```{r subMatExampleRcppParallel, eval=linking_ok}
```
```{r, echo=FALSE, eval=linking_ok}
knitr::kable(cbind.data.frame(M = M, Sp = Sp), row.names = TRUE)
```