1-SetupModel.Rmd

# Setup HBV model 
# ===============

Authors: Richard T. Gray, Neil Bretana

This script is used to create parameter sets for a project based on the 
input files and set-up scenario runs. 

```{r Setup}
# Clear workspace
rm(list = ls()) 

# Setup directories after setting working directory to source file 
# directory (so it is easy to move code around)
basePath <- getwd()
Rcode <- file.path(basePath, "code")

# Load useful libraries
source(file.path(Rcode, "LoadLibrary.R"))
LoadLibrary("readr")
LoadLibrary("dplyr")
LoadLibrary("tidyr")

# Project directory
project_directory <- file.path(basePath, "projects")
project_directory
```

## User inputs

The following chunk is where the user enters project details and
scenario specifications. 

```{r Project specifications}

# Project to load and setup
project_name <- "Australia_HBV_estimates-BM"

# Number of parameter sets to run
number_samples <- 10

# Current year to run future scenarios - simluation runs from the start
# of this year
current_year <- 2016

```

## Project set-up

```{r Update project specifications}

# Load current project specs and inputs
projectFolder <- file.path(project_directory, project_name)
load(file.path(projectFolder, paste0(project_name, ".rda")))

# Add some additional specs -----------------------------------------------

# Actual timesteps we will use for simulations - remove the final point
# so we go from begining of start_year to start of end_year. For 
# simulations will remove final point
pg$pts <- seq(pg$start_year, pg$end_year, by = pg$timestep)
pg$pts <- head(pg$pts, -1)
pg$npts <- length(pg$pts)

pg$parameter_samples <- number_samples

pg$current_year <- current_year

# Read in input files (which have been updated by the user)
initialPops <- read_csv(file.path(projectFolder, 
                                  "initial_populations.csv"))
constants <- read_csv(file.path(projectFolder,
                                "parameters_constants.csv"))
timeVarying <- read_csv(file.path(projectFolder,
                                  "parameters_time_varying.csv"))
timeFactors <- read_csv(file.path(projectFolder,
                                 "time_varying_ranges.csv"))
transitions <- read_csv(file.path(projectFolder,
                               "population_transitions.csv"))
waifw_matrix <- read_csv(file.path(projectFolder,
                                "population_waifw_matrix.csv"))

# Convert transistions and waifw_matrix into matrices
transitions <- as.matrix(transitions[, 2:(pg$npops + 1)])
rownames(transitions) <- colnames(transitions)

waifw_matrix <- as.matrix(waifw_matrix[, 2:(pg$npops + 1)])
rownames(waifw_matrix) <- colnames(waifw_matrix)

# Update inputs list
inputs <- list(initial_pops = initialPops, 
               constants = constants, time_varying = timeVarying,
               time_varying_ranges = timeFactors,
               transitions = transitions, waifw_matrix = waifw_matrix)

```
  
The following chunk contains the code to generate the best estimate 
parameter set and the ramdomly sampled parameter sets for the project.
  
