#######################
### SETUP WORKSPACE ###
#######################

library(maptools)
library(raster)
library(spatstat)
library(rgdal)
library(rgeos)

## Clear the workspace
rm(list=ls())

##############################################################
### NOTE: This line should be changed to point to the folder
###       where you have unzipped the contents of the
###       replication archive.
##############################################################

## Set working directory
setwd("C:/warren_jpr2015_rep")

## Import data from load_data.R 
load('load_data.RData')

## Import data from krig_data.R 
load('krig_data.RData')

###################################
### SET GLOBAL PARAMETER VALUES ###
###################################

#Set radii for obvservation windows 
#(maximum distance from survey location)
radius_list <- c(30000, 40000, 50000) #meters

#Set polygon size lower cut-off for obsevation windows
#(fraction of full circle with given radius)
small_size <- 0.01


##################################
### Create Observation Windows ###
##################################

#Convert survey locations to SpatialPoints
dhs.pnts <- as(dhs, "SpatialPoints")

#Loop over distances in radius list
for (radius in radius_list) {
  
  cat("Calculating buffer radius:", radius, "\n")  
  
  #Buffer points by radius
  buffer.list <- vector(mode="list", length=length(dhs.pnts))
  for (pnt in 1:length(dhs.pnts)) {
    buf <- gBuffer(dhs.pnts[pnt], width=radius)
    buf@polygons[[1]]@ID <- as.character(pnt)
    buffer.list[pnt] <- buf@polygons
  }
  rm(buf, pnt)
  
  #Merge buffered points from same surveys into single polygons
  buffer.clist <- vector(mode="list", length=(nrow(dhs.cowcode.year)-1))
  for (surv in c(1:nrow(dhs.cowcode.year))) {
    pnts <- buffer.list[ (dhs$cowcode == dhs.cowcode.year$cowcode[surv] &
                          dhs$dhsyear == dhs.cowcode.year$dhsyear[surv] ) ]
    sp <- SpatialPolygons(pnts)
    poly <- gUnaryUnion(sp)@polygons
    poly[[1]]@ID <- as.character(surv)    
    if (dhs.cowcode.year$dhsyear[surv] <= 2010) {
      if (dhs.cowcode.year$cowcode[surv] == 500) {
        #Pull out 1st Uganda survey
        uganda.poly1 <- poly
      }
      #All but 2nd Uganda survey recorded in full list
      buffer.clist[surv] <- poly
    } else {
        #Pull out 2nd Uganda survey
        uganda.poly2 <- poly
      }
  }
  rm(pnts, sp, poly, surv)
  
  #Convert buffered locations to SpatialPolygons
  dhs.sp <- SpatialPolygons(buffer.clist)
  proj4string(dhs.sp) <- proj.string
  
  #Restrict DV observation window to surveyed country shapes
  dhs.sp <- gIntersection(dhs.sp, dhs.cshapes, byid=TRUE)
  
  for (s in c('all', 'ug1', 'ug2')) {
     
    if (s %in% c('ug1', 'ug2')) {
      if (s == 'ug1') {
        #1st Uganda survey
        dhs.sp <- SpatialPolygons(uganda.poly1)
      } else {
        #2nd Uganda survey
        dhs.sp <- SpatialPolygons(uganda.poly2)       
      }      
      proj4string(dhs.sp) <- proj.string
    }
    
    #Combine polygons
    dhs.sp <- gUnaryUnion(dhs.sp)
    
    #Exclude very small polygons
    poly.list <- lapply( dhs.sp@polygons , slot , "Polygons" )
    row <- 0
    threshold <- pi * radius^2 * small_size
    small <- vector()
    for (i in c(1:length(poly.list))) {
      polys <- poly.list[[i]]
      for (j in c(1:length(polys))) {
        row <- row + 1
        area <- slot(polys[[j]], "area")
        hole <- slot(polys[[j]], "hole")
        if (area < threshold) {
          small <- c(small, row)
        }
      }
    }
    poly.list[[1]] <- poly.list[[1]][-small]
    dhs.sp <- Polygons(poly.list[[1]], 1)
    dhs.sp <- SpatialPolygons(list(dhs.sp))
    proj4string(dhs.sp) <- proj.string
    
    #Rasterize observation window
    dhs.owin <- rasterize(dhs.sp, r)
    dhs.owin <- as(dhs.owin, "SpatialGridDataFrame")
    dhs.owin <- as.im.SpatialGridDataFrame(dhs.owin)
    dhs.owin <- as.owin(dhs.owin)
    
    #Save results
    radi <- radius / 1000
    name.sp <- paste0("dhs.", s, ".sp.r", radi)
    assign(name.sp, dhs.sp)
    name.im <- paste0("dhs.", s, ".owin.r", radi)
    assign(name.im, dhs.owin)
    
    rm(dhs.owin, dhs.sp, poly.list, polys, 
       row, threshold, area, hole, i, j, small)
    
  } # end for s
  
  #Create 30km observation window for Uganda by 
  #finding the intersection between the polygons generated 
  #for the first survey year and the second survey year
  if (radi == 30) {
    name.out <- paste0("dhs.uganda.owin")
    name.owin1 <- paste0('dhs.ug1.owin.r', radi)
    name.owin2 <- paste0('dhs.ug2.owin.r', radi)
    
    owin <- intersect.owin(get(name.owin1),
                           get(name.owin2))
    
    #Restrict DV observation window to Uganda
    uganda.sp <- as(uganda.cshape, "SpatialPolygons")
    owin <- owin[as.owin(uganda.sp)]
    
    assign(name.out, owin)
  }
} # end for radius

