Time Dependent Arma Markov

R/RcppExports.R

@@ -99,6 +99,10 @@ ArmaSimulateMarkov <- function(sim, cohort, transition, duration, state_cost, di
     .Call(`_SpeedyMarkov_ArmaSimulateMarkov`, sim, cohort, transition, duration, state_cost, discounting, qalys, intervention_cost)
 }
 
+ArmaTDMarkovLoop <- function(m_TR, a_P) {
+    .Call(`_SpeedyMarkov_ArmaTDMarkovLoop`, m_TR, a_P)
+}
+
 #' @title Arrange Vectorised Matrix Samples using Rcpp
 #' 
 #' @description A convenience function used to arrange vectorised matrix samples into the correct

src/ArmaTDMarkovLoop.cpp

@@ -0,0 +1,21 @@
+// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-
+
+// we only include RcppArmadillo.h which pulls Rcpp.h in for us
+#include "RcppArmadillo.h"
+
+// via the depends attribute we tell Rcpp to create hooks for
+// RcppArmadillo so that the build process will know what to do
+//
+// [[Rcpp::depends(RcppArmadillo)]]
+//
+// [[Rcpp::export]]
+arma::mat ArmaTDMarkovLoop(arma::mat m_TR, arma::cube& a_P )
+    {
+        int rows = m_TR.n_rows;
+
+        for(int i = 1; i < rows; i++){
+            m_TR.row(i) = m_TR.row(i-1) * a_P.slice(i-1);
+        }
+
+        return m_TR;
+    }

src/RcppExports.cpp

@@ -38,6 +38,18 @@ BEGIN_RCPP
     return rcpp_result_gen;
 END_RCPP
 }
+// ArmaTDMarkovLoop
+arma::mat ArmaTDMarkovLoop(arma::mat m_TR, arma::cube& a_P);
+RcppExport SEXP _SpeedyMarkov_ArmaTDMarkovLoop(SEXP m_TRSEXP, SEXP a_PSEXP) {
+BEGIN_RCPP
+    Rcpp::RObject rcpp_result_gen;
+    Rcpp::RNGScope rcpp_rngScope_gen;
+    Rcpp::traits::input_parameter< arma::mat >::type m_TR(m_TRSEXP);
+    Rcpp::traits::input_parameter< arma::cube& >::type a_P(a_PSEXP);
+    rcpp_result_gen = Rcpp::wrap(ArmaTDMarkovLoop(m_TR, a_P));
+    return rcpp_result_gen;
+END_RCPP
+}
 // MatrixArrange
 Rcpp::List MatrixArrange(Rcpp::List samples);
 RcppExport SEXP _SpeedyMarkov_MatrixArrange(SEXP samplesSEXP) {
@@ -53,6 +65,7 @@ END_RCPP
 static const R_CallMethodDef CallEntries[] = {
     {"_SpeedyMarkov_ArmaMarkovLoop", (DL_FUNC) &_SpeedyMarkov_ArmaMarkovLoop, 4},
     {"_SpeedyMarkov_ArmaSimulateMarkov", (DL_FUNC) &_SpeedyMarkov_ArmaSimulateMarkov, 8},
+    {"_SpeedyMarkov_ArmaTDMarkovLoop", (DL_FUNC) &_SpeedyMarkov_ArmaTDMarkovLoop, 2},
     {"_SpeedyMarkov_MatrixArrange", (DL_FUNC) &_SpeedyMarkov_MatrixArrange, 1},
     {NULL, NULL, 0}
 };

tests/testthat/testArmaTDMarkov.R

@@ -0,0 +1,88 @@
+# Test that the ArmaMarkov function matches
+# the Rloop function for a simple example.
+
+Rloopfunc <- function(m_TR, a_P) {
+  # throughout the number of cycles
+  for (t in 1:(nrow(m_TR) - 1)) {
+ # estimate the Markov trace for cycle the next cycle (t + 1)
+    m_TR[t + 1,] <-
+      m_TR[t,] %*% a_P[, , t]
+  }
+  # return the trace
+  m_TR
+}
+
+# create example array function
+makeTransProbArray <- function(simDim = 4){
+  # need the random number stream to be consistent- fix at seed = 100
+  set.seed(100)
+  # create an array from slicing this matrix
+  array(
+    data = #matrix(
+      runif(simDim  * simDim * 100),
+      #byrow = T,
+      #nrow = simDim,
+      #ncol = simDim
+   # ),
+    dim = c(simDim, simDim, 100)
+  )
+}
+
+# create example markov trace function
+makeMarkovTrace <- function(simDim = 4) {
+  # make empty trace
+  m_TR <- matrix(data = 0,
+                 nrow = 100,
+                 ncol = simDim)
+  # initialise trace
+  m_TR[1,]  <- c(1, rep(0, simDim - 1))
+  # return trace
+  m_TR
+
+}
+
+
+
+#==================================================#
+# EQUAL
+#==================================================#
+test_that("ArmaMarkov equals BaseRloop", {
+  # test three different symmetric matrix dimensions:
+  for(s in c(3,12,21)){
+    # in each case test that the two approaches get equal answers
+    testthat::expect_identical(
+    ArmaTDMarkovLoop(m_TR = makeMarkovTrace(s), a_P = makeTransProbArray(s)),
+    Rloopfunc(m_TR = makeMarkovTrace(s), a_P = makeTransProbArray(s))
+  ) # end expect equal
+  }
+}) # end testthat.
+
+
+
+#==================================================#
+# FASTER
+#==================================================#
+test_that("ArmaMarkov is faster than BaseRloop", {
+  # test three different symmetric matrix dimensions:
+  for(s in c(3,12,21)){
+    mb <- microbenchmark::microbenchmark(
+      Cpp = ArmaTDMarkovLoop(m_TR = makeMarkovTrace(s), a_P = makeTransProbArray(s)),
+      Rloop = Rloopfunc(m_TR = makeMarkovTrace(s), a_P = makeTransProbArray(s))
+    )
+
+
+    # in each case test that the ArmaMarkov function is faster
+    # Mean
+    testthat::expect_false(
+      mean(mb[mb$expr == "Cpp" ,"time"]) > mean(mb[mb$expr == "Rloop" ,"time"])
+      ) # end expect equal
+
+    # Median
+    testthat::expect_false(
+      median(mb[mb$expr == "Cpp" ,"time"]) > median(mb[mb$expr == "Rloop" ,"time"])
+    ) # end expect equal
+
+  }
+}) # end testthat.
+
+