epidemic-model-tauseef_working_on.Rmd

---
title: "Modeling psychological impacts on epidemic spread"
output:
  word_document: default
  pdf_document: default
  html_document:
    df_print: paged
---
 
 
# Definitions 
 
 We will try to incorporate opinion dynamics, network modeling, and a bit of the game theory to understand epidemic spread.
 
 
 ## The agent
 
We will use a simplified task network model to represent the biological progression of the disease.
 First, let's suppose that an agent has two states: its psychological state and its biological state.   psychological state might be 'practicing distancing', 'believes conspiracy theory', 'in quarantine', and we can explore these later.  Let's just consider everyone is in a generic 'informed' state [1]. 
 
 The biological state has a few specific cases:
 

1. Unexposed
2. Asymptomatic but infected/contagious
3. Symptomatic and contagious
4. Symptomatic and not contagious
5. Post-COVID Immune
6. Naturally immune (will not contract)
7. Death


We could identify several others, but this might be fine.  So, we can define the agent according to just two numbers.  We might assume that initially, most people are in bio-state 1,  but some would be in state 6 already, which is essentially the same as state 5.

We will just define these as a named list.
```{r}
library(ggplot2)
library(dplyr)
library(sna)
library(knitr)
library(reshape2)
library(igraph)
library(GGally)
library(rlist)
#library(animation)

STATES <<- 8

STATENAMES <<-  c("Unexposed",
                  "Asymptomatic & contagious",
"Symptomatic and contagious",
"Symptomatic and not contagious",
"Post-COVID immune",
"Naturally immune",
"Death",
"Left")


STATELABELS <<-  c("Unexposed","Asymptomatic\n & contagious",
"Symptomatic \n& contagious",
"Symptomatic \n& not contagious",
"Post-COVID immune",
"Naturally immune",
"Death",
"Left network")


```

## Timecourse and Biological model

We need a way to transition the biological state in a reasonable way, much like a task network model that keeps track of time of sub-events.  This should consider ONLY the natural progression of the disease. We will model timecourse on a timecourse of 1-day units.


An easy way to do this is with ballistic events.  That is, we can keep track of the next transition point for any state, if it is programmed at the beginning.  This is similar to how I implemented this in the task network model, but a bit simpler because we only have to worry about changes each day.




```{r}


makeAgent <- function(psychstate,biostate,age=30,traveler=0,is_infected = 0)
{
  
  return (list(psychstate=psychstate,
               biostate=biostate,
               age=age,
               nextbiostate=NA,
               biostatecountdown=NA,traveler=traveler,is_infected = is_infected))
}
print(makeAgent(1,2))


```

This way, we can model the time it takes between becoming infected and showing symptoms.
```{r}
agent <- makeAgent(psychstate=1, biostate=2)

##infect the agent
agent$biostate <- 1   ##infected but asymptomatic
agent$nextbiostate <- 2  ##infected and symptomitac/contageous
agent$biostatecountdown <- ceiling(abs(rnorm(1,mean=8,sd=1)))
print(agent)
```




## Create transition matrix.

If on each timestep we could move between states probabilistically, this would be a markov process. We could represent this by a simple transition network. If we represented the possible transitions in a network, that would simplify the logic above.  But we need more like a task network because the time taken to transition matters.   To keep it simple,   Let's make all timing distributions uniform with a min and max parameter for time in each state.  We can program two pathways through the stages, with a couple branch points (death vs recovery; the possibility of recovering after acquiring with no symptoms).  The progression of the disease is completely specified by the data, and a generic update function will automatically progress the agent each day.


```{r}

# * 1. Unexposed
# * 2. Asymptomatic but infected/contagious
# * 3. Symptomatic and contagious
# * 4. Symptomatic and not contagious
# * 5. Post-COVID Immune
# * 6. Naturally immune (will not contract)
# * 7. Death
# * 8. Death


bioTransition <- matrix(0,STATES,STATES)
bioMin <- matrix(1,STATES)      #state time minimum
bioMax <- matrix(1,STATES)      #state time maximum



bioMin[2] <- 3            #infected but asymptomatic for 3 to 10 days
bioMax[2] <- 10         
bioTransition[2,3] <- .1  #transition to infected with symptoms
bioTransition[2,5] <- .90  #transition to no longer contagious/cured


bioMin[3] <-    3             #symptoms + contagion
bioMax[3] <- 8                #symptoms + contagion max
bioTransition[3,4] <-  .95    #transitioon to no longer contagious
bioTransition[3,7] <-  .05    #transitioon to death state 


bioMin[4] <- 1          #symptoms bot no longer contagiious
bioMax[4] <- 7
bioTransition[4,5] <- 1  #Transition to 'immune' cured state.




setAgentState<- function(agent, biostate)
{
  agent$biostate <- biostate
  if(sum(bioTransition[biostate,])>0) # this state transitions to something else.
  {
    ##which state do we go to?
     agent$biostatecountdown <- sample(x=seq(bioMin[biostate],bioMax[biostate]),1) #how long will we state in this state?
     agent$nextbiostate <- sample(1:STATES, prob=bioTransition[agent$biostate,],size=1)
     
  } else{
   agent$biostatecountdown <- NA
   agent$nextbiostate <- NA   ##just so we can tell if the agent is finished.
  }
  return(agent) 
}

transitionAgent<- function(agent)
{
   return(setAgentState(agent,agent$nextbiostate))
}

updateAgent<- function(agent)
{
  if(!is.na(agent$biostatecountdown))
  {
   agent$biostatecountdown <- agent$biostatecountdown -1
    if(agent$biostatecountdown <=0)  ##new state
    {
       agent <- transitionAgent(agent)

    }
  }
   return(agent)
}

setAgentStatusInfected<- function(agent,status=1)
{
   agent$is_infected <- status
   return(agent)
}


```

