From d1e9e8cfabfee86fd7ac9e108fab0d463da24ce8 Mon Sep 17 00:00:00 2001
From: Ben Green <bgreen@g.harvard.edu>
Date: Sat, 29 Aug 2015 16:30:46 -0400
Subject: made simulations much more efficient

---
 R Scripts/sim-analysis.R | 61 ++++++++++++++++++++++++++----------------------
 1 file changed, 33 insertions(+), 28 deletions(-)

(limited to 'R Scripts/sim-analysis.R')

diff --git a/R Scripts/sim-analysis.R b/R Scripts/sim-analysis.R
index 5010485..e368419 100755
--- a/R Scripts/sim-analysis.R	
+++ b/R Scripts/sim-analysis.R	
@@ -12,12 +12,28 @@ birthyears = as.numeric(format(as.Date(hyp_lcc_verts$dob[vic_ids]),'%Y'))
 # nbrs = neighborhood(graph, nodes=vic_ids, order=1)
 load('Results/vic-nbrs.RData')
 
-n = 150
-mean.time = matrix(0,1,n)
-med.time = matrix(0,1,n)
-mean.50 = matrix(0,1,n)
-mean.100 = matrix(0,1,n)
-n.vicpairs = matrix(0,1,n)
+# find all pairs of vic neighbors to check
+vic.nbrs = matrix(NA,14885,2)
+idx = 1
+for (i in 1:n.infections){
+  u = vic_ids[i]
+  nbhd = nbrs[[i]]
+  nbhd = intersect(vic_ids,nbhd)
+  nbhd = nbhd[u<nbhd] # only looks at pairs with u<v to avoid double-counting
+  nbhd = setdiff(nbhd,u)
+  nbhd = nbhd[hyp_lcc_verts$ir_no[u] != hyp_lcc_verts$ir_no[nbhd]] # don't count infections of the same person
+  if (length(nbhd)==0) next
+  vic.nbrs[idx:(idx+length(nbhd)-1),] = cbind(rep(match(u,vic_ids),length(nbhd)),match(nbhd,vic_ids))
+  idx = idx + length(nbhd)
+}
+
+# run simulation
+n = 10000
+mean.time = rep(NA,n)
+med.time = rep(NA,n)
+mean.50 = rep(NA,n)
+mean.100 = rep(NA,n)
+n.vicpairs = rep(NA,n)
 
 hyp_lcc = upgrade_graph(hyp_lcc)
 
@@ -27,28 +43,15 @@ for(q in 1:n){
   graph = hyp_lcc
   vics = vic_ids
   
+  # shuffle infection dates
   sim.dates = rep(0,length(vic_ids))
   for(y in unique(birthyears)){
     year_ids = which(birthyears==y)
     if (length(year_ids)>1) sim.dates[year_ids] = sample(inf.dates[year_ids])
   }
-#   {vics = vic_ids; sim.dates = inf.dates} # data
+#   vics = vic_ids; sim.dates = inf.dates # data
   
-  vic.time = rep(NA,14885)
-  idx = 1
-  for (i in 1:n.infections){
-    u = vics[i]
-    nbhd = nbrs[[i]]
-    nbhd = intersect(vic_ids,nbhd)
-    nbhd = nbhd[u<nbhd] # only looks at pairs with u<v to avoid double-counting
-    nbhd = setdiff(nbhd,u)
-    nbhd = nbhd[hyp_lcc_verts$ir_no[u] != hyp_lcc_verts$ir_no[nbhd]] # don't count infections of the same person
-    if (length(nbhd)==0) next
-    tu = sim.dates[i]
-    tvs = sim.dates[match(nbhd,vic_ids)]
-    vic.time[idx:(idx+length(nbhd)-1)] = abs(tu-tvs)
-    idx = idx + length(nbhd)
-  }
+  vic.time = abs(sim.dates[vic.nbrs[,1]] - sim.dates[vic.nbrs[,2]])
 
   n.vicpairs[q] = length(vic.time)
   mean.time[q] = mean(vic.time)
@@ -59,15 +62,17 @@ for(q in 1:n){
 }
 print(proc.time()-ptm)
 
+# save(mean.time,med.time,mean.50,mean.100,file='Results/sims-Aug28.RData')
 
 ####### plot simulation data #######
-data = 634 # 774.5483, 634, 0.07252472, 0.1176774
-simdata = med.time
-xl = c(0.9*min(min(simdata),data),1.1*max(max(simdata),data))
-yl = c(0,max(hist(simdata,15)$density))
-hist(simdata,15,xlim=xl,ylim=yl,col=rgb(0,1,0,1/4),freq=F,
-     xlab='Neighbors Infected Within 100 Days', main=NULL)
+data = 0.1179039 # 778.2056, 638, 0.07248908, 0.1179039
+simdata = mean.100
+xl = c(0.95*min(min(simdata),data),1.05*max(max(simdata),data))
+yl = c(0,max(hist(simdata,25)$density))
+hist(simdata,25,xlim=xl,ylim=yl,col=rgb(0,1,0,1/4),freq=F,
+     xlab='Neighbor Infections Within 100 Days', main=NULL)
 abline(v=data,lwd=4,col='red')
+box()
 
 
 ##### inf dates histogram
-- 
cgit v1.2.3-70-g09d2