myprojects-hg-reborn

Fork

rm(list=ls())

library(Cairo)
library(igraph)
library(lattice)
library(ggplot2)
library(tikzDevice)


get_frequence <- function(contact){

freq <- table(contact)

#convert table into a matrix

freq <- matrix(c(as.numeric(names(freq)), freq), ncol=length(freq), byrow=TRUE, dimnames=NULL) 

return (freq)

}


histogram_creation <- function(data, bin_structure){
my_hist <- hist(data, plot=F, bin_structure)

return (my_hist)

}

log_binning <- function(x_min,x_max,n_breaks=10){
delta_log <- log(x_max/x_min)/(n_breaks-1)
my_seq<-seq(0,n_breaks-1)
log_breaks <- x_min*exp(my_seq*delta_log)
}


calc_I_infinity <- 1 #whether I should also read the data about the final number of infected
#individuals vs the first infected individual.

n_bin <- 20 # number of bins I want to use in the following

delta_slice <- 20 #duration (in sec) of my timeslice

save_networks <-  38  #38  #31  #60


time_period <- 78 #selected time_period I want to look at

infected_diam <- seq(1)  # seq(length(time_period))
infected_mean_degree<- seq(1)  #seq(length(time_period))

mean_degree_after_infection <- seq(1)  #seq(length(time_period))
diam_after_infection <- seq(1)  #seq(length(time_period))

dist_seed_last_infected <- seq(1)
time_dist_seed_last_infected <- seq(1)

interval_list <- -seq(1) #initialization that will be ovewritten

saturation_time_list <- list()

visitor_list <- read.table("visitor_list.dat", header=FALSE)
visitor_list <- as.matrix(visitor_list)

seed_entry_time <- read.table("seed_entry_time.dat", header=FALSE)

seed_entry_time <- as.matrix(seed_entry_time)

number_networks <- length(visitor_list)


infection_paths <- list()

infected_number_list <- list()


#map_list <- list()

g1_box <- list()

social_paths <- list()

social_number_list <- list()


g2_box <- list()


