library(readstata13)
library(ggplot2)


masterdata.survival = NULL
for(scen in c("mortval_cohort_nomortadj","mortval_cohort_minimal_nomortadj","mortval_cohort_mortadj2004","mortval_cohort_minimal_mortadj2004")) {

	 survival = read.dta13(paste0("./data/",scen,"_cdc_bellmiller_survival.dta"))

	 survival$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(survival))
         survival$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(survival))

	 survival$sex = ifelse(survival$male==0, "Female", "Male")

	 masterdata.survival = rbind(masterdata.survival, survival) 

}

# make a separate plot for each 10-year range of birth cohorts
# 1910-1919
# 1920-1929
# etc.

for(decade in 1:5) {
#	for(male in 0:1) {
	   	
#		sexlab = ifelse(male==1,"male","female")

#		surv.plot = masterdata.survival[masterdata.survival$Model=="Full Model" & masterdata.survival$Mortality.Adjustment=="Mortality Adjustment 2004" & masterdata.survival$male==male & masterdata.survival$birthcohort >= 1900+10*decade & masterdata.survival$birthcohort <= 1900+10*decade+9,]
		surv.plot = masterdata.survival[masterdata.survival$Model=="Full Model" & masterdata.survival$Mortality.Adjustment=="Mortality Adjustment 2004" & masterdata.survival$birthcohort >= 1900+10*decade & masterdata.survival$birthcohort <= 1900+10*decade+9,]

		print(table(surv.plot$birthcohort))


		p = ggplot(surv.plot, aes(x=year,y=cumsurvprob)) + 
                geom_point(aes(color=factor(series), shape=factor(series)), alpha=0.8, size=0.75) + 
                labs(color="Source") + labs(shape="Source") +
                facet_wrap(~ sex + birthcohort, ncol=ifelse(decade==5,4,5)) +
		coord_cartesian(ylim = c(0.03, 1.0)) + 
		scale_x_continuous(breaks=seq(2004,2016,4)) + 
		scale_shape_manual(values=c(16, 17, 18)) +
		scale_color_manual(values=c("#E69F00","#009E73","#CC79A7")) +
		theme(axis.text=element_text(size=6.5))		
#        	p = p + xlab("Year") + ylab("Survival probability") + labs(title=paste0("Forecasted vs. actual survival probability for ",sexlab,"s by birth year"))
        	p = p + xlab("Year") + ylab("Survival probability") + labs(title=paste0("Forecasted vs. actual survival probability by birth year"))
        # print(p)
#        ggsave(filename=paste0("./plots/survival_",sexlab,"_byr",1900+10*decade,'-',1900+10*decade+9,"_full.pdf"),plot=p, device="pdf", height=ifelse(decade==5,4.5,8.5), width=11)
        ggsave(filename=paste0("./plots/survival_byr",1900+10*decade,'-',1900+10*decade+9,"_full.pdf"),plot=p, device="pdf", height=ifelse(decade==5,4.5,8.5), width=11)


#	}
}