---
title: "Examples on Simulating Migration Flows for the MicSim Package"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Examples on Simulating Migration Flows for the MicSim Package}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

The sample data on migration between Spain, Sweden and the Netherlands were prepared by Claudio Bosco (European Commission), Daniela Ghio (European Commission), Maurizio Teobaldelli (European Commission) and Sabine Zinn (German Socio-Economic Panel, Humboldt University Berlin). EUROSTAT data were used as the data source. The sample was intentionally kept small. MicSim can also handle larger numbers of cases, but to be efficient in terms of run times, this requires a bit more computing power with more CPUs. 
```{r}
# Load library
library(MicSim)
```

### Definition of Basic Simulation Frame 
```{r}
# Defining simulation horizon
startDate <- 20140101 # yyyymmdd
endDate   <- 20181231 # yyyymmdd
simHorizon <- c(startDate=startDate, endDate=endDate)

# Seed for random number generator
set.seed(12)

# Definition of maximal age (sharp, i.e. max age is 100.0)
maxAge <- 100

# Definition of state space
sex <- c("m","f")                     
country_R <- c("NL","ES","SE")
fert <- c("0","1","2","3","4+","999")           
stateSpace <- expand.grid(sex=sex,country_R=country_R,fert=fert)

# Definition of nonabsorbing and absorbing states
absStates <- c("dead","rest")   
```

### Definition of initial population and newborn characteristics
```{r setup}
# initial population included in the MicSim package
initPop <- initPopMigrExp
head(initPop)

# population of immigrants entering to virtual population over simulation time
immigrPop <- immigrPopMigrExp # included in the MicSim package
head(immigrPop)

# Definition of initial states for newborns
fixInitStates <- 2 # give indices for attribute/substate that will be taken over from the mother, here: nat (nationality)
varInitStates <- rbind(c("m","ES","0"), c("f","ES","0"), # to have possibility to define distinct sex ratios in distinct countries, 
                       c("m","NL","0"), c("f","NL","0"), # hold substate that are inherited by mother in the init state (i.e. nationality)
                       c("m","SE","0"), c("f","SE","0")) 
initStatesProb <- c(0.5151,0.4849, # probabilities for female / male newborns for each nationality separately
                    0.5124,0.4876,
                    0.5146,0.4854)
```