```{r, fig.width=8,fig.height=8}

par(xpd=NA,mar=c(5,5,8,5))
gplot(bioTransition,label=STATELABELS)
```






## Social model

This is not quite right, because we'd really expect interactions to be clumpy--you interact with the same people every day, but this model is like the opinion dynamics model in that anyone can interact with anyone.  Let's make a network using preferential attachment to represent that social world.  We can make a couple of them and add them together to represent how we each have different types of relationships and interactions--you may be the president of your drama club and so are a central member, but you are also in a pottery class where you don't talk to very many people (although your teacher does).

```{r}



makeNetwork<- function(numAgents=5,numsets=3,steps=1,power=1)
{
  ord <- sample(numAgents)
  
  #print("____________________ord")
  #print(ord)
  
 
  #old_agent_order <<- ord 
 
  tmp<-as_adjacency_matrix(sample_pa(numAgents,power=power),sparse=F)

  tmp <- tmp + t(tmp)

  tmp <- tmp[ord,ord]

  if(numsets>1)
  {
    for(i in 2:numsets)
  {
    ord <- sample(numAgents)
     
    sn2 <-as_adjacency_matrix(sample_pa(numAgents,power=power),sparse=F)[ord,ord]
    tmp <- tmp + sn2 + t(sn2)
    
   }
  }
  if(steps>1)
  {
   for(i in 2:steps)
   {
      tmp <- tmp + tmp %*% tmp  ## iterate to neighbors
   }
  }
  (tmp>0)+0
   

}

```

```{r}



    
    a <- sample(x, 3)
    print(a)
# innetwork <- sapply(.GlobalEnv$pool,function(y){!(y$biostate==7||y$biostate==8)})

```
```{r}
addAgentToNetwork<- function(numAgents=5,prev_network = socialnetwork,connection_number=1,self_connected=1,pool)
{
  
   if(connection_number < 1){
     
     connection_number <- 1
   }

   number_of_rows <- nrow(prev_network)
   
   starting_point <- number_of_rows+1
   end_point <- number_of_rows + numAgents

  
    for(i in starting_point:end_point) {


      prev_network <- rbind(prev_network, 0)
      prev_network <- cbind(prev_network, 0)

    
    
  
    }  

   
  if(number_of_rows>= connection_number ){
  
    #x <- 1:number_of_rows
    
    u <- 1:number_of_rows

    for (val in u) {
      
      if((pool[[val]]$biostate != 8) ){

        if((pool[[val]]$biostate != 7)){
          x<-u
        }
      }
     
    }    
    
    
     a <- sample(x, connection_number)
     # print("SSSSSSSSSSSSSSSSSSSS")
     # print(a)
     # 
   
    
    # repeat{
    #   
    #   if((pool[[a]]$biostate != 8) ){
    #     
    #     if((pool[[a]]$biostate != 7)){
    #       break
    #     }
    #     
    #   }
    # }
   # print(a)
  
    
     for(i in starting_point:end_point) {
          if(self_connected==1){
            
             prev_network[i,i] <- 1
          }
      
        current_new_agent <-i
        for (i in seq_along(a)) {
          
          #print(a[i])
          prev_network[a[i],current_new_agent] <- 1
          prev_network[current_new_agent,a[i]] <- 1
        }
   
      
      
    }     
    
    
    
  }
   else{
     
     #x <- 1:number_of_rows
     # print("SSSSSSSSSggfgfgSSSSSSSSSSS")
     # print(a)
     #     
    u <- 1:number_of_rows

    for (val in u) {
      
      if((pool[[val]]$biostate != 8) ){

        if((pool[[val]]$biostate != 7)){
          x<-u
        }
      }
     
    }      
     
     a <- sample(x, 1)
    
    
  
    
   # print(a)
  
    
     for(i in starting_point:end_point) {
          if(self_connected==1){
            
             prev_network[i,i] <- 1
          }
      
        current_new_agent <-i
        for (i in seq_along(a)) {
          
          #print(a[i])
          prev_network[a[i],current_new_agent] <- 1
          prev_network[current_new_agent,a[i]] <- 1
        }
   
      
      
    }    
     
     
     
   }
 
   

 
  return(prev_network)
    
  
}


```

```{r}

socialnetwork <- makeNetwork(5,numsets=1,power=.59,steps=1)
print(socialnetwork)
#socialnetwork <- addAgentToNetwork(5,numsets=3,power=.59,steps=1)
```
```{r}
 # socialnetwork <- addAgentToNetwork(numAgents = 1,prev_network = socialnetwork,
 #                                    connection_number = 1,self_connected=0)
 # print(socialnetwork)


#cc <- gplot(socialnetwork,edge.col="black")
```





```{r,fig.width=8,fig.height=8}
mygplot <- function(coord, network,states,main="")
{
  if(is.null(coord))
  {
    coord  <- gplot.layout.fruchtermanreingold(network,layout.par=list(niter=500))
  }
  
  newmin <- mean(coord[,2]) - (-min(coord[,2]) + mean(coord[,2])) * 1.4
   palette=c("white","yellow","red","green","darkgreen","blue","black")
   plot(coord,col="black",bty="n",pch=16,cex=2.7,xaxt="n",yaxt="n",main=main,xlab="",ylab="",axes=F,
        ylim=c(newmin,max(coord[,2])),type="n")
   
  for(i in 1:nrow(network))
     segments(coord[i,1],
            coord[i,2],
            coord[network[i,]==1,1],
            coord[network[i,]==1,2],col="grey40")

       points(coord,pch=16,cex=2.3,col= palette[states])
              
       points(coord,pch=1,cex=2.3,col="black")
      legend(mean(coord[,1]),min(coord[,2]),bty='n',y.intersp=.7,cex=.8,
              STATENAMES, pch=16,col=palette)
              
      return (coord)
  }

#cc <-  mygplot(coord=NULL,socialnetwork,rep(1,nrow(socialnetwork)),main="Initial state")
```
#########
###########
######
##########

