params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )

params = 
  list(
    t_names = c("lifetable"),         # Strategy names. 
    n_treatments = 1,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 0,                  # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt                   # Processed HMD life table data (2019, USA)
  )

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )

params = 
    list(
        t_names = c("lifetable"),         # Strategy names. 
        n_treatments = 1,                 # Number of treatments
        s_names  = c("Alive", "Dead"),    # State names
        n_states = 2,                     # Number of states
        n_cohort = 100000,                # Cohort size
        initial_age = 0,                  # Cohort starting age    
        n_cycles = 110,                   # Number of cycles in model.  
        cycle = 1,                        # Cycle length
        life_table = lt                   # Processed HMD life table data (2019, USA)
    )

sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}sim_cohort <- function(params) {
 params$t_names %>% map(~({ 
    tr_ <- t(c("alive" = params$n_cohort, "dead" = 0))
    
    res <- do.call(rbind,lapply(params$mP, function(tp) {
        tr_ <<- tr_ %*% matrix(unlist(tp[,,.x]),nrow=params$n_state)
    }))
    res <- rbind(c(params$n_cohort,0),res) 
    dimnames(res) <- list(paste0(c(0:params$n_cycles)), params$s_names)
    res
  })) %>% 
    set_names(params$t_names)
}

trace <- 
  sim_cohort(params)
trace[["lifetable"]][1:10,]

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = (trace[[1]] * (1 / params$n_cohort)) %*% payoff_life_exp

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = cycle_adj(life_exp_cycle)

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

gom_fit <- MortalityLaw(
    x  = ages,      # vector with ages
    Dx  = deaths,   # vector with death counts
    Ex  = exposure, # vector containing exposures
    law = "gompertz",
    opt.method = "LF2")

# Fitted coefficients
gom_fit$coefficients

# Fitted coefficients
gom_fit$coefficients

# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Fitted coefficients
gom_fit$coefficients

# Hazard rate of death for a 40 year old based on Gompertz model
# mx = A exp(Bx)
gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Can simply use the supplied function to calcualte 
MortalityLaws::gompertz(x = 40, par = gom_fit$coefficients)

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)

# Calculate death hazard rate at age 40 from fitted Gompertz model. 
mx = gom_fit$coefficients["A"] * exp(gom_fit$coefficients["B"] * 40 )

# Convert mortality rate to probability
-log(1-mx)

# Compare with life table probability (qx) at age 40:
lt[41,c("age","qx")] 

plot(gom_fit)

ages     <- lt$age[lt$age>=40 & lt$age<100]
deaths   <- lt$dx[lt$age>=40 & lt$age<100]
exposure <- lt$lx[lt$age>=40 & lt$age<100]

gom_fit40 <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "gompertz",
                opt.method = "LF2")

plot(gom_fit40)

ages     <- lt$age[lt$age<100]
deaths   <- lt$dx[lt$age<100]
exposure <- lt$lx[lt$age<100]

hp_fit <- MortalityLaw(
                x  = ages,      # vector with ages
                Dx  = deaths,   # vector with death counts
                Ex  = exposure, # vector containing exposures
                law = "HP",
                opt.method = "LF2")

plot(hp_fit)

HP(x = 40, par = hp_fit$coefficients)

x   = 40
mu1 = A^((x + B)^C) + G * H^x)
mu2 = D * exp(-E * (log(x/F_))^2))
eta = ifelse(x == 0, mu1, mu1 + mu2)
hx  = eta/(1 + eta)

gompertz(x = 40, par = gom_fit$coefficients)$hx

HP(x = 40, par = hp_fit$coefficients)$hx

u = 0.25

u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])

u = 0.25
qgompertz(p = u, shape = gom_fit$coefficients["B"], rate = gom_fit$coefficients["A"])

qHP <- 
  approxfun(
    y = 0:115,
    x = 1 - exp(-cumsum(HP(x = 0:115,
                           par = hp_fit$coefficients)$hx)))

qHP <- 
  approxfun(y = 0:115,
            x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(0.25)

# Sample death times for 10 individuals
u <- runif(10, min = 0, max = 1)
qHP <- 
  approxfun(y = 0:115,
  x = 1 - exp(-cumsum(HP(x = 0:115, par = hp_fit$coefficients)$hx)))
qHP(u)

params = 
  list(
    t_names = c("lifetable",          # Strategy names
                "heligman-pollard",
                "gompertz"),          
    n_treatments = 3,                 # Number of treatments
    s_names  = c("Alive", "Dead"),    # State names
    n_states = 2,                     # Number of states
    n_cohort = 100000,                # Cohort size
    initial_age = 40,                 # Cohort starting age    
    n_cycles = 110,                   # Number of cycles in model.  
    cycle = 1,                        # Cycle length
    life_table = lt,                  # Processed HMD life table data (2019, USA)
    gom_fit = gom_fit,                # Fitted Gompertz model object
    hp_fit = hp_fit                   # Fitted Heligman-Pollard model object
  )

# Add transition matrices into params.
params$mP <- 
  fn_mPt(1:params$n_cycles, params)

# Construct the Markov trace. 
trace <- 
  sim_cohort(params)

# Define the payoff for life expectancy
payoff_life_exp = c("Alive" = 1, "Dead" = 0)

# Calculate alive years in each cycle
life_exp_cycle = lapply(trace, function(tr) (tr*(1 / params$n_cohort)) %*%
                 payoff_life_exp)

# Integration via alternative Simpson's method (cycle correction)
alt_simp_coef <- function(i) c(17, 59, 43, 49, rep(48, i-8),
                               49, 43, 59, 17) / 48
cycle_adj     <- function(x) sum(alt_simp_coef(length(x)) * x)

total_life_exp = lapply(life_exp_cycle,cycle_adj)

total_life_exp

ihme_cvd <- 
  tibble::tribble(
        ~age_name,        ~val,
               1, 0.038771524,
               5, 0.038546046,
              10, 0.044403585,
              15, 0.033781126,
              20, 0.035856165,
              25, 0.053077797,
              30, 0.086001439,
              35, 0.130326551,
              40, 0.184310334,
              45,  0.21839762,
              50, 0.243705394,
              55, 0.256334637,
              60,  0.26828001,
              65, 0.272698709,
              70,  0.28529754,
              75, 0.310642009,
               0, 0.016750489,
              80, 0.353518012,
              85, 0.399856716,
              90, 0.447817792,
              95, 0.495305502
        ) %>% 
    mutate(age_ihme = cut(age_name,unique(c(0,1,seq(0,95,5),105)),right=FALSE))  %>% 
    select(age_ihme,  pct_cvd = val) 

lt_ <- 
  lt_ %>% 
    # Cause specific probability of death. 
    mutate(qi = q * dx_i / dx) %>% 
    mutate(Rd = (dx - dx_i) / dx) %>% 
    mutate(md = m * Rd) %>% 
    mutate(mi = m - md) %>% 
    mutate(dxd = dx - dx_i)

ages_     <- lt_$age[lt_$age<100 & lt_$age>=25]
deaths_   <- lt_$dxd[lt_$age<100 & lt_$age>=25] 
exposure_  <- lt_$lx[lt_$age<100 & lt_$age>=25]

mort_fit_CVDdeleted <- MortalityLaw(
    x  = ages_,
    Dx  = deaths_,   # vector with death counts
    Ex  = exposure_, # vector containing exposures
    law = "HP2",
    opt.method = "LF2")

trace <- 
  sim_cohort(params)
trace[["natural_history"]][1:10,]