library(readstata13)
library(ggplot2)

masterdata.popsize = NULL
for(scen in c("mortval_cohort_nomortadj_mig2004","mortval_cohort_minimal_nomortadj_mig2004","mortval_cohort_mortadj2004_mig2004","mortval_cohort_minimal_mortadj2004_mig2004")) {
    fem.census.popsize = read.dta13(paste0("./data/",scen,"_census_popsize.dta"))
    
    fem.census.popsize$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(fem.census.popsize))
    fem.census.popsize$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(fem.census.popsize))
#   fem.census.popsize$series = "FEM vs. Census forecast"

#		popsize.census = read.dta13("./data/census_proj2004_vs_est.dta")
#		popsize.census$excess_pop = popsize.census$census_diff / 100000.0
#		popsize.census$series = "Census estimate vs. Census forecast"
#		popsize.census$Mortality.Adjustment = rep(ifelse(grepl("nomortadj",scen),"Mortality Adjustment OFF",ifelse(grepl("mortadj2016",scen),"Mortality Adjustment 2016","Mortality Adjustment 2004")), nrow(popsize.census))
#    popsize.census$Model = rep(ifelse(grepl("minimal",scen),"Minimal Model","Full Model"), nrow(popsize.census))
#    popsize.census$pop_ratio = rep(NA, nrow(popsize.census))
#		popsize.census = popsize.census[,c("year","sex","excess_pop","pop_ratio","series","Mortality.Adjustment","Model")]
		

    masterdata.popsize = rbind(masterdata.popsize, fem.census.popsize) # , popsize.census)   
}


for(sex in c("M","F")) {

sexlab = ifelse(sex=="M","male","female")
  
#  for(model in c("Minimal Model","Full Model")) {

   # aes(color=factor(projectionyear)) + labs(color="Census Projections Year") + theme(legend.position = "top") + 

  	    p = ggplot(masterdata.popsize[masterdata.popsize$sex==sex,], aes(x=year,y=excess_pop)) + geom_point(aes(color=factor(series))) + 
  	    	labs(color="Projection model") +
                        geom_hline(yintercept=0.0) + 
                        facet_wrap(Model ~ Mortality.Adjustment)
        p = p + xlab("Year") + ylab("Percent difference relative to initial population size") + labs(title=paste0("Excess size of FEM population relative to Census for ",sexlab,"s"), subtitle="Difference in population size between 2004 FEM stock population, Census estimates, and 2004 Census projections as a percent of the initial population size")
        # print(p)
        ggsave(filename=paste0("./plots/excess_popsize_",sexlab,".pdf"),plot=p, device="pdf", height=8.5, width=11)


  	p = ggplot(masterdata.popsize[masterdata.popsize$sex==sex,], aes(x=year,y=pop_ratio)) + geom_point() + 
                        geom_hline(yintercept=1.0) + 
                        facet_grid(Model ~ Mortality.Adjustment)
        p = p + xlab("Year") + ylab("Population size ratio (FEM/Census)") + labs(title=paste0("Excess size of FEM 2004 stock population ",sexlab,"s"), subtitle="Ratio of FEM simulated population size and 2004 Census projection")
        # print(p)
        ggsave(filename=paste0("./plots/popsize_ratio_",sexlab,".pdf"),plot=p, device="pdf", height=8.5, width=11)


#  }
}