## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----setupMatrix2, echo=FALSE-------------------------------------------------
matrix(c(0, 0.2, 0.4, 
         0.5, 0, 0.6, 
         0.5, 0.8, 0),
       nrow = 3, ncol = 3,
       dimnames = list(c("A", "B", "C"), c("A", "B", "C")))

## ----setupMatrix, echo=FALSE, message=FALSE, warning=FALSE--------------------
if (!(requireNamespace("ggplot2", quietly = TRUE) &&
      requireNamespace("viridis", quietly = TRUE) &&
      requireNamespace("igraph", quietly = TRUE)  &&
      requireNamespace("ggnetwork", quietly = TRUE))) {
  message("Packages 'ggplot2', 'viridis', 'igraph' and 'ggnetwork'  are needed for plotting this figure.")
} else {
  library(ggplot2)
  library(viridis)
  library(igraph)
  library(ggnetwork)
  
  transition.matrix <- matrix(c(0, 0.2, 0.4, 0.5, 0, 0.6, 0.5, 0.8, 0),
                              nrow = 3,
                              ncol = 3,
                              dimnames = list(c("A", "B", "C"),
                                              c("A", "B", "C")))
  
  # melting the matrix go get from -> to in one line with probability
  melted.transition.matrix <- as.data.frame.table(transition.matrix,
                                                  stringsAsFactors = FALSE)
  colnames(melted.transition.matrix) <- c("from", "to", "prob")
  
  melted.transition.matrix <- subset(melted.transition.matrix,
                                     melted.transition.matrix$prob != 0)
  
  graph.Matrix <- igraph::graph.data.frame(melted.transition.matrix,
                                           directed = TRUE)
  
  graph.Matrix2 = igraph::layout_in_circle(graph.Matrix)
  
  
  if (utils::packageVersion("ggnetwork") > "0.5.12") {
    # unfortunately, the update to igraph 2.0.0 has broken the connection
    # between ggnetwork and igraph, see: 
    # https://github.com/briatte/ggnetwork/issues/74
    
    #using ggnetwork to provide the layout
    graph.matrix.network <- ggnetwork::ggnetwork(x = graph.Matrix,
                                                 layout = graph.Matrix2,
                                                 arrow.gap = 0.18)
    
    graph.matrix.network2 <- subset(graph.matrix.network,
                                    !is.na(graph.matrix.network$prob))
    
    #plotting the network
    ggplot(graph.Matrix, aes(x = x, y = y, xend = xend, yend = yend)) +
      geom_edges(color = "grey70",
                 arrow = arrow(length = unit(1, "lines"),
                               type = "closed"),
                 curvature = 0.2) +
      geom_nodes(aes(color = name) , size = 30) +
      geom_nodetext(aes(label = name),
                    color = "white",
                    fontface = "bold",
                    size = 15) +
      scale_color_viridis(guide = FALSE, discrete = TRUE) +
      theme_blank() +
      ylim(-0.5, 1.2) +
      xlim(-0.5, 1.2) +
      geom_text(data = graph.matrix.network2,
                aes(x = (x + xend) / 2,
                    y = (y + yend) / 2,
                    label = prob, 
                    color = name),
                size = 6)
  } else {
    plot(graph.Matrix,
         edge.label = melted.transition.matrix$prob,
         vertex.color = "white",
         edge.curved = seq(-0.5, 0.5, length = 1 + ecount(graph.Matrix)), 
         arrow.mode = 3,
         edge.arrow.size = 0.3)
  }
}

## ----setupA, eval=FALSE-------------------------------------------------------
#  SimulationSingle <- nosoiSim(type="single", popStructure="discrete", ...)

## ----setupB, eval=FALSE-------------------------------------------------------
#  SimulationSingle <- nosoiSim(type="single", popStructure="none",
#                               length.sim=300, max.infected=1000, init.individuals=1, ...)

## ----setupA-dual, eval=FALSE--------------------------------------------------
#  SimulationDual <- nosoiSim(type="dual", popStructure="discrete", ...)