##################################################Beginning###############################################################

A key step in contingency planning for a possible pandemic in a community is barring entry of travelers in the community. The potential for infectious diseases to spread in a community rapidly through well-connected infected travelers in the community is very high. Past researches looked into the effect of the international travel ban on reducing the spread of disease in countries. Partial restriction of travel helps in slowing the spread if the global number of cases of infection is low (Hollingsworth et al., 2006). According to Epstein et al. (2007), international air travel restrictions may provide a small but important delay in the spread of a pandemic, especially if other disease control measures are implemented during a timeframe. Using the mathematical model, Cooper et al. (2006) concluded that restrictions on air travel would achieve very little. This is probably because a virus may transmit from one person to another very quickly and affects many people. Once a major outbreak is underway, banning flights from affected cities would be effective at significantly delaying spread only if almost all travel between cities could be stopped almost as soon as an outbreak was detected in each city.

One of the main reasons against oversea travel ban is because it has a big negative economic impact in a country if we stop adults traveling oversea (Lam et al., 2011). But for small and moderate length communities whose economic activities circulated within the community, a travel restriction may be helpful during a pandemic. If the health care systems for these communities are not prepared to deal with a pandemic, it is better to enforce some sort of travel ban so that the health care system does not overwhelm with patients. This will help to keep a check on the number of patients. In the traditional SIR model (Smith & Moore, 2004), it is not defined how adding travelers to the scenario will affect the spread of disease. So, we created an agent-based social model for a community where we will be adding travelers to the local populations to see how a disease spreads in the whole community. We will be looking into the worst-case scenario in the model, where travelers will enter the community but nobody leaves within a timeframe from the community. We will vary different variables in the simulations to see the disease progress in the community within a specific timeframe.

References

Cooper, B. S., Pitman, R. J., Edmunds, W. J., & Gay, N. J. (2006). Delaying the international spread of pandemic influenza. PLoS Medicine, 3(6).
Epstein, J. M., Goedecke, D. M., Yu, F., Morris, R. J., Wagener, D. K., & Bobashev, G. V. (2007). Controlling pandemic flu: The value of international air travel restrictions. PloS One, 2(5).
Hollingsworth, T. D., Ferguson, N. M., & Anderson, R. M. (2006). Will travel restrictions control the international spread of pandemic influenza? Nature Medicine, 12(5), 497–499.
Lam, E. H., Cowling, B. J., Cook, A. R., Wong, J. Y., Lau, M. S., & Nishiura, H. (2011). The feasibility of age-specific travel restrictions during influenza pandemics. Theoretical Biology and Medical Modelling, 8(1), 44.
Smith, D., & Moore, L. (2004). The SIR Model for Spread of Disease: The Differential Equation Model. Loci.(Originally Convergence.) Https://Www. Maa. Org/Press/Periodicals/Loci/Joma/the-Sir-Model-for-Spread-of-Disease-the-Differential-Equation-Model.


##Travel ban lift and allowing traveler to the community
###social connection of the travelers




```{r,fig.width=8,fig.height=8}
mygplot_updated <- function(coord, network,states,main="")
{
  if(is.null(coord))
  {
    coord  <- gplot.layout.fruchtermanreingold(network,layout.par=list(niter=500))
      #gplot.layout.fruchtermanreingold(network,layout.par=list(niter=500))
      #ggnet2(network, mode = "fruchtermanreingold", layout.par=list(niter=500))
      #gplot.layout.fruchtermanreingold(network,layout.par=list(niter=500))
  }
  
   newmin <- mean(coord[,2]) - (-min(coord[,2]) + mean(coord[,2])) * 1.4
   palette=c("white","yellow","red","green","darkgreen","blue","black","grey")
   plot(coord,col="black",bty="n",pch=16,cex=2.7,xaxt="n",yaxt="n",main=main,xlab="",ylab="",axes=F,
        ylim=c(newmin,max(coord[,2])),type="n")
   
   
   travelers <- sapply(pool,function(x){x$traveler})
   
  for(i in 1:nrow(network))
  {
    if(sum(network[i,])>0)
    {
     segments(coord[i,1],
            coord[i,2],
            coord[network[i,]==1,1],
            coord[network[i,]==1,2],col="grey40")
    } 
  }
       points(coord,pch=c(16,15)[travelers+1],cex=2.3,col= palette[states])
       points(coord,pch=c(1,0)[travelers+1],cex=2.3,col="black")
       
      legend(mean(coord[,1]),min(coord[,2]),bty='n',y.intersp=.7,cex=.8,
              STATENAMES, pch=16,col=palette)
      
      if(sum(travelers)>0)
      {
        #print("plotting rectangle")
        rect(min(coord[travelers,1]),
             min(coord[travelers,2]),
             max(coord[travelers,1]),
             max(coord[travelers,2]))
          
      }
  
              
      return (coord)
}
#cc <-  mygplot_updated(coord=NULL,socialnetwork,rep(1,nrow(socialnetwork)),main="Initial state")
```


