### SHOULD BE RUN ON CLUSTER ###

rm(list=ls())
library(foreign)
library(rstanarm)
library(rstan)
library(Hmisc)
library(data.table)
library(parallel)

options("stringsAsFactors" = FALSE) 

setwd('~/research/gTiger/Article_citsci/open_code')
source('settings/parameters.r')
source('settings/paths.r')
source('settings/functions.r')

setwd(path_to_data)

ncores = 10

rstan_options(auto_write = TRUE)
options(mc.cores = ncores)

Apr1Julian = as.POSIXlt("01-04-2015", format = "%d-%m-%Y")$yday
Dec1Julian = as.POSIXlt("01-12-2015", format = "%d-%m-%Y")$yday

Sep1Julian = as.POSIXlt("01-09-2015", format = "%d-%m-%Y")$yday
Sep30Julian = as.POSIXlt("30-09-2015", format = "%d-%m-%Y")$yday

std2 = function(x) (x-mean(x,na.rm=T))/(2*sd(x, na.rm=T))
invstd2 = function(x, d) (x*2*sd(d, na.rm=T))+mean(d,na.rm=T)
p2o = function(x) x/(1-x)


std2_d = function(x, d) (x-mean(d, na.rm=T))/(2*sd(d, na.rm=T))

d = readRDS('mosquito_alert_model_dataset.Rds')

d$trials = 1

# temporarily taking out main model  and M0a since I've already run these and no need to do theme again
# list(name="mosquito_alert_model", vars=c("any_reliable_reports", "trials", "SE_adj_log_std2", "sea_days_std2",  "sea_days_std2_sq", "sea_days_std2_cu", "TigacellID")),  
# list(name="mosquito_alert_model_M0a", vars=c("any_reliable_reports", "trials", "sea_days_std2",  "sea_days_std2_sq", "sea_days_std2_cu", "TigacellID")),
# list(name="mosquito_alert_model_M0b", vars=c("any_reliable_reports", "trials", "SE_adj_log_std2", "sea_days_std2",  "sea_days_std2_sq", "TigacellID")), 

these_models = list(list(name="mosquito_alert_model_M0c", vars=c("any_reliable_reports", "trials", "SE_adj_log_std2", "sea_days_std2", "TigacellID")), list(name="mosquito_alert_model_M0d", vars=c("any_reliable_reports", "trials", "SE_adj_log_std2", "sea_days_std2",  "sea_days_std2_sq", "TigacellID")), list(name="mosquito_alert_model_M0e", vars=c("any_reliable_reports", "trials", "SE_adj_log_std2", "sea_days_std2",  "sea_days_std2_sq", "sea_days_std2_cu", "TigacellID")))


for(this_model in these_models){

M = readRDS(paste0('untracked/', this_model$name, '.Rds'))

modeled_cells = unique(d$TigacellID)
day_range = 1:365
dp = data.table(TigacellID = rep(modeled_cells, length(day_range)), sea_days = unlist(lapply(day_range, function(d) rep(d, length(modeled_cells)))))

dp$sea_days_std2 = std2_d(dp$sea_days, d$sea_days)
dp$sea_days_std2_sq = dp$sea_days_std2^2
dp$sea_days_std2_cu = dp$sea_days_std2^3

SE_1_std2 = std2_d(log(1), d$SE_adj_log)

dp$SE_adj_log_std2 = SE_1_std2

ntrials = 1000
sample_size = 1000

dp$trials = ntrials

dp$any_reliable_reports = 0

mf = readRDS('tigacell_muni_fracs_ine.Rds')

### daily predictions 1 April - 1 December
selected_sea_days = Apr1Julian:Dec1Julian

M_preds_dt_list = mclapply(selected_sea_days, function(this_sea_days){
	print(this_sea_days)
	flush.console()
	this_dp = dp [sea_days==this_sea_days, this_model$vars, with=FALSE]
	M_preds = posterior_predict(M, this_dp, draws=sample_size)
	colnames(M_preds) = this_dp$TigacellID
	M_preds_dt = as.data.frame(t(M_preds))
	setDT(M_preds_dt, keep.rownames=TRUE)
	M_preds_dt$sea_days = this_sea_days
	return(M_preds_dt)
}, mc.cores=ncores)

M_preds_dt = rbindlist(M_preds_dt_list)

saveRDS(M_preds_dt, file=paste0('untracked/M_preds_dt_', this_model$name, '.Rds'))

rm(M_preds_dt)
rm(M_preds_dt_list)
rm(M)
rm(dp)
}