## ----pExit1, eval=FALSE-------------------------------------------------------
#  p_Exit_fct  <- function(t,current.in){
#    if(current.in=="A"){return(0.02)}
#    if(current.in=="B"){return(0.05)}
#    if(current.in=="C"){return(0.1)}
#  }

## ----pMove1, eval=FALSE-------------------------------------------------------
#  p_Move_fct  <- function(t){return(0.1)}

## ----nContact1, eval=FALSE----------------------------------------------------
#  n_contact_fct = function(t){abs(round(rnorm(1, 0.5, 1), 0))}

## ----nContact2, echo=FALSE, message=FALSE-------------------------------------
if (!(requireNamespace("ggplot2", quietly = TRUE) &&
      requireNamespace("dplyr", quietly = TRUE))) {
  message("Packages 'ggplot2' and 'dplyr' are needed for plotting this figure.")
} else {
  library(ggplot2)
  library(dplyr)
  set.seed(900)
  data = data.frame(N=abs(round(rnorm(200, 0.5, 1), 0)))
  
  data = data %>% group_by(N) %>% summarise(freq=length(N)/200)
  
  ggplot(data=data, aes(x=as.factor(N), y=freq)) + geom_bar(stat="identity") + theme_minimal() + labs(x="nContact",y="Frequency")
}

## ----pTrans1, eval=FALSE------------------------------------------------------
#  p_Trans_fct <- function(t, p_max, t_incub){
#    if(t < t_incub){p=0}
#    if(t >= t_incub){p=p_max}
#    return(p)
#  }

## ----pTrans2, eval=FALSE------------------------------------------------------
#  t_incub_fct <- function(x){rnorm(x, mean=7, sd=1)}
#  p_max_fct <- function(x){rbeta(x, shape1=5, shape2=2)}

## ----pTrans3, echo=FALSE------------------------------------------------------
if (!(requireNamespace("ggplot2", quietly = TRUE) &&
      requireNamespace("dplyr", quietly = TRUE))) {
  message("Packages 'ggplot2' and 'dplyr' are needed for plotting this figure.")
} else {
  library(ggplot2)
  library(dplyr)
  
  set.seed(506)
  
  p_Trans_fct <- function(t, p_max, t_incub){
    if(t < t_incub){p=0}
    if(t >= t_incub){p=p_max}
    return(p)
  }
  
  t_incub_fct <- function(x){rnorm(x,mean = 7,sd=1)}
  p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}
  
  data = data.frame(t_incub = t_incub_fct(200),
                    p_max = p_max_fct(200),
                    host = paste0("H-", 1:200))
  
  t=c(0:12)
  data3=NULL
  for(t in 0:15){
    data2 = data %>% group_by(host) %>% 
      mutate(proba = p_Trans_fct(t = t, p_max = p_max, t_incub = t_incub))
    data2$t = t
    data3 = rbind(data3, data2)
  }
  
  ggplot(data = data3, aes(x = t, y = proba, group = host)) + 
    geom_line(color = "grey60") +
    theme_minimal() +
    labs(x = "Time since infection (t)", y = "pTrans")
}

## ----pTrans4, eval=FALSE------------------------------------------------------
#  t_incub_fct <- function(x){rnorm(x,mean = 7,sd=1)}
#  p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}
#  
#  param_pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct)

## ----setupF-------------------------------------------------------------------
library(nosoi)

#Transition matrix
transition.matrix <- matrix(c(0,0.2,0.4,0.5,0,0.6,0.5,0.8,0),nrow = 3, ncol = 3,dimnames=list(c("A","B","C"),c("A","B","C")))

#pExit
p_Exit_fct  <- function(t,current.in){
  if(current.in=="A"){return(0.02)}
  if(current.in=="B"){return(0.05)}
  if(current.in=="C"){return(0.1)}
}

#pMove
p_Move_fct  <- function(t){return(0.1)}

#nContact
n_contact_fct = function(t){abs(round(rnorm(1, 0.5, 1), 0))}

#pTrans
proba <- function(t,p_max,t_incub){
  if(t <= t_incub){p=0}
  if(t >= t_incub){p=p_max}
  return(p)
}