##############################
```{r fig.width=10,fig.height=10}
travel_ban_lift <- function( enable_outside_travel_restriction_day,second_wave_day,
                             connection_number = 1, will_affect_traveler = 0, numDays = 50,
                              people_leave_status = 1,
                             number_of_people_leave = 4,
                             naturalImmunity = .01,numInteractions_value = 14,
                              contagionProb = .05,####5% chance
                             sampleFromNetwork = .98,traveler_percentage_mean = 1,
                             traveler_percentage_ceiling = 2, traveler_increase = 0,
                             boundary_of_traveler_change = 2, traveler_interval_day = 5, 
                             enable_second_wave = 1,
                             enable_outside_travel_restriction =0,need_to_update_graph_cord = 0 ){ 



 
  socialnetwork <- original_network_of_agents
  numAgents <- numAgents_initial
  people_left <-0
  

  
  plotNetwork <-T
  
  # If i create the social network here, works ok
  #socialnetwork <-makeNetwork(numAgents,numsets=1,power=.5,steps=2)
  #print(socialnetwork)
  
  #pool <<-pool_initial
  
  .GlobalEnv$pool <<- pool_initial
  
  numInteractions <-  rep(numInteractions_value,numDays) ##interaction during no travel bans
  numInteractions[enable_outside_travel_restriction_day:second_wave_day] <- as.integer(numInteractions_value*(1/3))  ##interaction during travel restriction
  
  
  local_population_infection_number_local_variable <<- local_population_infection_number
  traveler_infection_number_local_variable <- 0

  if(plotNetwork)
  {
    cc <-mygplot_updated(coord=NULL,socialnetwork,rep(1,nrow(socialnetwork)),main="Initial state")
    
    xrange <-   range(cc[,1])
    xval <- xrange[2] + (xrange[2]-xrange[1])/3
    yval <- mean(cc[,2])
    sd <- (xrange[2]-xrange[1]) /10
  }  

  
  

  disthistory <- matrix(NA,ncol=STATES,nrow=numDays)
  
 # pool <<-pool_initial

  
  for(day in 1:numDays)
  {
      
      
       ###################agent leaving Start
    
      #innetwork <- sapply(.GlobalEnv$pool,function(y){!(y$biostate==7||y$biostate==8)})
      if(!enable_outside_travel_restriction){
        if(people_leave_status){
          ####Hard coding day for people leave
           #.GlobalEnv$pool <- list.remove(.GlobalEnv$pool,  c(1))
          
          if(day == (second_wave_day+7)){
                  
                  number_of_people_leave <- as.integer((1/100)*numAgents) ##1% will leave the system
                  
                      
                    u <- 1:numAgents
                
                    for (val in u) {
                      
                      if((.GlobalEnv$pool[[val]]$biostate != 8) ){
                
                        if((.GlobalEnv$pool[[val]]$biostate != 7)){
                          x<-u
                          
                        }
                      }
                     
                    }    
                    
                    a <- sample(x, number_of_people_leave)
                    
                            
                    for (i in seq_along(a)) {
                      
                      #print(a[i])
                      x <- a[i]
                      .GlobalEnv$pool[[x]] <- setAgentState( .GlobalEnv$pool[[x]],8)
          
                      socialnetwork[x,] <- 0
                      socialnetwork[,x] <-0
                    }                
                           
                  
                  
                 
          }          
          if(day == (second_wave_day+14)){
                  
                  number_of_people_leave <- as.integer((1/100)*numAgents) ##1% will leave the system
                  
                      
                    u <- 1:numAgents
                
                    for (val in u) {
                      
                      if((.GlobalEnv$pool[[val]]$biostate != 8) ){
                
                        if((.GlobalEnv$pool[[val]]$biostate != 7)){
                          x<-u
                          
                        }
                      }
                     
                    }    
                    
                    a <- sample(x, number_of_people_leave)
                    
                            
                    for (i in seq_along(a)) {
                      
                      #print(a[i])
                      x <- a[i]
                      .GlobalEnv$pool[[x]] <- setAgentState( .GlobalEnv$pool[[x]],8)
          
                      socialnetwork[x,] <- 0
                      socialnetwork[,x] <-0
                    }                
                           
                  
                  
                 
          }              

          # if(day == (second_wave_day+5)){
          #   
          #   d <- 1:numAgents
          #   number_traveler_found
          #   repeat{
          #     statements...
          #     if(condition){
          #       break
          #     }
          #   }            
          #   
          #   
          #   
          #       x <- 18
          #       .GlobalEnv$pool[[x]] <- setAgentState( .GlobalEnv$pool[[x]],8)
          # 
          #       socialnetwork[x,] <- 0
          #       socialnetwork[,x]<-0
          # }



        }


      }

      ###################agent leaving ENd
    


          innetwork <- sapply(.GlobalEnv$pool,function(y){!(y$biostate==7||y$biostate==8)})
  
  ##who are you going to talk to today.
      sneezers <- rep((1:numAgents)[innetwork],each=numInteractions[day])
      sneezedons <- rep(NA,length(sneezers))
      
      for(i in 1:length(sneezers))
        {
        if(runif(1)<(1-sampleFromNetwork))
        {
    
          sneezedons[i] <- sample((1:numAgents)[innetwork],1)
        }else{
          
          probs <- socialnetwork[sneezers[i],]* innetwork
          sneezedons[i] <- sample(1:numAgents,prob=probs,1)
          
        }
       }
      
        
    for(i in 1:length(sneezers))
    {
      
      agent1 <- .GlobalEnv$pool[[ sneezers[i] ]]
      agent2 <- .GlobalEnv$pool[[ sneezedons[i] ]]
    
      
      ##this constitutes the rules of infection.
      if((agent1$biostate==2 || agent1$biostate==3 ) & agent2$biostate==1 & runif(1)<contagionProb)
      {
           .GlobalEnv$pool[[ sneezedons[i] ]] <- setAgentState(agent2,2)##infect!
      }
        
    }
      ##increment each agent 1-day.
      for(i in 1:numAgents)
      {
         .GlobalEnv$pool[[i]] <- updateAgent(.GlobalEnv$pool[[i]])
         
         ###keep track number of local population infection
         
           
           if(.GlobalEnv$pool[[i]]["traveler"] == 0 ){
             if(.GlobalEnv$pool[[i]]["biostate"] == 2){
               if(.GlobalEnv$pool[[i]]["is_infected"] == 0){         
                 local_population_infection_number_local_variable <-local_population_infection_number_local_variable+1
                .GlobalEnv$pool[[i]] <- setAgentStatusInfected( .GlobalEnv$pool[[i]],1)
               }
               
             }
            
           
         }         
           if(.GlobalEnv$pool[[i]]["traveler"] == 1 ){
             if(.GlobalEnv$pool[[i]]["biostate"] == 2){
               if(.GlobalEnv$pool[[i]]["is_infected"] == 0){         
                 traveler_infection_number_local_variable <-traveler_infection_number_local_variable+1
                .GlobalEnv$pool[[i]] <- setAgentStatusInfected( .GlobalEnv$pool[[i]],1)
               }
               
             }
            
           
         }            
         
         
      }
    
      states <- sapply(.GlobalEnv$pool,FUN=function(x){x$biostate})
      distrib <- table(factor(states,levels=1:STATES))
      disthistory[day,] <- distrib
     # print(disthistory[day,])
      
      
     
  
      
       #####################################
      ############################
      ##################################
      ################
      ################################   
      ####################Adding New agents as travelers in the network
      
      #traveler_interval_day <-sample(c(3,2,5),1)
      
     # traveler_increase <- sample(c(0,1),1)
     
      
    
      if(traveler_increase == 0){
        if(day == second_wave_day){
          
          if(enable_second_wave == 1){
            traveler_percentage_mean <- as.integer(traveler_percentage_ceiling/1.5)
            
             
          }
          
          
        }
      }
      ##travel restriction on
  
      
      
      ###if second wave already happned due to premature removal of travel ban. traveler will come in future days
      if((enable_second_wave == 1) && (day == second_wave_day)){
        
        enable_outside_travel_restriction <-0
        
      }
      else if((enable_second_wave == 1) && (day > second_wave_day)){
        
        enable_outside_travel_restriction <-0 ###if we do not want to close down the travel after second wave
      }
      else if(enable_outside_travel_restriction_day == 0){
        
         enable_outside_travel_restriction <-0
         
      }else if(enable_outside_travel_restriction_day == day){
         enable_outside_travel_restriction <-1
        
      }else if(enable_outside_travel_restriction == 1){
        
        enable_outside_travel_restriction <-1
        
      }
      else{
        
        enable_outside_travel_restriction <-0
        
      }
      
      

      
  
      ####
      if(((!(day %% traveler_interval_day)) ||
           (day==second_wave_day)) && 
           (enable_outside_travel_restriction==0)){
        
            if(traveler_percentage_mean < 1){

              traveler_percentage_mean =  as.integer(traveler_percentage_ceiling/2)
            }
            else if(traveler_percentage_mean>traveler_percentage_ceiling){
              
              traveler_percentage_mean = traveler_percentage_ceiling
            }
     
        
          if(traveler_percentage_mean > 0){
            number_of_traveler <-as.integer((traveler_percentage_mean/100) * numAgents_initial)
            current_number_traveler_remaining <- number_of_traveler
            if(traveler_increase < 1){

              traveler_percentage_mean = traveler_percentage_mean - boundary_of_traveler_change


            }else{

              traveler_percentage_mean = traveler_percentage_mean + boundary_of_traveler_change

            

            }
            ###enabling travel restruition after second wave
            if(enable_second_wave == 1){
              if((day==second_wave_day) && (enable_outside_travel_restriction_day != 0)){
                
                enable_outside_travel_restriction<-1
              }
            }
            
            ###for evenly distributing travel to network I run a do while loop
            #print("number_of_traveler")
            #print(number_of_traveler)
            #print("-----------------------start-------------------")
            repeat{
              current_number_traveler <- sample(c(1,2,3,4,5),1)
                
                #
  
              
              if(current_number_traveler>current_number_traveler_remaining){
                
                current_number_traveler <- current_number_traveler_remaining
              }
              
             # print(current_number_traveler)
  
              
              if(current_number_traveler > 0)
              {
                starting_numagents <- numAgents + 1
                ending_numagent <- numAgents+current_number_traveler
                
                ####connection_number indicates new agent connected how many number of previous agents, currently it is a number between 1-3
                
                socialnetwork <- addAgentToNetwork(numAgents = current_number_traveler,
                                                   prev_network = socialnetwork,
                                                   connection_number = connection_number,
                                                   self_connected=1,.GlobalEnv$pool)
               # print(socialnetwork)
                
                newxy <- cbind(rnorm(current_number_traveler,mean=xval,sd=sd),
                               rnorm(current_number_traveler, mean=yval,sd=sd))
                
                cc <- rbind(cc,newxy)
                for(i in starting_numagents:ending_numagent)
                  {
                    ###printing to see how are travelers
                    # print("adding")
                     .GlobalEnv$pool[[i]] <<- makeAgent(psychstate=1,
                                            biostate=sample(c(1,6),
                                            p=c(1-naturalImmunity, naturalImmunity),1),
                                            traveler = 1)
                   # print("added")
                    ################to make randomize infect -> traveler
                     
                     will_be_infected <- sample(c(0,1), 1)
                     
                    
                     
                    # print(paste("infecting",i))
                     #####Tell shane about if setAgentState is 2, then number of local remains ok, but if state is 1 then it changes
                     if(will_affect_traveler == 1){
                       #print("poit a")
                       tmp <-  setAgentState(.GlobalEnv$pool[[i]],2) ##infect this person
                       #print("point b")
                       .GlobalEnv$pool[[i]] <<-tmp
                       
                     }else{
          
                      # print("point b1")
                      # tmp <- setAgentState(.GlobalEnv$pool[[i]],1) ##do not infect this person
                      # print("point b2")
                      # .GlobalEnv$pool[[i]] <<- tmp
                      # print("point b3")
                     }
                    # print("infected")
                     
                     # 
                     
                     
                 
                     
                     
                }
                numAgents <- ending_numagent
                need_to_update_graph_cord <- 1
              }
              current_number_traveler_remaining <- current_number_traveler_remaining-current_number_traveler
              if(current_number_traveler_remaining < 1){
                 #print("-----------------------end-------------------")
                break
              }
              
            }
        }
      }
      
     if(plotNetwork)
     {
        ## print(1)
       
       ##show this to SHane
        # cc <-mygplot(coord=NULL,socialnetwork,rep(1,nrow(socialnetwork)),main="Initial state")
        # mygplot(coord=NULL,socialnetwork,states,main=paste("Day",day))
       
       if(need_to_update_graph_cord){
        # print("UPDATE NEEDED")
        
         mygplot_updated(cc,socialnetwork,states,main=paste("Day",day, " ", "Traveler Added", number_of_traveler))
        
         
       }else{
         
         mygplot_updated(cc,socialnetwork,states,main=paste("Day",day))
       }
        need_to_update_graph_cord <- 0
       
        
         
         
         
  
     }    
      
      
  ###vary week for travel restriction
     
    
      
         
  }

  total_local <-0
  affected_locals <- 0
  affected_travelers <- 0
  local_death <- 0
  total_traveler <-0
  traveler_death <- 0
  ok_travler <- 0
  ok_local<- 0
  #print(disthistory)
     
   for(i in 1:numAgents)
   {
     
     if(pool[[i]]["traveler"] == 0){
       
       #print(i)
       
       total_local <-total_local+1
       if(pool[[i]]["biostate"] == 2 || pool[[i]]["biostate"] == 3 || pool[[i]]["biostate"] == 4)
         {
           affected_locals <-affected_locals + 1
       }
       else if(pool[[i]]["biostate"] == 7){
         
         local_death <-local_death+1
       }
       else{
         
         ok_local <- ok_local+1
       }
       
     }
     else{
       total_traveler <- total_traveler+1
       if(pool[[i]]["biostate"] == 2 || pool[[i]]["biostate"] == 3 || pool[[i]]["biostate"] == 4)
         {
           affected_travelers <-affected_travelers + 1
       }
       else if(pool[[i]]["biostate"] == 7){
         
         traveler_death <-traveler_death+1
       }  
       else{
         
         ok_travler <- ok_travler+1
       }       
       
       
     }
   }
  
    #.GlobalEnv$pool <- list.remove(.GlobalEnv$pool, 1)
    #print(.GlobalEnv$pool)
  print("TOTAL")
  print(numAgents)
  print("Total Local")
  print(total_local)
  # print("Local Recovered/Not affected")
  # print(ok_local)
  # print("Local Affected")
  # print(affected_locals)
  print("Local Death")
  print(local_death)
  print("-----")
  print("Total Traveler")
  print(total_traveler)  
  
  # print("Traveler Recovered/Not affected")
  # print(ok_travler)
  # print("Traveler Affected")
  # print(affected_travelers)
  print("Traveler Death")
  print(traveler_death)
  
  
  print("Local infected at somepoint")
  print(local_population_infection_number_local_variable)
  print("Traveler infected at somepoint")
  print(traveler_infection_number_local_variable) 
  print("Total infected at somepoint")
  print(traveler_infection_number_local_variable + local_population_infection_number_local_variable) 
  
 # print(traveler_percentage_mean)
  

  


  
  return(disthistory)
  
  


 
}


```

