From 9350ee2c6359562a23cf8efdefdd7de80b2a682e Mon Sep 17 00:00:00 2001
From: Ben Green <bgreen@g.harvard.edu>
Date: Sat, 1 Aug 2015 12:01:19 -0400
Subject: minor updates

---
 R Scripts/predict-victims.R | 98 +++++++++++++++++++++++----------------------
 1 file changed, 51 insertions(+), 47 deletions(-)

(limited to 'R Scripts/predict-victims.R')

diff --git a/R Scripts/predict-victims.R b/R Scripts/predict-victims.R
index 2bda7e2..c662651 100644
--- a/R Scripts/predict-victims.R	
+++ b/R Scripts/predict-victims.R	
@@ -1,8 +1,9 @@
 library(igraph)
+library(data.table)
 library(foreach)
 library(doMC)
 registerDoMC(cores=4)
-setwd('~/Documents/Cascade Project')
+setwd('~/Documents/Violence Cascades/')
 load('Raw Data/lcc.RData')
 load('Results/hyper-lcc.RData')
 load('Results/dag_dat_all.RData')
@@ -10,65 +11,68 @@ source('criminal_cascades/R Scripts/temporal.R')
 source('criminal_cascades/R Scripts/structural.R')
 
 ##### Initialize data
-formula = vic ~ sex + race + age + gang.member #+ gang.name
+formula = vic ~ sex + race + age + gang.member + gang.name
 lcc_verts$sex = as.factor(lcc_verts$sex)
 lcc_verts$race = as.factor(lcc_verts$race)
 lcc_verts$age = as.numeric(lcc_verts$age)
 lcc_verts$gang.name = as.factor(lcc_verts$gang.name)
-# sum(hyp_lcc_verts$vic)/length(days)
-df = data.frame(ir=lcc_verts$ir_no, dem=0, cas=0, comb=0)
+dt = data.table(ir=lcc_verts$ir_no, dem=0, cas=0, comb=0)
 
 alpha = 0.0028
 delta = 0.06
 days = sort(unique(hyp_lcc_verts$vic.day)) # 70:max(hyp_lcc_verts$vic.day, na.rm=T)
-lambdas = c(0,1)#c(0, exp(seq(log(0.0000001), log(.0005), length.out=150)), 1)
-nvics = sum(hyp_lcc_verts$vic.day %in% days)
-edges_all = dag_dat_all
+days = split(days, ceiling(seq_along(days)/456))
+lambdas = c(0, exp(seq(log(0.0000001), log(.95), length.out=50)), 1)
+nvics = sum(hyp_lcc_verts$vic)
+edges_all = data.table(dag_dat_all)
 
 ##### Loop through days
-writeLines(c(""), "Results/log.txt")
-ptm = proc.time()
-correct_rank = foreach (day = days, .combine=rbind) %dopar% {
-  if (which(day==days) %% 100 == 0){sink("Results/log.txt", append=TRUE);cat(paste("day:",day,"\n"))}
-
-  ##### Demographics model
-  vics = match(unique(hyp_lcc_verts$ir_no[which(hyp_lcc_verts$vic.day<day)]),lcc_verts$name)
-  victims = lcc_verts[,c('vic','sex','race','age','gang.member','gang.name')]
-  victims$vic[vics] = TRUE
-  victims$vic[-vics] = FALSE
-  fit = lm(formula, data=victims)
-#   fit = glm(formula, data=victims, family=binomial)
-#   fit = randomForest(formula, data=victims[,1:5], ntree=100)
-  probs = predict(fit, newdata=lcc_verts, type='response')
+correct_rank = c()
+for(i in 1:length(days)){
+  ptm = proc.time()
+  print(c(i,length(days)))
+  ds = unlist(days[i], use.names=F)
+  cr = foreach (day = ds, .combine=rbind) %dopar% {
   
-  ##### Cascade Model
-  edges = edges_all[which(edges_all$t1<day),]
-  f = temporal(edges$t1, day, alpha)
-  h = structural(delta,edges$dist)
-  weights = f*h
-  ids = edges$to
-  irs = hyp_lcc_verts$ir_no[ids]
-  risk = data.frame(id=ids, ir=irs, weight=weights)
-  risk = risk[order(weights, decreasing=T),]
-  risk = risk[match(unique(risk$ir),risk$ir),]
-# maybe need to change this to reflect new algorithm that accounts for \tilde{p}
-
-  ##### Combined Model
-  combined = df#data.frame(ir=attr(probs,'name'), dem=as.numeric(probs), cas=0, comb=0)
-  combined$dem[match(attr(probs,'name'), df$ir)] = as.numeric(probs)
-  combined$cas[match(risk$ir, attr(probs,'name'))] = risk$weight
+    ##### Demographics model
+    vics = match(unique(hyp_lcc_verts$ir_no[which(hyp_lcc_verts$vic.day<day)]),lcc_verts$name)
+    victims = lcc_verts[,c('vic','sex','race','age','gang.member','gang.name')]
+    victims$vic[vics] = TRUE
+    victims$vic[-vics] = FALSE
+#     fit = lm(formula, data=victims)
+    fit = glm(formula, data=victims, family=binomial)
+  #   fit = randomForest(formula, data=victims[,1:5], ntree=100)
+    probs = predict(fit, newdata=lcc_verts, type='response')
+    
+    ##### Cascade Model
+    edges = edges_all[which(edges_all$t1<day),]
+    f = temporal(edges$t1, day, alpha)
+    h = structural(delta,edges$dist)
+    weights = f*h
+    ids = edges$to
+    irs = hyp_lcc_verts$ir_no[ids]
+    
+    risk = data.table(id=ids, ir=irs, weight=weights)
+    if (dim(risk)[1]>0) risk = risk[, list(weight=sum(weight)), by=ir] # max or sum
   
-  ##### Gather results
-  infected_irs = hyp_lcc_verts$ir_no[which(hyp_lcc_verts$vic.day==day)]
-  crday = matrix(nrow=length(infected_irs), ncol=length(lambdas))
-  for (lambda in lambdas){
-    combined$comb = lambda*combined$dem + (1-lambda)*combined$cas
-    c_idx = which(lambdas==lambda)
-    crday[,c_idx] = rank(-combined$comb,ties.method='average')[match(infected_irs,combined$ir)]
+    ##### Combined Model
+    combined = dt
+    combined$dem[match(attr(probs,'name'), dt$ir)] = as.numeric(probs)
+    combined$cas[match(risk$ir, dt$ir)] = risk$weight
+    
+    ##### Gather results
+    infected_irs = hyp_lcc_verts$ir_no[which(hyp_lcc_verts$vic.day==day)]
+    crday = matrix(nrow=length(infected_irs), ncol=length(lambdas))
+    for (lambda in lambdas){
+      combined$comb = lambda*combined$dem + (1-lambda)*combined$cas
+      c_idx = which(lambdas==lambda)
+      crday[,c_idx] = rank(-combined$comb,ties.method='average')[match(infected_irs,combined$ir)]
+    }
+  
+    return(crday)
   }
-
-  return(crday)
+  correct_rank = rbind(correct_rank,cr)
+  print(proc.time()-ptm)
 }
-print(proc.time()-ptm)
 
 # save(correct_rank, file='Results/correct_rank_62815.RData') 
\ No newline at end of file
-- 
cgit v1.2.3-70-g09d2