From 23e50b1981ac91a448871741ff7fb39aa6032ac1 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sun, 22 Nov 2020 16:42:41 +0000 Subject: [PATCH 1/3] Added code for time dependant ArmaMarkovLoop. --- src/ArmaTDMarkovLoop.cpp | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) create mode 100644 src/ArmaTDMarkovLoop.cpp diff --git a/src/ArmaTDMarkovLoop.cpp b/src/ArmaTDMarkovLoop.cpp new file mode 100644 index 0000000..7fce9e3 --- /dev/null +++ b/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; + } From 23af0fa4f5fd69e42d27b8eb48314f30fcbbfdd5 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sun, 22 Nov 2020 16:44:50 +0000 Subject: [PATCH 2/3] also added some basic tests on TD markov. --- tests/testthat/testArmaTDMarkov.R | 88 +++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) create mode 100644 tests/testthat/testArmaTDMarkov.R diff --git a/tests/testthat/testArmaTDMarkov.R b/tests/testthat/testArmaTDMarkov.R new file mode 100644 index 0000000..a38c314 --- /dev/null +++ b/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. + + From 2da5b7bbef78feb01c1231f0b35bd48a0fc7eb03 Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sun, 22 Nov 2020 16:48:47 +0000 Subject: [PATCH 3/3] have left out Roxygen comments as not sure what you want to call this.. --- R/RcppExports.R | 4 ++++ src/RcppExports.cpp | 13 +++++++++++++ 2 files changed, 17 insertions(+) diff --git a/R/RcppExports.R b/R/RcppExports.R index 0d9b662..ebf36ca 100644 --- a/R/RcppExports.R +++ b/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 diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 08099be..e2af5c6 100644 --- a/src/RcppExports.cpp +++ b/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} };