We will run the model for a community of 700 local agents, for 100 days

The simulations will be divided into major parts. Interactions among agents will be one-third of  normal interactions during travel ban

In the first part, we will see how travlers affect the local populations when no from the local populaton if infeceted.




###First Part of the simulations (No locals are affected with a strain from beginning)


```{r,fig.width=15,fig.height=15}
  numAgents_initial <<- 700
  naturalImmunity <- .01
  socialnetwork <-makeNetwork(numAgents_initial,numsets=1,power=.5,steps=2)
  local_population_infection_number <-0
  
  infect_local <- 0
  
  original_network_of_agents <<- socialnetwork
  
  .GlobalEnv$pool <- list()
  for(i in 1:numAgents_initial)
  {
     .GlobalEnv$pool[[i]] <- makeAgent(psychstate=1,
                           biostate=sample(c(1,6),
                            p=c(1-naturalImmunity, naturalImmunity),1))
  }
  
  ##infect patient 0
  numInfected <- 3
  for(i in sample(numAgents_initial,numInfected))
  {
    if(infect_local){
        
       .GlobalEnv$pool[[i]] <- setAgentState(pool[[i]],2) ##infect this person
       
       
       if(.GlobalEnv$pool[[i]]["biostate"] == 2){
         if(.GlobalEnv$pool[[i]]["is_infected"] == 0){         
           local_population_infection_number <-local_population_infection_number+1
          .GlobalEnv$pool[[i]] <- setAgentStatusInfected( .GlobalEnv$pool[[i]],1)
         }
         
       }
    }
    
  } 
  
  pool_initial <<-.GlobalEnv$pool
  local_population_infection_number <<- local_population_infection_number
  

  

```




