rm(list=ls())
library(Hmisc)
library(data.table)

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

setwd(path_to_data)

spain_gpdf = readRDS('spain_grid.Rds')

spain_cells = spain_gpdf$TigacellID

season_day_0_2014 = as.POSIXct(strptime('2014-01-01 00:00:00.000000+0000', format="%Y-%m-%d %H:%M:%S%OS%z", tz="gmt"))
trs = readRDS('user_locations_propensities.Rds')

trs = trs[which(trs$TigacellID %in% spain_cells),]
length(which(spain_cells %in% unique(trs$TigacellID)))

biweek_SE = trs[ , list(SE=sum(propensity)), by=list(TigacellID, biweek)]

vrs = readRDS('mosquito_alert_cleaned_reports.Rds')

nrow(vrs)
vrs = vrs[which(vrs$TigacellID %in% spain_cells),]
nrow(vrs)

biweeks = seq(range(vrs$biweek)[1],range(vrs$biweek)[2],1)

ttt = data.table(TigacellID = rep(spain_cells, length(biweeks)), biweek = unlist(lapply(biweeks, function(x) rep(x, length(spain_cells)))))

biweek_vrs = vrs[, list(max_movelab_certainty=max(movelab_certainty_category), reliable_reports=sum(reliable_report)), by=list(TigacellID,biweek)]

ttt = merge(ttt, biweek_vrs, by=c('TigacellID', 'biweek'), all.x=TRUE, all.y=FALSE)
ttt$reliable_reports[which(is.na(ttt$reliable_reports))] = 0
ttt$any_reliable_reports = as.integer(ttt$reliable_reports >=1)
ttt$presence_prob12 = as.integer(!is.na(ttt$max_movelab_certainty) & ttt$max_movelab_certainty >=1)
ttt$presence_prob012 = as.integer(!is.na(ttt$max_movelab_certainty) & ttt$max_movelab_certainty >=0)
ttt = merge(ttt, biweek_SE, by=c('TigacellID', 'biweek'), all.x=TRUE, all.y=FALSE)

ttt$SE[which(is.na(ttt$SE))] = 0

ttt$biweek_start_date = as.Date(ttt$biweek*14, origin=as.Date("2014-01-01"))

ttt$sea_days = as.POSIXlt(ttt$biweek_start_date)$yday

ttt$year = format(ttt$biweek_start_date, '%Y')

ttt$y2015 = as.integer(ttt$year == "2015")

tigacell_total_SE = aggregate(ttt$SE, by=list(ttt$TigacellID), FUN=sum, na.rm=T)
names(tigacell_total_SE) = c('TigacellID', 'total_SE')
unsampled_cells = tigacell_total_SE$TigacellID[which(tigacell_total_SE$total_SE==0)]
sampled_cells = tigacell_total_SE$TigacellID[which(tigacell_total_SE$total_SE>0)]

tigacell_min_SE = aggregate(ttt$SE, by=list(ttt$TigacellID), FUN=min, na.rm=T)
names(tigacell_min_SE) = c('TigacellID', 'min_SE')

tigacell_any_prob12s = aggregate(ttt$presence_prob12, by=list(ttt$TigacellID), FUN=max, na.rm=T)
names(tigacell_any_prob12s) = c('TigacellID', 'any_prob12s')

ttt = merge(ttt, tigacell_min_SE, by="TigacellID", all=T)
ttt = merge(ttt, tigacell_any_prob12s, by="TigacellID", all=T)

ttt$SE_adj = ttt$SE + ttt$min_SE
describe(ttt$SE_adj)

ttt$keep = ttt$SE_adj > 0 | ttt$any_prob12s == 1

describe(ttt$any_reliable_reports)
describe(ttt$presence_prob12)
describe(ttt$presence_prob012)


describe(ttt$any_reliable_reports[which(ttt$keep)])
describe(ttt$presence_prob12[which(ttt$keep)])
describe(ttt$presence_prob012[which(ttt$keep)])

ttt = ttt[which(ttt$keep),]
ttt$SE_adj[which(ttt$SE_adj==0)] = 1e-10
ttt$SE_adj_log = log(ttt$SE_adj)
describe(ttt$SE_adj)
describe(ttt$SE_adj_log)


ttt$SE_adj_log_std2 = std2(ttt$SE_adj_log)
ttt$sea_days_std2 = std2(ttt$sea_days)
ttt$sea_days_std2_sq = ttt$sea_days_std2^2
ttt$sea_days_std2_cu = ttt$sea_days_std2^3

saveRDS(ttt, file="mosquito_alert_model_dataset.Rds")