### Rates Definition
Beware: Rates have to be given at least for age [0,maxAge) and for all years within the simulation horizon.
At this the exact value of *maxAge* is excluded, i.e. here 100.00 but not e.g. age=99.9999. 
```{r}
# Load rates from data included in the MicSim package (one column for one transition)
rates <- migrExpRates

# Define function to easily transform rates in function format required by MicSim 
require(glue)
for (i in 1:length(names(rates))) {
  script_var = names(rates[i])
  eval(parse(text = glue("{script_var} <- unlist(as.numeric(rates[,{i}]))")))
}

# --------------------------------------------------------------------------
# Fertility Rates
# --------------------------------------------------------------------------

fertRates_NL_0_1 <- function(age,calTime, duration){
  rate <- fert_NL_0_1[as.integer(age)+1]
  return(rate)  
}

fertRates_NL_1_2 <- function(age,calTime, duration){
  rate <- fert_NL_1_2[as.integer(age)+1]
  return(rate)  
}

fertRates_NL_2_3 <- function(age,calTime, duration){
  rate <- fert_NL_2_3[as.integer(age)+1]
  return(rate)  
}

fertRates_NL_3_4 <- function(age,calTime, duration){
  rate <- fert_NL_3_4[as.integer(age)-+1]
  return(rate)  
}

fertRates_ES_0_1 <- function(age,calTime, duration){
  rate <- fert_ES_0_1[as.integer(age)+1]
  return(rate)  
}

fertRates_ES_1_2 <- function(age,calTime, duration){
  rate <- fert_ES_1_2[as.integer(age)+1]
  return(rate)  
}

fertRates_ES_2_3 <- function(age,calTime, duration){
  rate <- fert_ES_2_3[as.integer(age)+1]
  return(rate)  
}

fertRates_ES_3_4 <- function(age,calTime, duration){
  rate <- fert_ES_3_4[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_0_1 <- function(age,calTime, duration){
  rate <- fert_SE_0_1[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_1_2 <- function(age,calTime, duration){
  rate <- fert_SE_1_2[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_2_3 <- function(age,calTime, duration){
  rate <- fert_SE_2_3[as.integer(age)+1]
  return(rate)  
}

fertRates_SE_3_4 <- function(age,calTime, duration){
  rate <- fert_SE_3_4[as.integer(age)+1]
  return(rate)  
}

# --------------------------------------------------------------------------
# Internal Migration Rates
# --------------------------------------------------------------------------

`%!in%` <- Negate(`%in%`)

for(i in 1:length(country_R)) {
  other_provinces = country_R[which(country_R %!in% country_R[i])]
  for(k in 1:length(other_provinces)) {
    eq = paste(sprintf('%s_%s_rates', glue("{country_R[i]}"),  glue("{other_provinces[k]}")), 
               '<- function(age,calTime, duration)', '{',
               sprintf('rate <- rate_%s_%s[as.integer(age)+1]', glue("{country_R[i]}"), 
               glue("{other_provinces[k]}")), "\n ", 'return(rate)','}')
    eval(parse(text = eq))
  }
}

# --------------------------------------------------------------------------
# Mortality Rates
# --------------------------------------------------------------------------

# Female mortality
for(i in 1:length(country_R)) {
  eq = paste(sprintf('mortRates_f_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- mort_f_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

# Male mortality
for(i in 1:length(country_R)) {
  eq = paste(sprintf('mortRates_m_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- mort_m_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

# ---------------------------------------------------------------------------
# Emigration rates
# ---------------------------------------------------------------------------

# Emigration rates for females
for(i in 1:length(country_R)) {
  eq = paste(sprintf('emigrRates_f_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- emig_f_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

# Emigration rates for males
for(i in 1:length(country_R)) {
  eq = paste(sprintf('emigrRates_m_%s', glue("{country_R[i]}")), '<- function(age,calTime, duration)',
             '{',
             sprintf('rate <- emig_m_%s[as.integer(age)+1]', glue("{country_R[i]}")),
             "\n ",
             'return(rate)',
             '}')
  eval(parse(text = eq))
}

```