###Travel ban is implemented from Day 7 and removed early (In this case from day 49). Each traveler will have only three connections with local agents

```{r,fig.width=10,fig.height=10}

  numDays = 100
  enable_outside_travel_restriction_day <- 7
  second_wave_day <- as.integer((49/100)*numDays)
    #

  #as.integer((numAgents_initial)*(1/100))
  disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
                  connection_number = 3,will_affect_traveler = 1,numDays = numDays  )

  disthist.df <-as.data.frame(disthistory)
  colnames(disthist.df) <- STATENAMES
  disthist.df$day <- 1:nrow(disthistory)
  histlong <- melt(disthist.df,id.vars="day")

  ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") +
    theme_bw()

  ##make the SIR plot:
  sir <- data.frame(day=disthist.df$day,
                    susceptible = disthist.df$Unexposed,
                    infected = disthist.df[,2]+disthist.df[,3],
                    recovered = rowSums(disthist.df[,4:7]))

  ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line(size=2) + theme_bw()
  ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line(size=2) + theme_bw()

  
```

###Travel ban is implemented from Day 7 and removed early (In this case from day 49). Each traveler will have only one connection with local agents


```{r,fig.width=10,fig.height=10}

  numDays = 100
  enable_outside_travel_restriction_day <- 7
  second_wave_day <- as.integer((49/100)*numDays)

  #as.integer((numAgents_initial)*(1/100))
  disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
                  connection_number = 1,will_affect_traveler = 1,numDays = numDays  )

  disthist.df <-as.data.frame(disthistory)
  colnames(disthist.df) <- STATENAMES
  disthist.df$day <- 1:nrow(disthistory)
  histlong <- melt(disthist.df,id.vars="day")

  ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") +
    theme_bw()

  ##make the SIR plot:
  sir <- data.frame(day=disthist.df$day,
                    susceptible = disthist.df$Unexposed,
                    infected = disthist.df[,2]+disthist.df[,3],
                    recovered = rowSums(disthist.df[,4:7]))

  ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line(size=2) + theme_bw()
  ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line(size=2) + theme_bw()
  
  
```