t_incub_fct <- function(x){rnorm(x,mean = 5,sd=1)}
p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}

param_pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct)

# Starting the simulation ------------------------------------

set.seed(846)
SimulationSingle <- nosoiSim(type="single", popStructure="discrete",
                             length.sim=300, max.infected=300, init.individuals=1, init.structure="A", 
                             
                             structure.matrix=transition.matrix,
                             
                             pExit = p_Exit_fct,
                             param.pExit=NA,
                             timeDep.pExit=FALSE,
                             diff.pExit=TRUE,
                             
                             pMove = p_Move_fct,
                             param.pMove=NA,
                             timeDep.pMove=FALSE,
                             diff.pMove=FALSE,
                             
                             nContact=n_contact_fct,
                             param.nContact=NA,
                             timeDep.nContact=FALSE,
                             diff.nContact=FALSE,
                             
                             pTrans = proba,
                             param.pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct),
                             timeDep.pTrans=FALSE,
                             diff.pTrans=FALSE,
                             
                             prefix.host="H",
                             print.progress=FALSE,
                             print.step=10)

## ----pExit1-dual, eval=FALSE--------------------------------------------------
#  p_Exit_fctB  <- function(t,prestime){(sin(prestime/(2*pi*10))+1)/16} #for a periodic function

## ----pExit2-dual, echo=FALSE--------------------------------------------------
p_Exit_fctx <- function(x){(sin(x/(2*pi*10))+1)/16} #for a periodic function
if (!requireNamespace("ggplot2", quietly = TRUE)) {
  message("Package 'ggplot2' is needed for plotting this figure.")
} else {
  ggplot(data=data.frame(x=0), aes(x=x)) + stat_function(fun=p_Exit_fctx) + theme_minimal() + labs(x="Absolute time (prestime)",y="pExit") + xlim(0,360)
}

## ----pMove.B, eval=FALSE------------------------------------------------------
#  p_Move_fct.B  <- NA

## ----nContact1.B, eval=FALSE--------------------------------------------------
#  n_contact_fct.B = function(t){sample(c(0,1,2),1,prob=c(0.6,0.3,0.1))}

## ----nContact2.B, echo=FALSE--------------------------------------------------
if (!(requireNamespace("ggplot2", quietly = TRUE) && 
      requireNamespace("dplyr", quietly = TRUE))) {
  message("Packages 'ggplot2' and 'dplyr' are needed for plotting this figure.")
} else {
  library(ggplot2)
  library(dplyr)
  set.seed(6059)
  data = data.frame(N=sample(c(0,1,2),200,replace=TRUE,prob=c(0.6,0.3,0.1)))
  
  data = data %>% group_by(N) %>% summarise(freq=length(N)/200)
  
  ggplot(data=data, aes(x=as.factor(N), y=freq)) + geom_bar(stat="identity") + theme_minimal() + labs(x="nContact.B",y="Frequency")
}

## ----pTrans1.B, eval=FALSE----------------------------------------------------
#  p_Trans_fct.B <- function(t, max.time){
#    dnorm(t, mean=max.time, sd=2)*5
#  }

## ----pTrans2.B, eval=FALSE----------------------------------------------------
#  max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}

## ----pTrans3.B, echo=FALSE----------------------------------------------------
if (!(requireNamespace("ggplot2", quietly = TRUE) && 
      requireNamespace("dplyr", quietly = TRUE))) {
  message("Packages 'ggplot2' and 'dplyr' are needed for plotting this figure.")
} else {
  library(ggplot2)
  library(dplyr)
  set.seed(6609)
  
  p_Trans_fct <- function(t, max.time){
    dnorm(t, mean=max.time, sd=2)*5
  }
  
  max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}
  
  data = data.frame(max.time=max.time_fct(200),host=paste0("H-",1:200))
  
  t=c(0:12)
  data3=NULL
  for(t in 0:15){
    data2 = data %>% group_by(host) %>% mutate(proba=p_Trans_fct(t=t,max.time=max.time))
    data2$t = t
    data3 = rbind(data3, data2)
  }
  
  ggplot(data=data3, aes(x=t, y=proba,group=host)) + geom_line(color="grey60",alpha=0.3) + theme_minimal() + labs(x="Time since infection (t)",y="pTrans")
}