rm(buffer.list, buffer.clist, name.im, name.sp,
   dhs.pnts, uganda.cshape, uganda.poly2, radius, radi, s)


#########################################
### Calculate Media Shifts for Uganda ###
#########################################

rad_diff.im <- eval.im(rad.ug2.im - rad.im)
cell_diff.im <- eval.im(cell.ug2.im - cell.im)

rm(rad.ug2.im, cell.ug2.im)

#######################################
### Create Dependent Variable (GED) ###
#######################################

#Convert GED data to SpatialPointsDataFrame
coordinates(ged) <- ~ lon + lat 
proj4string(ged) <- CRS("+init=epsg:4326")

#Re-project (using equal area projection)
ged <- spTransform(ged, CRS(proj.string))

#Keep only deadly events
ged <- ged[ged$best_est >= 1, ]

#Keep only geo-located events
ged <- ged[!is.na(ged$where_prec),]
ged <- ged[ged$where_prec <= 3, ]

#Jitter locations, according to 
#level of geographic precision
#(see GED codebook)
ged1 <- ged[ged$where_prec == 1, ]
ged2 <- ged[ged$where_prec == 2, ]
ged3 <- ged[ged$where_prec == 3, ]

ged1 <- as.ppp(ged1)
ged2 <- as.ppp(ged2)
ged3 <- as.ppp(ged3)

ged1 <- rjitter(ged1, 5000) #meters
ged2 <- rjitter(ged2, 12500)
ged3 <- rjitter(ged3, 25000)

ged <- superimpose(ged1, ged2, ged3)
ged <- as(ged, "SpatialPointsDataFrame")
names(ged) <- gsub("marks.", "", names(ged))
rm(ged1, ged2, ged3)

#Merge GED Events and Survey Country-Years
ged$GWCODE <- ged$gwno
ged.cdata <- as.data.frame(dhs.cshapes[, c("GWCODE", "dhsyear", "cowcode")])
names(ged.cdata)[1] <- 'gwno'
ged.cdata <- merge(ged, ged.cdata, by="gwno", all=FALSE)

#Exclude points over missing cells
missing.vars <- c("rad", "elec", "elev", "popdens", "urban", "gdp_ppp")
for (var in missing.vars) {
  varname <- paste0(var, ".im")
  rast.df <- as(get(varname), "SpatialGridDataFrame")
  ged.cdata <- ged.cdata[which(is.na(over(ged.cdata, rast.df)) == FALSE), ]
}
rm(var, varname, rast.df)


##Subset GED by year, to create
##two versions of dependent variable:

#1: Include observations following date of survey collection
ged.postsurv    <- ged.cdata[ged.cdata$year >= ged.cdata$dhsyear, ]

#2: Include obsevations occuring between 2005 and 2010
ged.post2005    <- ged.cdata[ged.cdata$year >= 2005, ]

ged.dv.sets <- c(
  'postsurv',
  'post2005')