###Travel ban is implemented from Day 7 and removed later (In this case from day 80). Each traveler will have only three connection with local agents


```{r,fig.width=10,fig.height=10}

  # numDays = 250
  # enable_outside_travel_restriction_day <- 7
  # second_wave_day <- as.integer((80/100)*numDays)
  # 
  # #as.integer((numAgents_initial)*(1/100))
  # disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
  #                 connection_number = 3,will_affect_traveler = 1,numDays = numDays  )
  #   
  # disthist.df <-as.data.frame(disthistory)
  # colnames(disthist.df) <- STATENAMES
  # disthist.df$day <- 1:nrow(disthistory)
  # histlong <- melt(disthist.df,id.vars="day")
  # 
  # ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
  #   theme_bw()
  # 
  # ##make the SIR plot:
  # sir <- data.frame(day=disthist.df$day,
  #                   susceptible = disthist.df$Unexposed,
  #                   infected = disthist.df[,2]+disthist.df[,3],
  #                   recovered = rowSums(disthist.df[,4:7]))
  # 
  # ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line(size = 2) + theme_bw()
  # ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line(size = 2) + theme_bw()
  
  
```



###Travel ban is implemented from Day 7 and removed later (In this case from day 80). Each traveler will have only one connection with local agent


```{r,fig.width=10,fig.height=10}

  # numDays = 250
  # enable_outside_travel_restriction_day <- 7
  # second_wave_day <- as.integer((80/100)*numDays)
  # 
  # #as.integer((numAgents_initial)*(1/100))
  # disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
  #                 connection_number = 1,will_affect_traveler = 1,numDays = numDays  )
  #   
  # disthist.df <-as.data.frame(disthistory)
  # colnames(disthist.df) <- STATENAMES
  # disthist.df$day <- 1:nrow(disthistory)
  # histlong <- melt(disthist.df,id.vars="day")
  # 
  # ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
  #   theme_bw()
  # 
  # ##make the SIR plot:
  # sir <- data.frame(day=disthist.df$day,
  #                   susceptible = disthist.df$Unexposed,
  #                   infected = disthist.df[,2]+disthist.df[,3],
  #                   recovered = rowSums(disthist.df[,4:7]))
  # 
  # ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line() + theme_bw()
  # ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line() + theme_bw()
  
  
```






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




###Second Part of the simulations (some locals are affected with a strain from beginning)


```{r,fig.width=15,fig.height=15}
  numAgents_initial <<- 700
  naturalImmunity <- .01
  socialnetwork <-makeNetwork(numAgents_initial,numsets=1,power=.5,steps=2)
  local_population_infection_number <-0
  
  infect_local <- 1
  
  original_network_of_agents <<- socialnetwork
  
  .GlobalEnv$pool <- list()
  for(i in 1:numAgents_initial)
  {
     .GlobalEnv$pool[[i]] <- makeAgent(psychstate=1,
                           biostate=sample(c(1,6),
                            p=c(1-naturalImmunity, naturalImmunity),1))
  }
  
  ##infect patient 0
  numInfected <- 3
  for(i in sample(numAgents_initial,numInfected))
  {
    if(infect_local){
        
       .GlobalEnv$pool[[i]] <- setAgentState(pool[[i]],2) ##infect this person
       
       
       if(.GlobalEnv$pool[[i]]["biostate"] == 2){
         if(.GlobalEnv$pool[[i]]["is_infected"] == 0){         
           local_population_infection_number <-local_population_infection_number+1
          .GlobalEnv$pool[[i]] <- setAgentStatusInfected( .GlobalEnv$pool[[i]],1)
         }
         
       }
    }
    
  } 
  
  pool_initial <<-.GlobalEnv$pool
  local_population_infection_number <<- local_population_infection_number
  

  

```


## IF no traveler is allowed and some locals are infected



