library(igraph)
setwd("~/Documents/Violence Cascades/")
load('Raw Data/lcc.RData')
library(foreach)
library(doMC)
registerDoMC(cores=4)

# gather data
vic_days = as.numeric(unlist(lcc_verts[,16:21]))
idxs = which(!is.na(vic_days)) %% vcount(lcc)
vic_days = vic_days[!is.na(vic_days)]
n.infected = sum(lcc_verts$vic)
birthyears = as.numeric(format(as.Date(lcc_verts$dob[idxs]),'%Y'))

# find all pairs of vic neighbors to check
## nbrs = neighborhood(lcc,order=1,nodes=vic_ids); save(nbrs,file='Raw Data/vic-nbrs.RData')
load('Raw Data/vic-nbrs.RData')
vic.nbrs = matrix(NA,9568,2)
idx = 1
for (i in 1:n.infected){
  u = vic_ids[i]
  nbhd = as.numeric(nbrs[[i]])
  nbhd = intersect(vic_ids,nbhd)
  nbhd = nbhd[u<nbhd] # only looks at pairs with u<v to avoid double-counting
  if (length(nbhd)==0) next
  vic.nbrs[idx:(idx+length(nbhd)-1),] = cbind(rep(u,length(nbhd)),nbhd)
  idx = idx + length(nbhd)
}

# run simulation
n = 10000

ptm = proc.time()
simdata = foreach(q=1:n, .combine=rbind) %dopar% {  
  # shuffle infection dates
  sim.dates = rep(0,length(vic_days))
  for(y in unique(birthyears)){
    year_ids = which(birthyears==y)
    if (length(year_ids)>1) sim.dates[year_ids] = sample(vic_days[year_ids])
  }
#   sim.dates = vic_days # data
  
  vic.dt = rep(0,dim(vic.nbrs)[1])
  for(i in 1:dim(vic.nbrs)[1]){
    mindt = 5000
    tus = sim.dates[which(vic.nbrs[i,1]==idxs)]
    tvs = sim.dates[which(vic.nbrs[i,2]==idxs)]
    for(tu in tus){
      for(tv in tvs){
        dt = abs(tu-tv)
        if(dt<mindt) mindt = dt
      }
    }
    vic.dt[i] = mindt
  }
  
#   mean.time[q] = mean(vic.dt)
#   med.time[q] = median(vic.dt)
#   mean.30[q] = mean(vic.dt<=30)
#   mean.60[q] = mean(vic.dt<=60)
  return(c(mean(vic.dt),median(vic.dt),mean(vic.dt<=30),mean(vic.dt<=100)))
}
print(proc.time()-ptm)
simdata = as.data.frame(simdata,row.names=1:dim(simdata)[1])
colnames(simdata)=c('mean.time','med.time','mean.30','mean.100')
# save(simdata,file='Results/sims-Sept21.RData')

####### plot simulation data #######
# mean=662.6, med=488, m30=0.076, m60=0.119, m100=0.167
d = c(662.6,488,0.076,0.167)
xlabs = c('Mean time between infections',
          'Median time between infections',
          'Infections within 30 days',
          'Infections within 100 days')
par(mfrow=c(2,2))
for(i in 1:4){
  data = d[i] 
  sdata = simdata[,i]
  xl = c(0.99*min(min(sdata),data),1.01*max(max(sdata),data))
  h = hist(sdata,50,xlim=xl,col='#1f78b4',freq=T,border=NA,axes=F,
       xlab=xlabs[i],ylab='Relative Frequency',main=NULL)
  axis(1,at=pretty(xl,5))
  axis(2,at=pretty(h$counts,3),lab=pretty(h$counts/n,3))
  abline(v=data,lwd=4,col='#e41a1c')
  box(lwd=1.1)
}

##### inf dates histogram
vic_dates = as.Date(unlist(lcc_verts[,10:15]))
vic_dates = vic_dates[!is.na(vic_dates)]
hist(vic_dates, breaks='months',freq=T,col='lightblue',format='%Y',
     xlab='Month',ylab='Number of Shootings',main='')

### data vs sim hist comparison
hist(vic.time.data,breaks=seq(0,3000,by=10),col=rgb(0,0,1,1/2),freq=F,
     xlim=c(0,max(vic.time,vic.time.data)),
     xlab='Days between infections',main='')
hist(vic.time,breaks=seq(0,3000,by=10),col=rgb(1,0,0,1/2),add=T,freq=F)
legend("topright", inset=0.05, 
       c("Data", "Simulation"), 
       fill=c(rgb(0,0,1,1/2),rgb(1,0,0,1/2)))
box(which='plot')