```{r Generate parameter sets}

# Create function to sample from time varying factors and constants range
randomFactors <- function(param, paramFactors) { 
  paramValues <- filter(paramFactors, parameter == param)
  paramStart <- runif(1, paramValues$start_lower, 
                      paramValues$start_upper)
  paramEnd <- runif(1, paramValues$end_lower, 
                      paramValues$end_upper)
  
  return(seq(paramStart, paramEnd, length = pg$nyears))
}

randomParams <- function(param, constants) { 
  paramValues <- filter(constants, parameter == param)
  paramSample <- runif(1, paramValues$lower, 
                      paramValues$upper)
  return(paramSample)
} 

# Create a best estimate parameter set -----------------------------------
bestEstimateYears <- timeVarying

# Reshape constants so they can be appended to timevarying
constantsSpread <- constants %>%
  select(parameter, value)
parameters <- constantsSpread$parameter
constantsSpread <- as.data.frame(t(constantsSpread[,-1]))
colnames(constantsSpread) <- parameters
rownames(constantsSpread) <- NULL

# Append to time varyingCC
constantsDf <- constantsSpread[rep(1,nrow(timeVarying)),]
bestEstimateYears <- cbind(timeVarying, constantsDf)

# Expand to all points by linear extrapolation bewteen years
best_estimates <- as.data.frame(matrix(0, ncol = ncol(bestEstimateYears), 
                        nrow = pg$npts))
colnames(best_estimates) <- colnames(bestEstimateYears)

for (var in colnames(bestEstimateYears)) {
  tempValues <- bestEstimateYears[, var]
  for (year in 1:(pg$nyears-1)) {
    indices <- ((year-1)* 1/pg$timestep + 1): (year * 1/pg$timestep + 1)
    yearValues <- seq(tempValues[year], tempValues[year+1], 
                      length = (1 + 1/pg$timestep))

    # Note each time step is given the annual value 
    # This is adjusted in the model code
    best_estimates[indices, var] <- yearValues 
  }  
}

# Set up best estimate initial population
init_pop <- filter(initialPops, parameter == "init_pop")$value

initial_pop <- initialPops %>%
  filter(parameter != "init_pop")

initialPopMat <- as.data.frame(matrix(initial_pop$value, 
                                      ncol = pg$nstates + 1,  
                      nrow = pg$npops))
colnames(initialPopMat) <- c(pg$state_names, "pop_prop")
rownames(initialPopMat) <- pg$population_names

popProp <- init_pop * initialPopMat$pop_prop

best_initial_pop <- apply(initialPopMat[, 1:pg$nstates], 2, 
                          function(x) x * popProp)

# Create parameter set samples --------------------------------------------

param_sets <- list()
param_initial_pops <- list()

for (set in 1:number_samples) {

  paramSetYear <- bestEstimateYears
  
  # Loop through time varying and sample relative factors
  for (var in colnames(timeVarying)) {
    tempFactors <- randomFactors(var, timeFactors)
    tempParam <- tempFactors * timeVarying[, var]
    paramSetYear[, var] <- tempParam
  }
  
  # Loop through constants and sample
  for (var in colnames(constantsSpread)) {
    tempSample <- randomParams(var, constants)
    # tempParam <- tempSample * constantsSpread[, var]
    paramSetYear[, var] <- tempSample
  }
  
  # Expand to all points by linear extrapolation bewteen years
  paramSet <- as.data.frame(matrix(0, ncol = ncol(paramSetYear), 
                                        nrow = pg$npts))
  colnames(paramSet) <- colnames(paramSetYear)
  
  for (var in colnames(paramSetYear)) {
    tempValues <- paramSetYear[, var]
    for (year in 1:(pg$nyears-1)) {
      indices <- ((year-1)* 1/pg$timestep + 1): (year * 1/pg$timestep + 1)
      yearValues <- seq(tempValues[year], tempValues[year+1], 
                        length = (1 + 1/pg$timestep))
      paramSet[indices, var] <- yearValues
    }  
  }
  
  # Set up initial population size for each parameter set
  tempInitPop <- randomParams("init_pop", initialPops)
  
  tempInitialPop <- c()
  for (param in initial_pop$parameter) {
    paramSample <- randomParams(param, initialPops)
    tempInitialPop <- c(tempInitialPop, paramSample)
  }
  
  tempPopMat <- as.data.frame(matrix(tempInitialPop, 
                                        ncol = pg$nstates + 1,  
                        nrow = pg$npops))
  colnames(tempPopMat) <- c(pg$state_names, "pop_prop")
  rownames(tempPopMat) <- pg$population_names
  
  tempPopProp <- tempInitPop * tempPopMat$pop_prop
  
  paramPop <- apply(tempPopMat[, 1:pg$nstates], 2, 
                            function(x) x * tempPopProp)
  
  # Save final parameter set and initial populations in a list
  setStr <- paste0("paramSet", toString(set))
  param_sets[[setStr]] <- paramSet
  param_initial_pops[[setStr]] <- paramPop
  
}

```

The final chunk updates the project file for running in the model. 

```{r Save file project file}
# Re-save projects file with updated specs, inputs and append best
# estimates and parameter sets
save(pg, inputs, best_estimates, param_sets, best_initial_pop,
     param_initial_pops, transitions, waifw_matrix,
     file = file.path(projectFolder,
                      paste0(project_name, ".rda")))
```