#Convert event data to PPP format
#Add observation windows for each radius distance
for (set in ged.dv.sets) {
  varname <- paste0('ged.', set)
  cat(paste0('Creating [', varname, '.ppp] \n'))
  ppp <- as(get(varname), "SpatialPointsDataFrame")
  ppp <- as.ppp(ppp)
  for (radius in radius_list) {
    radi <- radius / 1000
    varname.r <- paste0(varname, ".r", radi)
    owin.name <- paste0("dhs.all.owin.r", radi)
    owin <- get(owin.name)
    assign(varname.r, ppp[owin])
  }
  rm(ppp, owin, list=varname)
}
rm(set, varname, varname.r, owin.name, radius, radi)

##########################
### Generate Lagged DV ###
##########################

##Two versions of "prior violence", corresponding to DV versions:

#1: Include obsevations for the 5 years preceding each survey
ged.pre5surv   <- ged.cdata[ged.cdata$year <  ged.cdata$dhsyear &
                            ged.cdata$year >= (ged.cdata$dhsyear - 5), ]

#2: Include observations for the 5 years between 2000 and 2004
ged.pre2005 <- ged.cdata[ged.cdata$year >= 2000 &
                         ged.cdata$year <= 2004, ]

ged.iv.sets <- c(
  "pre5surv",
  'pre2005')

#Create IV observation window for surveyed countries
iv.owin <- rasterize(dhs.cshapes, r)
iv.owin <- as(iv.owin, "SpatialGridDataFrame")
iv.owin <- as.im.SpatialGridDataFrame(iv.owin)
iv.owin <- as.owin(iv.owin)

# Calculate point densities
for (set in ged.iv.sets) {
  varname <- paste0('ged.', set)
  varname.im <- paste0(varname, '.im')
  cat(paste0('Creating [', varname.im, '] \n'))
  ppp <- as(get(varname), "SpatialPointsDataFrame")
  ppp <- as.ppp(ppp)
  ppp <- ppp[iv.owin]
  unmark(ppp)
  dens <- density(
    ppp, 
    leaveoneout=FALSE, 
    sigma=500000,
    eps=20000,
    diggle=TRUE)
  dens[eval.im(dens < 0)] <- 0
  assign(varname.im, dens)
  rm(ppp, dens, list=varname)
}
rm(ged, ged.cdata, set, varname, varname.im, r)


#################################
### Create Quadrature Schemes ###
#################################

for (radius in c(radius_list, 999)) {
  
  n <- 10000
  radi <- radius / 1000
  sp.name <- paste0("dhs.all.sp.r", radi)
  quad.name <- paste0("quad.r", radi)   

  #Flag to create Uganda sample
  if (radius == 999) {
    n <- 300
    sp.name <- paste0("dhs.uganda.owin")
    quad.name <- paste0("quad.ug")
  }
  
  cat(paste0('Creating [', quad.name, '] \n'))
  
  #Select observation window for given radius
  dhs.sp <- as(get(sp.name), "SpatialPolygons")
  proj4string(dhs.sp) <- proj.string 
  
  #Randomly sample locations within observation window
  samp <- spsample(dhs.sp, n=n, type="random")
  
  #Exclude points in missing cells
  missing.vars <- c("rad", "elec", "elev", "popdens", 
                    "urban", "gdp_ppp", "ged.pre2005")
  for (var in missing.vars) {
    varname <- paste0(var, ".im")
    rast.df <- as(get(varname), "SpatialGridDataFrame")
    proj4string(rast.df) <- proj.string 
    samp <- samp[is.na(over(samp, rast.df)) == FALSE]
  }
  
  samp <- as.ppp(samp)
  assign(quad.name, samp)
  
} # end for radius

rm(radi, radius, sp.name, dhs.sp, var, varname, rast.df, samp, quad.name)

####################
### Save Results ###
####################
radi_list <- radius_list / 1000
varlist.out <- c(paste0('quad.r', radi_list), 'quad.ug',
                 paste0('dhs.all.owin.r', radi_list), 'dhs.uganda.owin',
                 paste0('ged.', ged.iv.sets, '.im'),
                 paste0(rep('ged.', 2), ged.dv.sets, '.r', rep(radi_list, 2)),
                 'rad.im', 'rad_diff.im', 'cell.im', 'cell_diff.im', 
                 'elec.im', 'elev.im', 'gdp_ppp.im', 'popdens.im',
                 'urban.im', 'shrpower.im', 'cowcode.im')

save(list=varlist.out, file="ppm_data_new.RData")

