#############Executable file ######################################
#########################################
### Scenario: We need to choose the treatment scenario that we want to apply 
### Simulation starts at time of HCV infection
### This file include all the required packages, and call all the required sources
### In this file the Cohort Size and the duration of simulation are defined. 
#####################################################################
cat("\014")  
cat('\n')
cat('***********************************************************\n')
cat('Mathematical Modeling of Hepatitis C Virus Progression Disease\n')
cat('\t \t 2017')
cat('\n')
cat ('\n \t Version: 14')
cat('\n Requires files: \n')
cat('simulate_cohort.r: To do the simulation \n')
cat('baseline.r: Identify the baseline chracteristics \n')
cat('hazard_fun.r: Calls file specifying the hazard functions \n')
cat('params_cov.r: To choose parameters file \n')
cat('scenario_treatFx.r: To choose the treatment scenario \n')
cat('simulate_cohort.r: To do the simulation \n')
cat('TransitionTime.r: Shows all the possible transitions\n')
cat('***********************************************************\n')
#############################################################
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
rm(list =  ls())
wd =  getwd()
######################### Required Packages #################
require(MASS)
require(msm)
require(mstate)
require(survival)
require(splines)
require(plyr)
require(gems)
#########################################################
start.time<-Sys.time()
########################################################
cohortSize <- 100 # specify how many patients to simulate
maxTime<-100 # specify duration of simulation in years
#data <- read.table("baseline_data4.txt", header=T)
setwd("../Required_Data")
data <- read.table("baseline_data5.txt", header=T)
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")
case = 1 # the case study
#####################Source Codes ###########################
source("baseline.r") # call file specifying baseline characteristics
setwd("~/HCV_project_2017/hcv_git/hcv_modeling")

cat('***********************************************************\n')
ans <- readline(prompt="Please enter the treatment scenario you are interested to study: Treating patients in (and after) stage F0, F1, F2, F3, F4 or BL (baseline) \t")
source(sprintf("scenario_treat%s.r",ans))
#source("scenario_treatF3.r")
source("hazard_fun.r") # call file specifying the hazard functions
source("params_cov.r") #choose parameters file
source("TransitionTime.r")
source("sim_cohort.r") # call the cohort simulator
#######################################################
end.time <- Sys.time()
print("The simulation duration time: ")
print (end.time - start.time)
cat('***********************************************************\n')
tt <- readline(prompt="Do you want to save the results as a table (y/n)?")
if (tt == "y")
{
  setwd("../Outcomes")
ans <- readline(prompt="Please Enter a File Name:")
write.table(cohort@time.to.state,ans) # save the cohort
Total = cbind(cohort@time.to.state, bl)
write.table(Total, sprintf("Total_results_%s.txt",ans))
}
print("Do you want to do the post processing?")
answer <- readline(prompt="Do you want to do the post processing(y/n)?")
if (answer == "y")
{
source("outcomes.r")
source("out_bc.R")
source("barPlot.r")
source("out_baseline.r")
}

########################################################
# Age_in_2017 = 2017 - floor(bl[,"BirthYear"])
# hist(Age_in_2017, col="blue",xlab="Age in 2017")
hist(floor(bl[,"BirthYear"]+bl[,"Age"]+Cohort@time.to.state[,49]), col="blue",xlab="Year of Death",xlim = c(2017,2030))
######## censoring: ##################################
# 
#  aa =  floor(bl[,"BirthYear"]+bl[,"Age"]+cohort@time.to.state[,49])
#  aa[aa <= 2000] = NA
#  Death_since_2011 = aa
#  hist(Death_since_2011, col="green",xlab = "Year of Death", xlim = c(2017,2030))

#######################################################
Patient = Cohort
PatientAll = Patient
PatientAll[is.na(Patient)] <- 0
 # Kernel Density Plot
 d <- density( (bl[,"BirthYear"]+bl[,"Age"]+PatientAll[,49] + PatientAll[,50] + PatientAll[,51] + PatientAll[,52] ))
 plot(d, main = "Density for Year of Death", xlab = "Year of Death")
