Run Model from R

OpenM++ integration with R

It is a convenient to use GNU R to prepare model parameters and analyze output values. There are two different R APIs which we can use for openM++ models:

• openMpp package: simple and convenient specially for desktop users, upstream and downstream analysis;
• oms JSON web-service API: preferable choice to run models on computational clusters and in cloud.

Below is a simple loop example to run NewCaseBased model on desktop using openMpp R package. There is similar example how to run model in cloud and save results in CSV file using oms JSON web-service.

OpenM++ provides R package openMpp to simplify access to openM++ database for R developers. To find out more about openMpp R package please check:

• openMpp package documentation
• installation ReadMe and source code

There is also an excelent R package created by Matthew T. Warkentin available at: mattwarkentin.github.io/openmpp.

Following R example is running openM++ "NewCaseBased" test model with 16 subsamples using mortality hazard data:

mortalityData <- data.frame( 
  value = seq(from = 0.014, by = 0.005, length.out = 20)
)

As result Mortality Hazard increases about eight times in the range of [0.014, 0.109] and we can see eight time decrease of Duration of Life from initial 72 years down to 9 years.

R script

# use openMpp library for openM++ database access
library(DBI)
library("openMpp")
library("RSQLite")

#
# R integration example using NewCaseBased model
#   loop over MortalityHazard parameter 
#   to analyze DurationOfLife
#

##################################################################
# To run this example please uncomment and/or change values below
# to match your hardware and file system environment:
#
# model_exe    <- path to the model executable, i.e.: "./NewCaseBased" or "NewCaseBased.exe"
# model_sqlite <- path to the model.sqlite database:  "NewCaseBased.sqlite"
# model_args   <- optional arguments to control model run, for example:
#       -OpenM.SubValues 16 <- number of simation members
#       -OpenM.Threads 4    <- number of computatinal threads
#
### For running on a local machine using the working directory in R 
#
# For the following values to work, you must first set the R Working directory
# to the directory containing the NewCaseBased executable and the SQLite database.
# In RStudio Session > Set Working Directory > Choose Directory, 
# then navigate to location, e.g.: /OM_ROOT/models/NewCaseBased/ompp/bin
# Alternatively, one may use setwd(), e.g.: setwd("/OM_ROOT/models/NewCaseBased/ompp/bin")
#
model_exe = "./NewCaseBased"
model_sqlite = "NewCaseBased.sqlite"
model_args = " -OpenM.SubValues 16 -OpenM.Threads 4" 
# model_args = " "  # default: 1 simulation member and 1 thread
#
### For running on a local machine using explicit paths
#
# model_exe = "/path/to/executable/model/NewCaseBased"
# model_sqlite = "/path/to/SQLite/database/NewCaseBased.sqlite"
#
### For running on cluster (change to match your cluster)
#
# model_exe = "mpiexec"
# model_sqlite = "/mirror/NewCaseBased.sqlite"
# model_args = "-n 8 /mirror/NewCaseBased -OpenM.SubValues 16 -OpenM.Threads 2"
##################################################################

#
# NewCaseBased model parameters:
#   Mortality hazard: double number
#   Simulation cases: number of simulation cases
#   Simulation seed:  random seed
#
MortalityHazard <- list(
  
  name = "MortalityHazard",   # model parameter name
  value = 0.014,              # value of parameter
  
  txt = data.frame(
    lang = c("EN", "FR"),
    note = c("An arbitrarily selected value, chosen to produce a life expectancy of about 70 years", NA),
    stringsAsFactors = FALSE
  )
)

SimulationCases <- list(name = "SimulationCases", value = 5000L)
SimulationSeed <- list(name = "SimulationSeed", value = 16807)

#
# name, description and notes for this set of model parameters
#
inputSet <- data.frame(
  name = "LifeVsMortality",
  lang = c("EN", "FR"),
  descr = c("NewCaseBased working set of parameters", "(FR) NewCaseBased working set of parameters"),
  note = c(NA, NA),
  stringsAsFactors = FALSE
)

# 
# connect to database and find NewCaseBased model
#
theDb <- dbConnect(RSQLite::SQLite(), model_sqlite, synchronous = "full")
invisible(dbGetQuery(theDb, "PRAGMA busy_timeout = 86400"))   # recommended

defRs <- getModel(theDb, "NewCaseBased")      # find NewCaseBased model in database

# create new working set of model parameters based on existing model run results
#
firstRunId <- getFirstRunId(theDb, defRs)

setId <- createWorksetBasedOnRun(
  theDb, defRs, firstRunId, inputSet, 
  MortalityHazard, SimulationCases, SimulationSeed
  )
if (setId <= 0L) stop("workset creation failed")

setReadonlyWorkset(theDb, defRs, TRUE, setId)  # workset must be read-only to run the model

#
# analyze NewCaseBased model varying mortality hazard values
#
mortalityData <- data.frame( 
  value = seq(from = 0.014, by = 0.005, length.out = 20)
)

for (mortality in mortalityData$value)
{
  print(c("Mortality hazard: ", mortality))

  system2(
    model_exe,
    paste(
      model_args,
      " -Parameter.MortalityHazard ", mortality,
      " -OpenM.SetId ", setId, 
      " -OpenM.LogToConsole false",
      " -OpenM.LogToFile true",
      " -OpenM.ProgressPercent 100",
      sep = ""
    )
  )
}

#
# read final results from database
#   average duration of life: DurationOfLife.Expr3
#
runIdRs <- getWorksetRunIds(theDb, setId)   # get result id's

lifeDurationData <- NULL
for (runId in runIdRs$run_id)
{
  lifeDurationData <- rbind(
    lifeDurationData,
    selectRunOutputValue(theDb, defRs, runId, "DurationOfLife", "Expr3")
  )
}

dbDisconnect(theDb)   # close database connection

#
# display the results
#
plot(
  mortalityData$value, 
  lifeDurationData$expr_value, 
  type = "o",
  xlab = "Mortality Hazard", 
  ylab = "Duration of Life", 
  col = "red"
)