for (i in seq(number_networks)-1){

fn <- paste("dynamic_infectious_contacts_binary_",time_period,"_",i,"_.dat",sep="")

  
  print ("fn is, ")
  print (fn)

network_frame <- read.table(fn,header=FALSE)

g_infected <- graph.data.frame(network_frame,dir=FALSE)

g_infected <- simplify(g_infected)

infected_diam[i+1] <- length(get.diameter(g_infected,weights=NA))-1
infected_mean_degree[i+1] <- mean(degree(g_infected))

fn <- paste("dynamic_infectious_contacts_binary_and_time_",time_period,"_",i,"_.dat",sep="")

contacts_and_time <- read.table(fn,header=FALSE)
contacts_and_time <- as.matrix(contacts_and_time)

my_dim <- dim(contacts_and_time)

last_contact <- contacts_and_time[my_dim[1],3]

fn <- paste("network_after_infection_",time_period,"_",i,"_.dat",sep="")

network_frame_after_infection <- read.table(fn,header=FALSE)

g_after_infection <- graph.data.frame(network_frame_after_infection[ ,2:3],dir=FALSE)

g_after_infection <- simplify(g_after_infection)

diam_after_infection[i+1] <- length(get.diameter(g_after_infection,weights=NA))-1

mean_degree_after_infection[i+1] <- mean(degree(g_after_infection))


seed_inf <- visitor_list[i+1]

fn <- paste("dynamic_infectious_contacts_",time_period,"_",i,"_.dat",sep="")

infected_arr <- read.table(fn,header=FALSE)

infected_arr <- as.matrix(infected_arr)

last_infected <- as.numeric(tail(infected_arr,1))

#Now ensure this is really the last infected individual

if (length(which(infected_arr==last_infected))>1){

last_infected <- as.numeric(tail(infected_arr,2)[1])

}


if (last_infected==seed_inf){

last_infected <- as.numeric(tail(infected_arr,2)[1])

}

if (i==11){


isave <- last_infected

}

## fn <- paste("network_after_infection_unique",time_period,"_",i,"_.dat",sep="")

## network <- read.table(fn,header=FALSE)

## g1 <- graph.data.frame(network,dir=FALSE)

## g1 <- simplify(g1)

## V(g1)$number  <- as.numeric(strsplit(V(g1)$name, ","))

## network <- as.matrix(network)



g1 <-g_after_infection

V(g1)$number <- as.numeric(strsplit(V(g1)$name, ",")) #this converts a string into a numeric array 
V(g1)$id <- seq(vcount(g1))-1


g1_box[[i+1]] <- g1

seed_pos <- which(V(g1)$number==seed_inf)

last_infected_pos <- which(V(g1)$number==last_infected)

fn <- paste("dynamic_infectious_contacts_binary_unique_",time_period,"_",i,"_.dat",sep="")

  
network_inf <- read.table(fn,header=FALSE)

g2 <- graph.data.frame(network_inf,dir=FALSE)

g2 <- simplify(g2)

V(g2)$number <- as.numeric(strsplit(V(g2)$name, ","))
V(g2)$id <- seq(vcount(g2))-1


g2_box[[i+1]] <- g2

last_infected_pos_g2 <- which(V(g2)$number==last_infected)
seed_pos_g2 <- which(V(g2)$number==seed_inf)

mapping_nodes <- seq(vcount(g2))


for (m in seq(vcount(g2))){

mapping_nodes[m] <- which(V(g1)$number==V(g2)$number[m])-1


}

#map_list[[i+1]] <- mapping_nodes

#infection_paths is a list whose entries are other lists of shortest paths from the seed to the
#last 

infection_paths[[i+1]] <- get.shortest.paths(g2, from=(seed_pos_g2-1))

print("seed_pos -1 is, ")
print(seed_pos-1)

social_paths[[i+1]] <- get.shortest.paths(g1, from=(seed_pos-1), to=mapping_nodes) 



if (i==11){

gsave <- g2

li_save <- last_infected

s_save <- seed_inf

li_pos <- last_infected_pos_g2
s_pos <- seed_pos_g2

}



d3 <- get.shortest.paths(g2,(seed_pos_g2-1),(last_infected_pos_g2-1),weights=NA)

network_inf <- as.matrix(network_inf)




dist_seed_last_infected[i+1] <- length(d3[[1]])-1

time_dist_seed_last_infected[i+1] <- last_contact- seed_entry_time[i+1]








if (i==save_networks){



set.seed(1234)

pdf("Aggregated_and_infected_network.pdf")
par(bg="white")
l <- layout.fruchterman.reingold(g1)
l <- layout.norm(l, -1,1, -1,1) 

## l2 <- l[ d, ]

d <- get.diameter(g1,weights=NA)

V(g1)$color <- "#ff000033"
#V(g1)[seed_pos-1]$color <- "blue"

plot(g1, layout=l,
     #vertex.label.dist=0.5,   vertex.color="#ff000033",
     vertex.label.dist=0.5,
     vertex.frame.color="#ff000033", edge.color="#55555533",vertex.label=NA, vertex.size=4)





g1d <- subgraph(g1, d)
plot(g1d, layout=l[ V(g1d)$id+1, ],
#vertex.color="#ff000033", vertex.frame.color="#ff000033",
   edge.color="darkorange",vertex.label=NA, edge.width=2, vertex.size=4,
   add=TRUE, rescale=FALSE)


#The following 3 commented lines are Gabor"s suggestion. They generate a network equivalent to
#mine, but it does not look so good since it does change the overall graph layout.

last_infected_pos_g2 <- which(V(g2)$number==last_infected)
seed_pos_g2 <- which(V(g2)$number==seed_inf)



V(g2)$color <- "#ff0000"
V(g2)[seed_pos_g2-1]$color <- "blue"
V(g2)[last_infected_pos_g2-1]$color <- "black"


V(g2)$frame.color <- "#ff0000"
V(g2)[seed_pos_g2-1]$frame.color <- "blue"
V(g2)[last_infected_pos_g2-1]$frame.color <- "black"

l2 <- l[ V(g1)$name %in% V(g2)$name, ]
plot(g2, layout=l2,
#     vertex.color="#ff0000",
#     vertex.frame.color="#ff0000",
edge.color="#555555", vertex.label=NA, vertex.size=4, add=TRUE,rescale=FALSE)

d <- get.diameter(g2,weights=NA)




g2d <- subgraph(g2, d)

seed_pos_g2d <- which(V(g2d)$number==seed_inf)


V(g2d)[seed_pos_g2d-1]$color <- "blue"
V(g2d)[seed_pos_g2d-1]$frame.color <- "blue"

V(g2d)[last_infected_pos_g2-1]$color <- "black"
V(g2d)[last_infected_pos_g2-1]$frame.color <- "black"


plot(g2d, layout=l2[ V(g2d)$id+1, ],
#vertex.color="#ff0000",
   #  vertex.frame.color="#ff0000",
   edge.color="blue",vertex.label=NA, edge.width=2, vertex.size=4,
   add=TRUE, rescale=FALSE)




## d3 <- get.shortest.paths(g2,(seed_pos_g2-1),(last_infected_pos_g2-1),weights=NA)

## g3d <- subgraph(g2, d3[[1]])
## plot(g3d, layout=l2[ V(g3d)$id+1, ],
## #vertex.color="#ff0000", vertex.frame.color="#ff0000",
##    edge.color="darkgreen",vertex.label=NA, edge.width=2, vertex.size=4,
##    add=TRUE, rescale=FALSE)






dev.off()


set.seed(1234)


pdf("Aggregated_and_infected_network_2.pdf")
par(bg="white")
l <- layout.fruchterman.reingold(g1)
l <- layout.norm(l, -1,1, -1,1) 

## l2 <- l[ d, ]

d <- get.diameter(g1,weights=NA)

V(g1)$color <- "#ff000033"
#V(g1)[seed_pos-1]$color <- "blue"

plot(g1, layout=l,
     #vertex.label.dist=0.5,   vertex.color="#ff000033",
     vertex.label.dist=0.5,
     vertex.frame.color="#ff000033", edge.color="#55555533",vertex.label=NA, vertex.size=4)


V(g1)$id <- seq(vcount(g1))-1


## g1d <- subgraph(g1, d)
## plot(g1d, layout=l[ V(g1d)$id+1, ],
## #vertex.color="#ff000033", vertex.frame.color="#ff000033",
##    edge.color="darkorange",vertex.label=NA, edge.width=2, vertex.size=4,
##    add=TRUE, rescale=FALSE)


#The following 3 commented lines are Gabor"s suggestion. They generate a network equivalent to
#mine, but it does not look so good since it does change the overall graph layout.



V(g2)$color <- "#ff0000"
V(g2)[seed_pos_g2-1]$color <- "blue"
V(g2)[last_infected_pos_g2-1]$color <- "black"


V(g2)$frame.color <- "#ff0000"
V(g2)[seed_pos_g2-1]$frame.color <- "blue"
V(g2)[last_infected_pos_g2-1]$frame.color <- "black"



l2 <- l[ V(g1)$name %in% V(g2)$name, ]
plot(g2, layout=l2,
#     vertex.color="#ff0000",
#     vertex.frame.color="#ff0000",
edge.color="#555555", vertex.label=NA, vertex.size=4, add=TRUE,rescale=FALSE)

## d <- get.diameter(g2,weights=NA)

## V(g2)$id <- seq(vcount(g2))-1


## g2d <- subgraph(g2, d)

## seed_pos_g2d <- which(V(g2d)$number==seed_inf)


## V(g2d)[seed_pos_g2d-1]$color <- "blue"
## V(g2d)[seed_pos_g2d-1]$frame.color <- "blue"


## plot(g2d, layout=l2[ V(g2d)$id+1, ],
## #vertex.color="#ff0000",
##    #  vertex.frame.color="#ff0000",
##    edge.color="blue",vertex.label=NA, edge.width=2, vertex.size=4,
##    add=TRUE, rescale=FALSE)






g3d <- subgraph(g2, d3[[1]])
plot(g3d, layout=l2[ V(g3d)$id+1, ],
#vertex.color="#ff0000", vertex.frame.color="#ff0000",
   edge.color="#D40000",vertex.label=NA, edge.width=2, vertex.size=4,
   add=TRUE, rescale=FALSE)






dev.off()


pdf("Infected_network.pdf")

l <- layout.fruchterman.reingold(g2)
l <- layout.norm(l, -1,1, -1,1) 

## l2 <- l[ d, ]

V(g2)$color <- "#ff0000"
#V(g1)[select_nodes-1]$color <- 
E(g2)$color <- "#555555"
#E(g1)[select_edges-1]$color <- 


plot(g2, layout=l,
     vertex.label.dist=0.5, 
     vertex.frame.color="#ff000033", vertex.label=NA, vertex.size=4)


## V(g_patches)$id <- seq(vcount(g_patches))-1
## g2 <- subgraph(g_patches, d_patches)
## plot(g2, layout=l[ V(g2)$id+1, ],
## vertex.color="#ff0000", vertex.frame.color="#ff0000",
##    edge.color="#555555",vertex.label=NA, edge.width=2, vertex.size=4,
##    add=TRUE, rescale=FALSE)



dev.off()



  

}


fn <- paste("dynamic_infectious_times_unique_",time_period,"_",i,"_.dat",sep="")

time_list <- read.table(fn,header=FALSE)
time_list <- as.matrix(time_list)

saturation_time_list[[i+1]] <- time_list

time_list <- diff(time_list) #now I actually have a list of intervals between consecutive infections

interval_list <- c(interval_list,time_list)

}