### Transition pattern and assignment of functions specifying transition rates
```{r}
# ---------------------------------------------------------------------------
# Transition matrix for fertility
# ---------------------------------------------------------------------------
fertTrMatrix <- cbind(c("f/ES/0->f/ES/1","f/ES/1->f/ES/2","f/ES/2->f/ES/3","f/ES/3->f/ES/4+",
                        "f/SE/0->f/SE/1","f/SE/1->f/SE/2","f/SE/2->f/SE/3","f/SE/3->f/SE/4+",
                        "f/NL/0->f/NL/1","f/NL/1->f/NL/2","f/NL/2->f/NL/3","f/NL/3->f/NL/4+"),
                      c("fertRates_ES_0_1", "fertRates_ES_1_2", "fertRates_ES_2_3", "fertRates_ES_3_4",
                        "fertRates_SE_0_1", "fertRates_SE_1_2", "fertRates_SE_2_3", "fertRates_SE_3_4",
                        "fertRates_NL_0_1", "fertRates_NL_1_2", "fertRates_NL_2_3", "fertRates_NL_3_4"))

# ---------------------------------------------------------------------------
# Transition matrix for migration
# ---------------------------------------------------------------------------
# First: make names for transition matrix, i.e. stateOfOrigin->stateOfDestination 
testo1 <- "intmigTrMatrix <- cbind(c("
for(i in 1:length(country_R)) {
  for(m in 1:length(country_R)) {
    if(m != i) {
      eq1 = paste(sprintf("'%s->%s'", glue("{country_R[i]}"), glue("{country_R[m]}")))
      if (i == length(country_R) & m == (i-1)) {
        testo1 = paste(testo1,eq1)
      } else {
        testo1 = paste(testo1 ,eq1, ",")
      }
    }
    if (m == i & m == length(country_R)){
      testo1 = paste(testo1, "),")
      
    }
  }
}

#Second: set names for transition functions
testo1 = paste (testo1,"c(")
for(i in 1:length(country_R)) {
  for(m in 1:length(country_R)) {
    if(m != i) {
      eq1 = paste(sprintf("'%s_%s_rates'", glue("{country_R[i]}"), glue("{country_R[m]}")))
      if (i == length(country_R) & m == (i-1)) {
        testo1 = paste(testo1,eq1)
      } else {
        testo1 = paste(testo1 ,eq1, ",")
      }
    }
    if (m == i & m == length(country_R)){
      testo1 = paste(testo1, "))")
    }
  }
}
eval(parse(text = testo1))

# ---------------------------------------------------------------------------
# Transition matrix for mortality and emigration
# ---------------------------------------------------------------------------

testo <- "absTransitions <- rbind("
for(i in 1:length(country_R)) {
  for(m in 1:length(sex)) {
    eq1 = paste(sprintf("c('%s/%s/dead', 'mortRates_%s_%s')", glue("{sex[m]}"), glue("{country_R[i]}") ,
                        glue("{sex[m]}"), glue("{country_R[i]}")),',',
                sprintf("c('%s/%s/rest', 'emigrRates_%s_%s')", glue("{sex[m]}"),glue("{country_R[i]}") ,
                        glue("{sex[m]}"), glue("{country_R[i]}")))
    if(i == length(country_R) & m == length(sex)) {
      testo = paste(testo, eq1, ")")
    } else {
      testo = paste(testo,eq1, ",")
    }
  }
}
eval(parse(text = testo))

# ---------------------------------------------------------------------------
# Combine all
# ---------------------------------------------------------------------------

allTransitions <- rbind(fertTrMatrix, intmigTrMatrix)
transitionMatrix <- buildTransitionMatrix(allTransitions=allTransitions,
                                          absTransitions=absTransitions, 
                                          stateSpace=stateSpace)

# ---------------------------------------------------------------------------
# Define transitions triggering a birth event
# ---------------------------------------------------------------------------

fertTr <- fertTrMatrix[,1]
```

### Run Simulation (using one core)
For illustration purpose, the subsequent run is limited to the first 500 people and to 100 migrants.
However, just remove the restriction to run the whole sample, i.e. 
using *initPop* instead of *initPop[1:500,]* and *immigrPop* instead of *immigrPop[1:100,]*.

```{r}
pop <- micSim(initPop=initPop[1:500,], immigrPop=immigrPop[1:100,],
              transitionMatrix=transitionMatrix, absStates=absStates,
              varInitStates=varInitStates, initStatesProb=initStatesProb,
              fixInitStates=fixInitStates,
              maxAge=maxAge, simHorizon=simHorizon,fertTr=fertTr)
head(pop)

```

### Convert to long format

```{r}
popLong <- convertToLongFormat(pop, migr=TRUE)
head(popLong)

```

### Convert to wide format

```{r}
popWide <- convertToWideFormat(pop)
head(popWide)
```

### Run Simulation with parallel computing on several cores
Try this to speed up your simulation run. The more cores you can use the faster the simulation will be executed.
This example uses three cores. Give as many seeds as you use cores. That way you can replicate your results.
```{r}

cores <- 3
seeds <- c(34,145,97)
pop <- micSimParallel(initPop=initPop, immigrPop=immigrPop,
                      transitionMatrix=transitionMatrix, absStates=absStates,
                      varInitStates=varInitStates, initStatesProb=initStatesProb,
                      fixInitStates=fixInitStates,
                      maxAge=maxAge, simHorizon=simHorizon,fertTr=fertTr,
                      cores=cores, seeds=seeds)
head(pop)
```