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

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

padron2015_mun = readRDS(paste0(path_to_external_data, "padron2015_mun.Rds"))

##########################
### Province Codes ##
##########################

prov_codes = read.csv('aps_province_codes.csv')

###################
## APS Commuting ##
###################

aps_file_names = c('EPAT2014Q01','EPAT2014Q02','EPAT2014Q03', 'EPAT2014Q04', 'EPAT2015Q01','EPAT2015Q02','EPAT2015Q03','EPAT2015Q04') 
hw_list = mclapply(aps_file_names, function(fn)return(as.data.table(read.fwf(paste0(path_to_external_data, fn, '.txt'), col.names = c('home', 'work', 'Nn'), widths=c(-5,2,-76,2, -84, 7), colClasses='character'))), mc.cores=ncores)

hw = rbindlist(hw_list)
hw$Nn = (as.numeric(substr(hw$Nn, 1, 5)) + as.numeric(substr(hw$Nn, 6, 7))/100)/length(aps_file_names)

hw$work[hw$work == "  "] = NA
hw$home[which(hw$home == "  ")] = NA

nrow(hw[!is.na(work & !is.na(home)])])



setkey(hw, home, work)
ods  = hw[, sum(Nn), by=list(home, work)]
names(ods) = c('o', 'd', 'n')
ods$o = as.integer(ods$o)
ods$d = as.integer(ods$d)

APS_ods = merge(ods, prov_codes, by.x='o', by.y='code', all.x=T, all.y=F)
names(APS_ods)[4] = 'home_prov'
APS_ods = merge(APS_ods, prov_codes, by.x='d', by.y='code', all.x=T, all.y=F)
names(APS_ods)[5] = 'work_prov'

APS_ods = APS_ods[!is.na(o) & !is.na(d)]

# saving full commuting table
saveRDS(APS_ods, file="APS_ods.Rds")

# adding work-home part of commute
temp_aps_ods = APS_ods[, c('o', 'd', 'n'), with=F]
names(temp_aps_ods) = c( "d", "o", "n_wh")
APS_ods_dos = merge(APS_ods, temp_aps_ods, by=c("d", "o"))

APS_ods_dos$n2 = APS_ods_dos$n + APS_ods_dos$n_wh

muni_sept_probs = readRDS("muni_sept_probs.Rds")

muni_sept_probs = merge(muni_sept_probs, padron2015_mun, by='PROVMUN', all=T)
muni_sept_probs$prov_code = as.integer(substr(muni_sept_probs$PROVMUN, 1,2))
muni_sept_probs = muni_sept_probs[!is.na(PROVMUN) & !is.na(Total) & !prov_code%in%c(35, 38,51,52) ]
muni_sept_probs[, prov_pop:=sum(Total), by=prov_code]
muni_sept_probs$pop_frac = muni_sept_probs$Total/muni_sept_probs$prov_pop
muni_sept_probs[, tiga_prob_prov_mean:=mean(tiga_prob, na.rm=T), by=prov_code]
muni_sept_probs$tiga_prob_imp = muni_sept_probs$tiga_prob
muni_sept_probs$tiga_prob_imp[which(is.na(muni_sept_probs$tiga_prob))] = muni_sept_probs$tiga_prob_prov_mean[which(is.na(muni_sept_probs$tiga_prob))]
muni_sept_probs$tiga_prob_frac = muni_sept_probs$tiga_prob_imp * muni_sept_probs$pop_frac
prov_sept_probs = muni_sept_probs[, list(tiga_prob=sum(tiga_prob_frac)), by=prov_code]

APS_ods_dos = APS_ods_dos[!o%in%c(35, 38,51,52) & !d%in%c(35, 38,51,52)]
APS_ods_dos = merge(APS_ods_dos, prov_sept_probs, by.x='o', by.y='prov_code', all=T)
APS_ods_dos$transport_risk_od = APS_ods_dos$n2 * APS_ods_dos$tiga_prob

saveRDS(APS_ods_dos, file='APS_ods_dos.Rds')

APS_ods_dos_interprov = APS_ods_dos[o != d & !is.na(d) & !is.na(o), list(home_prov, work_prov, o,d,transport_risk_od)]

saveRDS(APS_ods_dos_interprov, file='APS_ods_dos_interprov.Rds')


sources = APS_ods_dos_interprov[, list(source_n=sum(transport_risk_od)), by=list(home_prov, o)]

sinks = APS_ods_dos_interprov[, list(sink_n=sum(transport_risk_od)), by=list(work_prov,d)]

source_sink = merge(sources, sinks, by.x='o', by.y='d')
source_sink$net = source_sink$source_n - source_sink$sink_n

saveRDS(source_sink, file='source_sink.Rds')


sst = source_sink[, list(home_prov, source_n, sink_n, net)]
names(sst) = c("Province", "Outflow (N)", "Inflow (N)", "Net Outflow (N)")
saveRDS(sst, file='APS_source_sink_table.Rds')

ssd = source_sink[, list(o, home_prov, source_n, sink_n, net)]
saveRDS(ssd, file='APS_source_sink_data.Rds')

bcn_out = APS_ods_dos_interprov[home_prov=="Barcelona", list(o,d, home_prov, work_prov, transport_risk_od)]

saveRDS(bcn_out, file='APS_bcn_outflows.Rds')