interval_list <- interval_list[2:length(interval_list)] #get rid of the initial negative value

data_interval <- get_frequence(interval_list)


pdf("Interval_distribution_network_infectious_contacts.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot( 
data_interval[1, ]/60.,data_interval[2, ]/sum(data_interval[2, ]), "p" , col="blue",ylab="P(delta t)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="delta t between infections [min]", main="Infectious Network",
     lwd=2)
dev.off()



pdf("Interval_distribution_network_infectious_contacts-log-log.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot( 
data_interval[1, ]/60.,data_interval[2, ]/sum(data_interval[2, ]), "p" ,log="xy", col="blue",ylab="P(delta t)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="delta t between infections [min]", main="Infectious Network",
     lwd=2)
dev.off()


#Now I create a histogram

interval_grid <- seq(min(data_interval[1, ]/60.), max(data_interval[1, ]/60.),length=n_bin)

h <- histogram_creation(interval_list/60., interval_grid)

pdf("Interval_distribution_network_infectious_contacts_binned_semilog.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
i <- seq(1,length(h$breaks)-1)
plot(h$mids,h$counts/sum(h$counts), "b" ,log="y", col="blue",
     ylab="P(delta t)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="delta t between infections [min]", 
     lwd=2)
rect(h$breaks[i],min(h$counts[which(h$counts>0)]/sum(h$counts )) , h$breaks[i+1],
     h$counts/sum(h$counts ))
dev.off()



#Define a new grid

interval_grid <-log_binning(min(data_interval[1, ]/60.), max(data_interval[1, ]/60.),n_bin)

h <- histogram_creation(interval_list/60., interval_grid)


pdf("Interval_distribution_network_infectious_contacts_log_binning.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
i <- seq(1,length(h$breaks)-1)
plot(h$mids,h$counts/sum(h$counts), "b" ,log="xy", col="blue",
     ylab="P(delta t)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="delta t between infections [min]", 
     lwd=2)
rect(h$breaks[i],min(h$counts[which(h$counts>0)]/sum(h$counts )) , h$breaks[i+1],
     h$counts/sum(h$counts ))
dev.off()




data_infection <- as.data.frame(infected_diam)
names(data_infection) <- "inf_diam"

gpl <- ggplot(data_infection,aes(x=inf_diam,y=..count../sum(..count..)) )+
  geom_histogram(binwidth=2, colour="darkgreen", fill="white")+
  ## scale_x_log10() +
  ## scale_y_log10()+
  #scale_y_continuous(limits=c(0,.3))+
  ## scale_y_continuous(trans = "log10")+

  xlab(expression(paste("Diameter Length")))+
  ylab(expression(paste("P(diameter length)")))+
  opts(title = expression(paste("Infection Network in Dublin: July, ",14^th)))+
  opts(plot.title = theme_text(size = 20))+
  opts(axis.title.x = theme_text(size = 20))+
  opts(axis.title.y = theme_text(size = 20, angle=90))+
  opts(axis.text.x = theme_text(size=15, colour="gray",vjust=1))+
  opts(axis.text.y = theme_text(size=15, colour="gray", hjust=1))



fn <- paste("histogram-infection-network-diameter.pdf")

CairoPDF(fn)
print(gpl)
dev.off()





data_infected_diam <- get_frequence(infected_diam)



pdf("Diameter_distribution_network_infectious_contacts.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot( 
data_infected_diam[1, ],data_infected_diam[2, ]/sum(data_infected_diam[2, ]), "p" , col="blue",ylab="P(diameter length)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="diameter length", main="Infectious Network",
     lwd=2)
dev.off()




data_infection <- as.data.frame(diam_after_infection)
names(data_infection) <- "inf_diam"

gpl <- ggplot(data_infection,aes(x=inf_diam,y=..count../sum(..count..)) )+
  geom_histogram(binwidth=2, colour="darkgreen", fill="white")+
  ## scale_x_log10() +
  ## scale_y_log10()+
  #scale_y_continuous(limits=c(0,.3))+
  ## scale_y_continuous(trans = "log10")+

  xlab(expression(paste("Diameter Length")))+
  ylab(expression(paste("P(diameter length)")))+
  opts(title = expression(paste("Aggregated Network in Dublin: July, ",14^th)))+
  opts(plot.title = theme_text(size = 20))+
  opts(axis.title.x = theme_text(size = 20))+
  opts(axis.title.y = theme_text(size = 20, angle=90))+
  opts(axis.text.x = theme_text(size=15, colour="gray",vjust=1))+
  opts(axis.text.y = theme_text(size=15, colour="gray", hjust=1))



fn <- paste("histogram-aggregated-network-diameter.pdf")

pdf(fn)
print(gpl)
dev.off()



data_infected_diam <- get_frequence(diam_after_infection)



pdf("Diameter_distribution_network_social_contacts_after_infection.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot( 
data_infected_diam[1, ],data_infected_diam[2, ]/sum(data_infected_diam[2, ]), "p" , col="blue",ylab="P(diameter length)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="diameter length", main="Network after infection",
     lwd=2)
dev.off()



k_grid <- seq(min(infected_mean_degree), max(infected_mean_degree), length=n_bin)

h <- hist(infected_mean_degree, plot=F, k_grid)





pdf("Distribution_k_aver_infected_binned.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
i <- seq(1,length(k_grid)-1)
plot(h$mids,h$counts/sum(h$counts), "b" , col="blue",
     ylab="P(<k>)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="<k>", 
     lwd=2,main="Infectious Network")
rect(h$breaks[i], 0, h$breaks[i+1],h$counts/sum(h$counts ))
dev.off()


k_grid <- seq(min(mean_degree_after_infection), max(mean_degree_after_infection), length=n_bin)

h <- hist(mean_degree_after_infection, plot=F, k_grid)





pdf("Distribution_k_aver_after_infection_binned.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
i <- seq(1,length(k_grid)-1)
plot(h$mids,h$counts/sum(h$counts), "b" , col="blue",
     ylab="P(<k>)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="<k>", 
     lwd=2,main="Network after infection")
rect(h$breaks[i], 0, h$breaks[i+1],h$counts/sum(h$counts ))
dev.off()




if (calc_I_infinity==1){

I_infinity_distr <-   read.table("infected_number_distr.dat", header=FALSE)
I_infinity_distr <- as.matrix(I_infinity_distr)

I_freq <- get_frequence(I_infinity_distr)

pdf("Distribution_number_infected.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
plot( 
I_freq[1, ],I_freq[2, ]/sum(I_freq[2, ]), "p" , col="blue",ylab="P(# infected individuals)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="# infected individuals", main="Infectious Network",
     lwd=2)
dev.off()


number_grid <- seq(min(I_infinity_distr), max(I_infinity_distr), length=n_bin)

h <- hist(I_infinity_distr, plot=F, number_grid)





pdf("Distribution_number_infected_binned.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)
i <- seq(1,length(number_grid)-1)
plot(h$mids,h$counts/sum(h$counts), "b" , col="blue",
     ylab="P(# infected individuals)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="# infected individuals", 
     lwd=2,main="Infectious Network")
rect(h$breaks[i], 0, h$breaks[i+1],h$counts/sum(h$counts ))
dev.off()


}

#Now some work on the saturation time list

last_infection_time <- -1 #initialize to impossible value 

for (i in seq(length(saturation_time_list))){

  last_infection_time <- max(last_infection_time,max(saturation_time_list[[i]]))
#saturation_time_list[[i]] 

}

time_seq <- seq(0,last_infection_time,by=delta_slice)

progress_infection <- seq(length(saturation_time_list)*length(time_seq))

dim(progress_infection) <- c(length(saturation_time_list),length(time_seq))


#Now calculate all the trajectories giving the spread of the infection

for (i in seq(length(saturation_time_list))){
  print ("i is, ")
  print(i)
for (j in seq(length(time_seq))){
progress_infection[i,j] <- length(which(saturation_time_list[[i]]<=time_seq[j]))

}
}


aver_infection <- seq(length(time_seq))
std_infection <- seq(length(time_seq))


for (i in seq(length(time_seq))){
aver_infection[i] <-mean(progress_infection[ ,i]) 
std_infection[i] <-sqrt(var(progress_infection[ ,i])) 


}

write.table(aver_infection, "average_number_infected_individuals.dat",col.names=FALSE,row.names=FALSE )
write.table(time_seq, "time_sequence_for_average_number_infected_individuals.dat",col.names=FALSE,row.names=FALSE )
write.table(std_infection, "standard_deviation_number_infected_individuals.dat",col.names=FALSE,row.names=FALSE )



write.table(progress_infection, "number_infected_individuals_all_trajectories.dat",col.names=FALSE,row.names=FALSE )



pdf("Progress_average_infection.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot(time_seq/60.,aver_infection, "b" , col="blue",
     ylab="Average number infected individuals",
     ,cex.axis=1.4,cex.lab=1.6, xlab="Time from first new infection [min]", 
     lwd=2)

dev.off()

pdf("Progress_average_infection_normalized.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot(time_seq/60.,aver_infection/max(aver_infection), "b" , col="blue",
     ylab="Average number infected individuals (norm)",
     ,cex.axis=1.4,cex.lab=1.6, xlab="Time from first new infection [min]", 
     lwd=2)

dev.off()



xx <- c(time_seq/60., rev(time_seq/60.))
yy <- c((aver_infection+std_infection)/max(aver_infection),rev(aver_infection-std_infection)/max(aver_infection)) 

pdf("shaded_normalized.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot(xx ,yy,type="n", col="blue",ylab="Normalized # infected",
     ,cex.axis=1.4,cex.lab=1.6, xlab="Time from first new infection [min]", main="Dublin")

polygon(xx, yy, col="gray", border = "black")

lines(time_seq/60.,aver_infection/max(aver_infection) , col="blue", lwd=2, lty=6)

#lines(time_inf_ht,aver_inf_ht-std_inf_ht , col="black", lwd=2, lty=1)


## legend("bottomright",cex=1.2, c(expression("Dublin"),
##                    expression("HT2009")), 
## lwd=c(2,2), lty=c(2,6),pch = c(-1,-1),col=c("blue","red"),box.lwd=0,box.lty=0,
## ,xjust = 1, yjust = 1)

legend("bottomright",cex=1.2, c(expression("Normalized Mean Number of Infected")) 
,lwd=c(2), lty=c(6),pch = c(-1),col=c("blue"),box.lwd=0,box.lty=0,
,xjust = 1, yjust = 1)



dev.off()

yy <- c((aver_infection+std_infection),rev(aver_infection-std_infection)) 

pdf("shaded.pdf")
par( mar = c(4.5,5, 2, 1) + 0.1)

plot(xx ,yy,type="n", col="blue",ylab="Number of infected",
     ,cex.axis=1.4,cex.lab=1.6, xlab="Time from first new infection [min]", main="Dublin")

polygon(xx, yy, col="gray", border = "black")

lines(time_seq/60.,aver_infection , col="blue", lwd=2, lty=6)

#lines(time_inf_ht,aver_inf_ht-std_inf_ht , col="black", lwd=2, lty=1)


## legend("bottomright",cex=1.2, c(expression("Dublin"),
##                    expression("HT2009")), 
## lwd=c(2,2), lty=c(2,6),pch = c(-1,-1),col=c("blue","red"),box.lwd=0,box.lty=0,
## ,xjust = 1, yjust = 1)

legend("bottomright",cex=1.2, c(expression("Mean Number of Infected")) 
,lwd=c(2), lty=c(6),pch = c(-1),col=c("blue"),box.lwd=0,box.lty=0,
,xjust = 1, yjust = 1)



dev.off()

###########################################################
###########################################################


data_new <- cbind(as.data.frame(time_dist_seed_last_infected/60.),as.data.frame(dist_seed_last_infected))

names(data_new) <- c("time_dist","net_dist")


p_cor <- cor(data_new$time_dist,data_new$net_dist, method='pearson')

print("p_cor is, ")
print(p_cor)


write.table(data_new, "time_and_network_dist.dat",row.names=FALSE )

gpl <- ggplot(data_new,aes(x=time_dist,y=net_dist)) +
  geom_point(size=3, colour="blue")+
   opts(panel.background=theme_rect(fill="white"))+
  ## scale_x_log10() +
  ## scale_y_log10()+
  scale_y_continuous(limits=c(0,18), breaks=seq(0,18, by=2))+
  scale_x_continuous(limits=c(0,500), breaks=seq(0,500,by=100))+
  opts(panel.grid.minor = theme_blank())+
opts(panel.grid.major = theme_blank())+
  xlab(expression(paste("Time Distance (min)")))+
  ylab(expression(paste("Transmission ", n[d])))+
  opts(title = expression(paste("Dublin: July, ",14^th)))+
  opts(plot.title = theme_text(size = 35))+
  opts(axis.title.x = theme_text(size = 35))+
  opts(axis.title.y = theme_text(size = 35, angle=90))+
    opts(axis.ticks = theme_segment(colour = "black", size=0.5),axis.ticks.length = unit(0.3, "cm"))+
  opts(axis.text.x = theme_text(size=30, colour="black",vjust=1))+
  opts(axis.text.y = theme_text(size=30, colour="black", hjust=1))



fn <- paste("scatterplot_time_and_network_distance.pdf")

CairoPDF(fn)
print(gpl)
dev.off()


gpl2 <- ggplot(data_new,aes(x=time_dist,y=net_dist)) +
  geom_point(size=2, colour="blue")+
   opts(panel.background=theme_rect(fill="white",size=1.5))+
  ## scale_x_log10() +
  ## scale_y_log10()+
  scale_y_continuous(limits=c(0,18), breaks=seq(0,18, by=2))+
  scale_x_continuous(limits=c(0,500), breaks=seq(0,500,by=100))+
  opts(panel.grid.minor = theme_blank())+
opts(panel.grid.major = theme_blank())+
  xlab(expression(paste("Time Distance (min)")))+
  ylab(expression(paste("Transmission  $n_d$")))+
  opts(title = expression(paste("SG: July, $14^{\\rm th}$")))+
  opts(plot.title = theme_text(size = 25))+
  opts(axis.title.x = theme_text(size = 20))+
  opts(axis.title.y = theme_text(size = 20, angle=90))+
    opts(axis.ticks = theme_segment(colour = "black", size=1),axis.ticks.length = unit(0.3, "cm"))+
  opts(axis.text.x = theme_text(size=18, colour="black",vjust=1))+
  opts(axis.text.y = theme_text(size=18, colour="black", hjust=1))

tikz('scatterplot_time_and_network_distance_tik.tex', standAlone = TRUE, width=5,height=5)
print(gpl2)
dev.off() 

tools::texi2dvi('scatterplot_time_and_network_distance_tik.tex',pdf=T)




################################################################

count <- 1



for (i in seq(length(infection_paths))){

for (j in seq(length(infection_paths[[i]]))){

v_numbers <- unlist(infection_paths[[i]][j])+1

v_numbers_2 <- unlist(social_paths[[i]][j])+1

infected_number_list[[count]] <- V(g2_box[[i]])$number[v_numbers]
social_number_list[[count]] <- V(g1_box[[i]])$number[v_numbers_2]

count <- count+1

}


}


length_inf <- seq(1)
length_soc <- seq(1)
overlap <- seq(1)

count <- 1

for (i in seq(length(social_number_list))){
length_inf[i] <- length(infected_number_list[[i]])-1
length_soc[i] <- length(social_number_list[[i]])-1

if (length_inf[i]>0){
my_inf <- unlist(infected_number_list[[i]])
my_soc <- unlist(social_number_list[[i]])

overlap[count] <- length(intersect(my_inf, my_soc))/(length_soc[i]+1)

  count <- count+1

}


}

over_data <- as.data.frame(overlap)

names(over_data) <- c("overlap")

gpl <- ggplot(over_data,aes(x=overlap))+
geom_histogram(aes(y=..count../sum(..count..)),binwidth=0.05, colour="black", fill="blue",origin=0)+
  opts( panel.background=theme_rect(fill='white'))+
    opts(panel.grid.minor = theme_blank())+
opts(panel.grid.major = theme_blank())+
   opts(axis.ticks = theme_segment(colour = "black", size=0.5),axis.ticks.length = unit(0.3, "cm"))+
scale_x_continuous(limits=c(0.2,1),breaks=seq(0.2,1,by=0.2))+
scale_y_continuous( limits=c(0,0.3), breaks=seq(0,.3,by=.1))+
##scale_y_continuous(trans = "log10")+
## scale_x_log10()+
## scale_y_log10()+
xlab("Overlap")+
ylab("P(overlap)")+
opts(axis.title.x = theme_text(size = 35))+
opts(axis.title.y = theme_text(size = 35, angle=90, hjust=0.3))+
  opts(title = expression(paste("Dublin, July,", 14^th)))+
  opts(plot.title = theme_text(size = 35))+
opts(axis.text.x = theme_text(size=30, colour="black",vjust=1))+
opts(axis.text.y = theme_text(size=30, colour="black", hjust=1))


CairoPDF("overlap_distribution.pdf")
print(gpl)
dev.off()



sel <- which(length_inf>0)
length_inf <- length_inf[sel]
length_soc <- length_soc[sel]


freq_inf <- t(get_frequence(length_inf))
freq_soc <- t(get_frequence(length_soc))

freq_inf[,2] <- freq_inf[,2]/sum(freq_inf[,2])
freq_soc[,2] <- freq_soc[,2]/sum(freq_soc[,2])

freq_inf <- as.data.frame(freq_inf)
freq_soc <- as.data.frame(freq_soc)


names(freq_inf) <- c("length", "count")
names(freq_soc) <- c("length", "count")

data_path <- rbind(transform(freq_inf, set="Transmission network"), transform(freq_soc, set="Aggregated network"))


lbls <- levels(data_path$set)



gpl <- ggplot(data_path, aes(x=length, y=count, colour=set, shape=set)) +
 ## geom_vline(xintercept = seq(0,80,by=20), col="white", size=1)+
 ## geom_vline(xintercept = seq(0,80,by=10), col="white", size=0.5)+
 ## geom_hline(yintercept = seq(0,25,by=5), col="white", size=1)+
 ## geom_hline(yintercept = seq(0,25,by=2.5), col="white", size=0.5)+
  geom_point(size=4.5) +
#geom_line() +
      opts( panel.background=theme_rect(fill='white'))+
#scale_colour_hue("Data from", breaks=lbls, labels=lbls) +
scale_shape_manual("Network", breaks=lbls, labels=lbls, values=c(1,2)) +
scale_colour_manual("Network", breaks=lbls, labels=lbls, values=c(1,2)) +

  #extra stuff
  xlab(expression(paste(n[d])))+
#  ylab(expression(paste("<s>/k")))+
#  ylab(expression(""*symbol("\341")*s*symbol("\361")/k*""))+
 ylab(expression(paste("P(",n[d],")")))+
  scale_x_continuous(limits=c(1,19), breaks=seq(1,19, by=2))+
scale_y_continuous(limits=c(0,.3), breaks=seq(0,.3,by=.1))+
opts(title = expression(paste("Dublin: July, ",14^th)))+
   opts(plot.title = theme_text(size = 35))+
  opts(axis.title.x = theme_text(size = 35))+
  opts(legend.text = theme_text(size = 14, vjust=0.8))+
   opts(legend.title = theme_text(size = 18,hjust=0))+
   #opts(legend.title = theme_blank())+
    opts(panel.grid.minor = theme_blank())+
opts(panel.grid.major = theme_blank())+
   opts(axis.ticks = theme_segment(colour = "black", size=0.5),axis.ticks.length = unit(0.3, "cm"))+
  opts(legend.position = c(0.65, 0.7), legend.background=theme_rect(fill='white'),
       legend.key = theme_rect(colour = NA))+
  opts(axis.title.y = theme_text(size = 35, angle=90))+
  opts(axis.text.x = theme_text(size=30, colour="black",vjust=1))+
  opts(axis.text.y = theme_text(size=30, colour="black", hjust=1))


fn <- paste("path_comparison.pdf")

#pdf(fn)
CairoPDF(fn)
print(gpl)
dev.off()


gpl2 <- ggplot(data_path, aes(x=length, y=count, colour=set, shape=set)) +
 ## geom_vline(xintercept = seq(0,80,by=20), col="white", size=1)+
 ## geom_vline(xintercept = seq(0,80,by=10), col="white", size=0.5)+
 ## geom_hline(yintercept = seq(0,25,by=5), col="white", size=1)+
 ## geom_hline(yintercept = seq(0,25,by=2.5), col="white", size=0.5)+
  geom_point(size=4) +
    geom_point(size=3.8) +
    geom_point(size=3.6) +
#geom_line() +
      opts( panel.background=theme_rect(fill='white',size=1.5))+
#scale_colour_hue("Data from", breaks=lbls, labels=lbls) +
scale_shape_manual("", breaks=lbls, labels=lbls, values=c(1,2)) +
scale_colour_manual("", breaks=lbls, labels=lbls, values=c(1,2)) +

  #extra stuff
  xlab(expression(paste("$n_d$")))+
#  ylab(expression(paste("<s>/k")))+
#  ylab(expression(""*symbol("\341")*s*symbol("\361")/k*""))+
 ylab(expression(paste("$P(n_d)$")))+
  scale_x_continuous(limits=c(1,19), breaks=seq(1,19, by=2))+
scale_y_continuous(limits=c(0,.3), breaks=seq(0,.3,by=.1))+
opts(title = expression(paste("SG: July, $14^{\\rm th}$")))+
   opts(plot.title = theme_text(size = 25))+
  opts(axis.title.x = theme_text(size = 20))+
  opts(legend.text = theme_text(size = 14, vjust=0.4))+
   opts(legend.title = theme_text(size = 18,hjust=0))+
   #opts(legend.title = theme_blank())+
    opts(panel.grid.minor = theme_blank())+
opts(panel.grid.major = theme_blank())+
   opts(axis.ticks = theme_segment(colour = "black", size=1),axis.ticks.length = unit(0.15, "cm"))+
  opts(legend.position = c(0.65, 0.7), legend.background=theme_rect(col=0),
       legend.key = theme_rect(colour = NA))+
  opts(axis.title.y = theme_text(size = 20, angle=90))+
  opts(axis.text.x = theme_text(size=18, colour="black",vjust=1))+
  opts(axis.text.y = theme_text(size=18, colour="black", hjust=1))

tikz('path_comparison_tik.tex', standAlone = TRUE, width=5,height=5)
print(gpl2)
dev.off() 

tools::texi2dvi('path_comparison_tik.tex',pdf=T)



print("So far so good")