# ae_functions.R
AE_outside_pr <-function(Rate_cancer.RDS,
                         Rate_cardio.RDS,
                         Rate_stroke.RDS,
                         life_exp_minus.RDS,
                         life_exp_smn.RDS,
                         cycles_inyear){
  
  # input: Rate_cancer.RDS created by AE-data-processing.Rmd ( in 01_data_raw)
  #        Rate_cardio.RDS created by AE-data-processing.Rmd ( in 01_data_raw)
  #        Rate_stroke.RDS created by AE-data-processing.Rmd ( in 01_data_raw)
  #         life_exp_minus.RDS created by AE-data-processing.Rmd ( in 01_data_raw)
  #        cycles_inyear: number of cycles in a year (must be min of 1)
  # output: return_mat a matrix, with the first collumn being seq(0, 110, by = cycle_inyear )
  #         the second collumn being probability of stroke during that cycle,
  #         3rd probaiblity of cardiovascular, and 4th probability of cancer event, 5th probability
  #         of outside risk of death
  
  
  require(dplyr)
  ## Currently all data is rate per 100000, between age X and X +1 
  ## Turning into probabilities
  ## First by dividing rate by 100000 to get a rate per person year
  ## Dividing by 12 to get a rate by person month
  ## Turning into a probability
  
  vect_time <- seq(0, 110, by =1/cycles_inyear)
  
  Pr_cancer <- Rate_cancer.RDS %>% 
    mutate(rate_person = Rate_100000/100000,
           rate_cylce = rate_person/cycles_inyear,
           pr_cycle  = 1-exp(-rate_cylce))
  
  Pr_cancer_vect <- vector(mode = "numeric", length = length(vect_time))
  
  for( i in 1:dim(Pr_cancer)[1]){
    Pr_cancer_vect[ vect_time >= Pr_cancer$AgeStart[i] & vect_time <Pr_cancer$AgeEnd[i]] <- Pr_cancer$pr_cycle[i]
  }
  
  pr_cardio <- Rate_cardio.RDS %>% 
    mutate(AgeEnd = Age_group+1,
           AgeStart = Age_group,
           rate_person = Rate/100000,
           rate_cylce = rate_person/cycles_inyear,
           pr_cycle  = 1-exp(-rate_cylce)
    )
  
  Pr_cardio_vect <- vector(mode = "numeric", length = length(vect_time))
  
  for( i in 1:dim(pr_cardio)[1]){
    Pr_cardio_vect[ vect_time >= pr_cardio$AgeStart[i] & vect_time <pr_cardio$AgeEnd[i]] <- pr_cardio$pr_cycle[i]
  }
  
  
  pr_stroke <- Rate_stroke.RDS %>% 
    mutate(AgeEnd = AgeEnd_exclusive, 
           rate_person = Rate/100000,
           rate_cylce = rate_person/cycles_inyear,
           pr_cycle  = 1-exp(-rate_cylce))
  
  Pr_stroke_vect <- vector(mode = "numeric", length = length(vect_time))
  
  for( i in 1:dim(pr_stroke)[1]){
    Pr_stroke_vect[ vect_time >= pr_stroke$AgeStart[i] & vect_time < pr_stroke$AgeEnd[i]] <- pr_stroke$pr_cycle[i]
  }
  
  
 pr_deathoutside <-  life_exp_minus.RDS %>% 
    mutate(AgeEnd = Ending_Age,
           AgeStart = Starting_Age,
           rate_person = Allcause_minus/100000,
           rate_cylce = rate_person/cycles_inyear,
           pr_cycle  = 1-exp(-rate_cylce))
 Pr_death_vect <- vector(mode = "numeric", length = length(vect_time))
 
 for( i in 1:dim(pr_deathoutside)[1]){
   Pr_death_vect[ vect_time >= pr_deathoutside$AgeStart[i] & vect_time < pr_deathoutside$AgeEnd[i]] <- pr_deathoutside$pr_cycle[i]
 }
 
 
 pr_deathoutside_SMN <- life_exp_smn.RDS %>% 
   mutate(AgeEnd = Ending_Age,
          AgeStart = Starting_Age,
          rate_person = Rate_1000000/100000,
          rate_cycle = rate_person/cycles_inyear,
          pr_cycle = 1-exp(-rate_cycle))
 
 Pr_death_smn <- vector(mode = "numeric", length = length(vect_time))
 
 for( i in 1:dim(pr_deathoutside_SMN)[1]){
   Pr_death_smn[ vect_time >= pr_deathoutside_SMN$AgeStart[i] & vect_time < pr_deathoutside_SMN$AgeEnd[i]] <- pr_deathoutside_SMN$pr_cycle[i]
 } 
 
  
  return_mat  <-  cbind(vect_time,Pr_stroke_vect,Pr_cardio_vect,Pr_cancer_vect,Pr_death_vect,Pr_death_smn)
  return(return_mat)
}
# AE_outside_pr(Rate_cancer.RDS = readRDS(here::here("02_data", "Rate-cancer.RDS")),
#              Rate_cardio.RDS = readRDS(here::here("02_data", "Rate-cardio.RDS")),
#              Rate_stroke.RDS = readRDS(here::here("02_data", "Rate-stroke.RDS")),
#             life_exp_minus.RDS = readRDS(here::here("02_data","life-exp-minus.RDS")),
#               life_exp_smn.RDS =  readRDS(here::here("02_data","life-exp-smn.RDS"))
#              cycles_inyear = 12)


### SRS_gen
# Creating a function that creates an SRS vector as a function of cycle length