```{r,fig.width=10,fig.height=10}

  # numDays = 100
  # enable_outside_travel_restriction_day <- 1
  # second_wave_day <- 101
  # 
  # #as.integer((numAgents_initial)*(1/100))
  # disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
  #                 connection_number = 3,will_affect_traveler = 1,numDays = numDays,people_leave_status = 0  )
  #   
  # disthist.df <-as.data.frame(disthistory)
  # colnames(disthist.df) <- STATENAMES
  # disthist.df$day <- 1:nrow(disthistory)
  # histlong <- melt(disthist.df,id.vars="day")
  # 
  # ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
  #   theme_bw()
  # 
  # ##make the SIR plot:
  # sir <- data.frame(day=disthist.df$day,
  #                   susceptible = disthist.df$Unexposed,
  #                   infected = disthist.df[,2]+disthist.df[,3],
  #                   recovered = rowSums(disthist.df[,4:7]))
  # 
  # ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line(size=2) + theme_bw()
  # ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line(size=2) + theme_bw()
  
  
```

###Travel ban is implemented from Day 7 and removed early (In this case from day 49). Each traveler will have only three connections with local agents

```{r,fig.width=10,fig.height=10}

  numDays = 100
  enable_outside_travel_restriction_day <- 7
  second_wave_day <- as.integer((49/100)*numDays)
  
  #as.integer((numAgents_initial)*(1/100))
  disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
                  connection_number = 3,will_affect_traveler = 1,numDays = numDays  )
    
  disthist.df <-as.data.frame(disthistory)
  colnames(disthist.df) <- STATENAMES
  disthist.df$day <- 1:nrow(disthistory)
  histlong <- melt(disthist.df,id.vars="day")
  
  ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
    theme_bw()
  
  ##make the SIR plot:
  sir <- data.frame(day=disthist.df$day,
                    susceptible = disthist.df$Unexposed,
                    infected = disthist.df[,2]+disthist.df[,3],
                    recovered = rowSums(disthist.df[,4:7]))
  
  ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line(size=2) + theme_bw()
  ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line(size=2) + theme_bw()
  
  
```

###Travel ban is implemented from Day 7 and removed early (In this case from day 49). Each traveler will have only one connection with local agents


```{r,fig.width=10,fig.height=10}

  # numDays = 100
  # enable_outside_travel_restriction_day <- 7
  # second_wave_day <- as.integer((49/100)*numDays)
  # 
  # #as.integer((numAgents_initial)*(1/100))
  # disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
  #                 connection_number = 1,will_affect_traveler = 1,numDays = numDays  )
  #   
  # disthist.df <-as.data.frame(disthistory)
  # colnames(disthist.df) <- STATENAMES
  # disthist.df$day <- 1:nrow(disthistory)
  # histlong <- melt(disthist.df,id.vars="day")
  # 
  # ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
  #   theme_bw()
  # 
  # ##make the SIR plot:
  # sir <- data.frame(day=disthist.df$day,
  #                   susceptible = disthist.df$Unexposed,
  #                   infected = disthist.df[,2]+disthist.df[,3],
  #                   recovered = rowSums(disthist.df[,4:7]))
  # 
  # ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line() + theme_bw()
  # ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line() + theme_bw()
  
  
```


###Travel ban is implemented from Day 7 and removed later. Each traveler will have only three connection with local agents


```{r,fig.width=10,fig.height=10}

  numDays = 250
  enable_outside_travel_restriction_day <- 7
  second_wave_day <- as.integer((80/100)*numDays)
  
  #as.integer((numAgents_initial)*(1/100))
  disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
                  connection_number = 3,will_affect_traveler = 1,numDays = numDays  )
    
  disthist.df <-as.data.frame(disthistory)
  colnames(disthist.df) <- STATENAMES
  disthist.df$day <- 1:nrow(disthistory)
  histlong <- melt(disthist.df,id.vars="day")
  
  ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
    theme_bw()
  
  ##make the SIR plot:
  sir <- data.frame(day=disthist.df$day,
                    susceptible = disthist.df$Unexposed,
                    infected = disthist.df[,2]+disthist.df[,3],
                    recovered = rowSums(disthist.df[,4:7]))
  
  ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line(size=2) + theme_bw()
  ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line(size=2) + theme_bw()
  
  
```



###Travel ban is implemented from Day 7 and removed later (In this case from day 80). Each traveler will have only one connection with local agents


```{r,fig.width=10,fig.height=10}

  # numDays = 250
  # enable_outside_travel_restriction_day <- 7
  # second_wave_day <- as.integer((80/100)*numDays)
  # 
  # #as.integer((numAgents_initial)*(1/100))
  # disthistory <- travel_ban_lift(enable_outside_travel_restriction_day,second_wave_day,
  #                 connection_number = 1,will_affect_traveler = 1,numDays = numDays  )
  #   
  # disthist.df <-as.data.frame(disthistory)
  # colnames(disthist.df) <- STATENAMES
  # disthist.df$day <- 1:nrow(disthistory)
  # histlong <- melt(disthist.df,id.vars="day")
  # 
  # ggplot(histlong,aes(x=day,y=value,fill=variable)) + geom_bar(stat="identity",position="stack") + 
  #   theme_bw()
  # 
  # ##make the SIR plot:
  # sir <- data.frame(day=disthist.df$day,
  #                   susceptible = disthist.df$Unexposed,
  #                   infected = disthist.df[,2]+disthist.df[,3],
  #                   recovered = rowSums(disthist.df[,4:7]))
  # 
  # ggplot(melt(sir,id.vars="day"),aes(x=day,group=variable,y=value,color=variable)) + geom_line() + theme_bw()
  # ggplot(histlong,aes(x=day,y=value,group=variable,color=variable)) + geom_line() + theme_bw()
  
  
```

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


```{r}

```













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