# Stochastic population model in R, Rcpp, and Fortran

The stochastic logistic population model described in this blog post has become my default exercise when I'm exploring a new programming language (Racket, F#, Rust). I ended that Rust post by noting that I was interested in Rust as a modern alternative for scientific computing and thought it would be a good learning exercise to re-write small legacy Fortran programs in Rust. In the process of looking for Fortran programs to translate to Rust, though, I found myself becoming more interested in the idea of learning Fortran than Rust, particularly after learning about the efforts to improve the tooling around Fortran (e.g., here and here). So, here we are...exploring Fortran via the stochastic population model exercise.

### Standalone Fortran Program

First, I wrote a standalone Fortran program of the stochastic population model. The main thing that tripped me up was reading arguments from the command line. I initialized a 1D array for holding the values of the arguments.

```
character(len=20), dimension(5) :: args
```

I initially missed the `len`

argument for `character`

and was confused by the results I was getting until I realized that the program was only reading the first character provided for each argument.

We can straightforwardly fill the `args`

array with a loop.

```
do i=1,command_argument_count()
call get_command_argument(i, args(i))
end do
```

But the syntax for converting the arguments from strings to the correct types is unusual.

```
read(args(1), *) t
```

`t`

is the first value in `args`

(indexed with `args(1)`

). The `*`

indicates the default format of the value, which the compiler infers is an `integer`

because that is the type of `t`

.

### Fortran Subroutines

I also created a file that only contains subroutines for calling from R. Two of the subroutines (`rstduniform`

and `rnormal`

) are mostly the same as from the main program, but the names are changed because the R help for the `.Fortran`

function recommended not using underscores in the subroutine names and the type of several parameters is `double precision`

, not `real`

, for compatibility with R. One of the subroutines, `rnorm`

, is only included as a test of `rnormal`

.

```
subroutine rnorm(reps, mu, sigma, arr)
implicit none
integer :: i
integer, intent(in) :: reps
double precision, intent(in) :: mu, sigma
double precision, intent(inout) :: arr(reps)
do i=1,reps
call rnormal(mu, sigma, arr(i))
end do
end subroutine rnorm
```

The `arr`

parameter is declared as `intent(inout)`

because on the R side we will pass in a vector that is filled on the Fortran side and then returned again to R. There are performance pitfalls with this approach if you are passing large vectors back-and-forth between Fortran and R, but those can be partly mitigated with the use of the dotCall64 package.

We compile this file into a shared object for use in R with

```
R CMD SHLIB stochastic-logistic-subroutines.f90
```

### Calling Fortran from R

In the R code, we load the shared object with

```
dyn.load("stochastic-logistic-subroutines.so")
```

To test that the `rnormal`

subroutine is giving reasonable results, we call the `rnorm`

subroutine as the first argument to `.Fortran`

. In this example, we draw 1 million random variates from a normal distribution with a mean of 7.8 and standard deviation of 8.3.

```
> n = 1000000L
> rnorm_result <- .Fortran("rnorm", n, 7.8, 8.3, vector("numeric", n))
# .Fortran returns a list; 4th element of list in this case contains the vector
> mean(rnorm_result[[4]])
[1] 7.798748
> sd(rnorm_result[[4]])
[1] 8.302705
```

Following the original blog post, we benchmark the R, Rcpp, and Fortran versions.

```
> library(microbenchmark)
> microbenchmark(
Rcpp = logmodc(t, yinit, r, k, thetasd),
Fortran = .Fortran("logmodf", t, yinit, r, k, thetasd, vector("numeric", t)),
R = logmodr(t, yinit, r, k, thetasd),
times = 500L
)
Unit: microseconds
expr min lq mean median uq max neval
Rcpp 13.235 14.6260 17.11903 15.5280 16.8100 86.612 500
Fortran 21.950 22.8965 26.76053 23.8760 25.3160 72.528 500
R 56.994 58.8775 67.25860 59.9085 62.3655 243.965 500
```

The original blog post (from 2017) observed a 19x speed up of the Rcpp version over the R version. Here we only get a 4x improvement presumably related to performance improvements in R. Our Fortran version is 2.5x faster than the R version.

Again, following the original blog post (with a slight modification), we benchmark running multiple simulations.

```
> rseq <- seq(1.1, 2.2, length.out = 10)
> microbenchmark(
Rcpp = lapply(rseq, function(x) logmodc(t, yinit, x, k, thetasd)),
Fortran = lapply(rseq, function(x) .Fortran("logmodf", t, yinit, x, k, thetasd, vector("numeric", t))),
R = lapply(rseq, function(x) logmodr(t, yinit, x, k, thetasd)),
times = 500L
)
Unit: microseconds
expr min lq mean median uq max neval
Rcpp 159.578 166.4295 214.4782 198.2240 202.9475 8265.555 500
Fortran 233.397 239.9620 264.1335 256.8515 262.1665 2994.857 500
R 621.029 635.1115 676.0114 665.7525 674.2735 2952.635 500
```

The relative performance gains for multiple simulations are similar to the single simulation.

### Conclusion

This post doesn't make a resounding case for using Fortran over C++ for speeding up simulations in R. The C++ code is shorter and runs faster. My C++ avoidance has centered on the complexity of the language. Perhaps there is a simpler subset of C++ that would go a long way towards speeding up my R code. However, this post makes the case that Fortran is still easier to learn than C++ for physics simulations.