## ----pTrans4.B, eval=FALSE----------------------------------------------------
#  max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}
#  
#  param_pTrans.B = list(max.time=max.time_fct)

## ----setupF.B-----------------------------------------------------------------
library(nosoi)

#Transition matrix
transition.matrix <- matrix(c(0,0.2,0.4,0.5,0,0.6,0.5,0.8,0),nrow = 3, ncol = 3,dimnames=list(c("A","B","C"),c("A","B","C")))

#Host A -----------------------------------

#pExit
p_Exit_fct  <- function(t,current.in){
  if(current.in=="A"){return(0.02)}
  if(current.in=="B"){return(0.05)}
  if(current.in=="C"){return(0.1)}
}

#pMove
p_Move_fct  <- function(t){return(0.1)}

#nContact
n_contact_fct = function(t){abs(round(rnorm(1, 0.5, 1), 0))}

#pTrans
proba <- function(t,p_max,t_incub){
  if(t <= t_incub){p=0}
  if(t >= t_incub){p=p_max}
  return(p)
}

t_incub_fct <- function(x){rnorm(x,mean = 5,sd=1)}
p_max_fct <- function(x){rbeta(x,shape1 = 5,shape2=2)}

param_pTrans = list(p_max=p_max_fct,t_incub=t_incub_fct)

#Host B -----------------------------------

#pExit
p_Exit_fct.B  <- function(t,prestime){(sin(prestime/(2*pi*10))+1)/16}

#pMove
p_Move_fct.B  <- NA

#nContact
n_contact_fct.B = function(t){sample(c(0,1,2),1,prob=c(0.6,0.3,0.1))}

#pTrans
p_Trans_fct.B <- function(t, max.time){
  dnorm(t, mean=max.time, sd=2)*5
}

max.time_fct <- function(x){rnorm(x,mean = 5,sd=1)}

param_pTrans.B = list(max.time=max.time_fct)

# Starting the simulation ------------------------------------

set.seed(60)
SimulationDual <- nosoiSim(type="dual", popStructure="discrete",
                           length.sim=300, 
                           max.infected.A=100,
                           max.infected.B=200,
                           init.individuals.A=1,
                           init.individuals.B=0,
                           init.structure.A="A",
                           init.structure.B=NA,
                           structure.matrix.A=transition.matrix,
                           structure.matrix.B=transition.matrix,
                           
                           pExit.A = p_Exit_fct,
                           param.pExit.A=NA,
                           timeDep.pExit.A=FALSE,
                           diff.pExit.A=TRUE,
                           
                           pMove.A = p_Move_fct,
                           param.pMove.A=NA,
                           timeDep.pMove.A=FALSE,
                           diff.pMove.A=FALSE,
                           
                           nContact.A=n_contact_fct,
                           param.nContact.A=NA,
                           timeDep.nContact.A=FALSE,
                           diff.nContact.A=FALSE,
                           
                           pTrans.A = proba,
                           param.pTrans.A = list(p_max=p_max_fct,t_incub=t_incub_fct),
                           timeDep.pTrans.A=FALSE,
                           diff.pTrans.A=FALSE,
                           prefix.host.A="H",
                           
                           pExit.B = p_Exit_fct.B,
                           param.pExit.B=NA,
                           timeDep.pExit.B=TRUE,
                           diff.pExit.B=FALSE,
                           
                           pMove.B = p_Move_fct.B,
                           param.pMove.B=NA,
                           timeDep.pMove.B=FALSE,
                           diff.pMove.B=FALSE,
                           
                           nContact.B=n_contact_fct.B,
                           param.nContact.B=NA,
                           timeDep.nContact.B=FALSE,
                           diff.nContact.B=FALSE,
                           
                           pTrans.B = p_Trans_fct.B,
                           param.pTrans.B = param_pTrans.B,
                           timeDep.pTrans.B=FALSE,
                           diff.pTrans.B=FALSE,
                           prefix.host.B="V",
                           
                           print.progress=FALSE)