# plot(floor(bl[,"BirthYear"]+bl[,"Age"]+Patient[,49]), col="blue",xlab="Year of Death")
##### Average time to each state ################################
 Average_Time = 12 * apply(PatientAll,2, mean)
 print("Person time at each state: (Average time(months) at each stage)")
 print(Average_Time)
 write.table(Average_Time, "Average_Time") # save the average time to each state
 ###################################################################
 PatientSum = apply(cohort@time.to.state,2,sum)
 c = apply (who, 2, sum)
 AverageTimeToState = PatientSum / c
 ###################################################################
 for (i in 1: cohortSize)
 if (is.na(Patient[i,49])) Patient[i,49] = maxTime
 PatientAgeDie = Patient[,49] + bl [, "Age"]
 PatientYearDie = Patient[,49] + bl [, "Age"] + bl[,"BirthYear"]
 hist(PatientAgeDie, breaks=10, col="red",xlab = "Age of Death", xlim = c(0,150),lwd=2)
 hist(PatientYearDie, breaks=10, col="red",xlab = "Year of Death",xlim = c(1995,2014),lwd=2)
 xx = c(1995:2014)
 plot(PatientYearDie,col = "red")
 
 plot(xx,y/sd(y),col = "red",type= "l", ylim = c(0,9),lwd=2)
 grid(NA, 5, lwd = 2)
 # Kernel Density Plot
 #d <- density( PatientAgeDie)
 #plot(d)
  h<-hist(PatientAgeDie, breaks=10, col="red", xlab="Age of Death", xlim = c(0,150),main="Histogram with Normal Curve") 
 xfit<-seq(min(PatientAgeDie),max(PatientAgeDie),length=40)
 yfit<-dnorm(xfit,mean=mean(PatientAgeDie),sd=sd(PatientAgeDie))
 yfit <- yfit*diff(h$mids[1:2])*length(PatientAgeDie)
 lines(xfit, yfit, col="blue", lwd=2)
 
# ###################################################
 #PatientAgeDC = Patient[,31:35] + bl [, "Age"]
# hist(PatientAgeDC, breaks=10, col="blue",xlab = "Age of HCC", xlim = c(0,100))
# h<-hist(PatientAgeDC, breaks=10, col="blue", xlab="Age of DC",main="Histogram with Normal Curve") 
# xfit<-seq(min(PatientAgeDC),max(PatientAgeDC),length=40)
# yfit<-dnorm(xfit,mean=mean(PatientAgeDC),sd=sd(PatientAgeDCC))
# yfit <- yfit*diff(h$mids[1:2])*length(PatientAgeDC)
# lines(xfit, yfit, col="blue", lwd=2)
# ####################################################
#PatientAgeHCC = Patient[,36:42] + bl [, "Age"]
 message("Make sure to have enough data to graph Age at HCC!")
 message("==================================================================================")
 hist(PatientAgeHCC, breaks=10, col="black",xlab = "Age of HCC", xlim = c(0,100))
 h<-hist(PatientAgeHCC, breaks=10, col="black", xlab="Age of HCC",main="Histogram with Normal Curve") 
 # xfit<-seq(min(PatientAgeHCC),max(PatientAgeHCC),length=40)
# yfit<-dnorm(xfit,mean=mean(PatientAgeDC),sd=sd(PatientAgeDC))
# yfit <- yfit*diff(h$mids[1:2])*length(PatientAgeDC)
# lines(xfit, yfit, col="blue", lwd=2)
# #####################################################
setwd(wd)
#########################################################
end.time<-Sys.time()
########################################################
 cinc2 <-  cumulativeIncidence(cohort,times = seq(0,maxTime,1))   
 plot(cinc2[,49],main = "Cumulative incidence of death state", ci = TRUE)
 post2  <- transitionProbabilities(cohort,times = seq(0,100,1))
 plot(post2[,49],main = "Probability of transition to Death", ci = TRUE)
 cincF0 = cinc2[,1] + cinc2[,6] + cinc2[,11] +cinc2[,16]+ cinc2[,21]+cinc2[26]
 cincF1 = cinc2[,2] + cinc2[,7] + cinc2[,12] +cinc2[,17]+ cinc2[,22]+cinc2[27]
 cincF2 = cinc2[,3] + cinc2[,8] + cinc2[,13] +cinc2[,18]+ cinc2[,23]+cinc2[28]
 cincF3 = cinc2[,4] + cinc2[,9] + cinc2[,14] +cinc2[,19]+ cinc2[,24]+cinc2[29]
 cincF4 = cinc2[,5] + cinc2[,10] + cinc2[,15] +cinc2[,20]+ cinc2[,25]+cinc2[30]
 cincDC = cinc2[,31] + cinc2[,32] + cinc2[,33] +cinc2[,34]+ cinc2[,35]+cinc2[36]
 cincHCC = cinc2[,37] + cinc2[,38] + cinc2[,39] +cinc2[,40]+ cinc2[,41]+cinc2[42]
 cincLT = cinc2[,43] + cinc2[,44] + cinc2[,45] +cinc2[,46]+ cinc2[,47]+cinc2[48]