RSR_gen <- function(cycles_year,SRS_small ){
  # input: cycles_year: number of cycles in a year
  #        SRS_small: RDS generated by AE-data-processing.Rmd 
  # output: return_mat a matrix with two collumns, (time since diagnosis and SRS linearly extrapolated)
  t2 <- SRS_small %>% 
    mutate(TimeFromDx2 = lead(TimeFromDx),
           SRS2 = lead(SRS)) %>% 
    na.omit()
  
  
  time_vect <- vector(mode = "list")
  SRS_vect <- vector(mode = "list")
  for(i in 1:dim(t2)[1]){
    time_vect_temp <-   c(seq(from = t2$TimeFromDx[i],to = t2$TimeFromDx2[i], by = 1/cycles_year))
    SRS_vect_temp <-   seq(from = t2$SRS[i],to = t2$SRS2[i], length.out = length(time_vect_temp) )
    
    time_vect[[i]] <-time_vect_temp 
    SRS_vect[[i]] <- SRS_vect_temp 
  }
  
  return_mat <- cbind(time_vect = unlist(time_vect),SRS_vect=unlist(SRS_vect))
  return_mat <- return_mat[-c(which(return_mat[,1]-lead(return_mat[,1]) == 0)),]
  
  
 # Creating SRS(t + 1)/STS(t) 
  return_mat <-cbind(return_mat, SRS_dif = c(1, return_mat[-1,2]/return_mat[-dim(return_mat)[1],2]))

  return(return_mat)
}

# SRS_gen(cycles_year = 52,SRS_small = readRDS(here::here("02_data","SRS_cancer.RDS")))

# PrDeath_SM probability of death for secondary neoplasm
#
PrDeath_SM <- function(CurAge, VectSRSdif,MortRiskAge){
  # CurAge: a vector or numeric indicating the current age (floor) of those with cancer
  # MortRiskAge: Probability of dying from non SMN reasons for current Age
  # VectSRSdif: output SRS ratios , indexed from time from dx
  # output:   PrMortOnc (probability of stroke mortality)
  
  
  # if time from diagnosis assum that the SRS ratio is SRS(10)/SRS(10-cycle_length) 

  PrMortOnc <- 1-VectSRSdif*(1-MortRiskAge)
  return(PrMortOnc)
}

#PrDeath_SM(CurAge =  c(60,61,87),TimeFromDx =c(1,9,6),VectSRSdif =  SRS_mat[fun_index(TimeFromDx, cycles_year = 12), "SRS_dif"] , 
#           AE_outside_mat[fun_index(CurAge, cycles_year = 12) ,"Pr_death_smn"] )


# Sourcing a 1 yr mortality probability from Saposnik, Gustavo, et al. "IScore: a risk score to predict death early after hospitalization for an acute ischemic stroke." Circulation 123.7 (2011): 739-749.
# Assuming no mortality after year 1. Function of covaraites (Age + Male*5 ) + Baseline Average Score from population

PrDeath_Stroke <- function(AgeStroke,Gender, cycles_year){
  # input Agestroke (round)
  # Gender (1 is equal to male)
  # cycles_year: number of cycles in a year
  
  PrMort <- vector(length = length(AgeStroke),mode = "numeric")  
  MortScore  <- AgeStroke + (Gender == 1)*5 + 41.765
  
  PrMort[MortScore < 59 ] <-  0  
  PrMort[MortScore >= 59 & MortScore  <	71	] <-  0.0178
  PrMort[MortScore >= 71  & MortScore <	81	] <-  0.0293
  PrMort[MortScore >= 81  & MortScore <	91	] <-  0.0436
  PrMort[MortScore >= 91  & MortScore <	101	] <-  0.0656
  PrMort[MortScore >= 101 & MortScore <	111	] <-  0.0976
  PrMort[MortScore >= 111 & MortScore <	121	] <-  0.144
  PrMort[MortScore >= 121 & MortScore <	131	] <-  0.204
  PrMort[MortScore >= 131 & MortScore <	141	] <-  0.29
  PrMort[MortScore >= 141 & MortScore <	151	] <-  0.382
  PrMort[MortScore >= 151 & MortScore <	161	] <-  0.489
  PrMort[MortScore >= 161 & MortScore <	171	] <-  0.593
  PrMort[MortScore >= 171 & MortScore <	181	] <-  0.748
  PrMort[MortScore >= 181 & MortScore <	191	] <-  0.757
  PrMort[MortScore >= 191 & MortScore <	201	] <-  0.873
  PrMort[MortScore >= 201 & MortScore <	211	] <-  0.938
  PrMort[MortScore >= 211 & MortScore <	221	] <-  0.953
  PrMort[MortScore >= 221 & MortScore <	231	] <-  0.973
  PrMort[MortScore >= 231 & MortScore <	240	] <-  0.977
  PrMort[MortScore >= 240 ] <-  1 
  
  # turning yearly probaibliteis into cycle probabilities
  
  PrRate <- -log(1-PrMort)/1
  PrRate <-  PrRate/cycles_year
  ProbMort <-  1-exp(-PrRate)
  return(ProbMort)
}



# PrDeath_Stroke(AgeStroke = c(10,15,80),Gender = c(1,1,0),cycles_year = 12)

MortCardio_est <- function(est_par, cycleinyear){
  # input:
  #est_par: output from CardioMortPar.RDS
  #n umber of cycles in a year
  
  time_vect <- seq(0,100, by = 1/cycleinyear)
  est_mort <- 1-pweibull( q=time_vect, est_par$par1,est_par$par2)
  PrMortCardio <- 1-est_mort[-1]/est_mort[-length(est_mort)]
  PrMortCardio <- c(0,PrMortCardio)
  
  PrMortCardio_mat <- cbind(time_vect, PrMortCardio)
  return(PrMortCardio_mat)
}