COVID19-variable-SD-Testing

;-------------------------------------------------------------------------------------------------------------------
;
;
;
; percentage_SD:      % of agents who keep SD
;
; SD_mean, SD_sigma:  mean value and variance of SD (normal distribution)
;
; movement_type:      "away from mean" o "metropolis T=0"
;
; test_period:        ticks between tests
;
; num_to_test:        number of tests performed each test_period ticks
;
; max_inf_radius:     maximum radius of infection
; transmissibility:   probability of becoming infected when there is a patient less than max_inf_radius
;
;
;-------------------------------------------------------------------------------------------------------------------

extensions [array]     ; to use arrays

globals [acum acum2 sigma n_experiments t s1 s2]

breed [infecteds infected]
breed [susceptibles susceptible]
breed [recovereds recovered]

infecteds-own [ clock tested? ]

turtles-own [ neighbours keep_SD SD_own xcor_old ycor_old ]




;-------------------------------------------------------------------------------------------------------------------
; General setup and setup to restart an experiment
;-------------------------------------------------------------------------------------------------------------------

to setup
  set acum array:from-list n-values t_max [0]   ; fill the array with zeros
  set acum2 array:from-list n-values t_max [0]
  set sigma array:from-list n-values t_max [0]
  set n_experiments 1        ; set number of experiments initially to zero

  setup-experiment

  let LX ( max-pxcor * 2 + 1)
  let LY ( max-pycor * 2 + 1)

  clear-output
  output-print "DENSITY:"
  output-write (number / (LX * LY)) output-print " agents/m2"
  output-write ((LX * LY) / number) output-print " m2/agent"
  output-write sqrt ((LX * LY) / number) output-print " average distance"
end 

; setup for a new experiment

to setup-experiment

  clear-turtles
  clear-ticks
  random-seed new-seed    ; to start each simulation with a different seed
  create-turtles number
  [setxy random-xcor random-ycor
    set size 1
    set keep_SD ((random-float 1) < (percentage_SD / 100.0))  ; is a fixed attribute of each agent
    (ifelse distri_SD = "uniforme"
    [
      set SD_own SD_mean
    ]
    distri_SD = "normal"
    [
      set SD_own random-normal SD_mean SD_sigma
      if (SD_own < 0 )  [set SD_own 0 ]
    ]
    distri_SD = "poisson"
    [
      set SD_own random-poisson SD_mean
    ])

    ifelse (who < number_sick_initial)
    [set breed infecteds
    set clock 0
    set tested? False]
    [set breed susceptibles]
    recolor
  ]
  reset-ticks

  file-close ;cierra el lammps_file
  if file-exists? lammps_file     
      [ file-delete lammps_file ] 
  file-open lammps_file
end 

to recolor
  if (breed = susceptibles) [
    set color green
  ]
  if (breed = infecteds) [
    ifelse (tested? = True)
    [set color white]
    [set color red]
  ]
  if (breed = recovereds) [
    set color blue
  ]
end 


;-------------------------------------------------------------------------------------------------------------------
; move
;-------------------------------------------------------------------------------------------------------------------

to move [dist]
  ifelse(keep_SD) ; that an agent save SD is an attribute that is defined in setup with some probability
  [ ; aqui guarda SD...

    (ifelse
    movement_type = "metropolis T=0"
    [
      ; type of Metropolis movement at temperature 0 (move-count-accept-reject)
      set xcor_old xcor
      set ycor_old ycor
      rt random-float 360
      fd dist
      if (count(other turtles in-radius SD_own) >= 1 )
      [ set xcor xcor_old
        set ycor ycor_old
      ]

    ]
    movement_type = "away from mean"
    [
      ; type of movement runs away from the center of mass of neighbors (counts-turns 180-moves)
      set neighbours other turtles in-radius SD_own
      if (count(neighbours) >= 1) ;Move only if you have neighbours
      [
        facexy (mean [xcor] of neighbours)
               (mean [ycor] of neighbours)
        rt 180  ; Turn away from the mean x and y co_ordinates of your neighbours
        fd dist
      ]

    ]
    [
      error 111
    ])

  ]
  [
    ; random walk...
    rt random-float 360  ;Turn to a random direction
    fd dist

  ]
end 

to get_infected ;S _> I
  let sick_neighbours infecteds in-radius max_inf_radius
  if (breed = susceptibles and count sick_neighbours  >= 1)
   [if ((random-float 1) <  (1 - ((1 - transmissibility_rate)^(count(sick_neighbours)))))
    [set breed infecteds
      set clock 0
      set tested? False]
      ]
end 

to recover ;I _> R
  if (clock >= infection_period and (random-float 1 < recovery_rate))
    [set breed recovereds]
end 

to lose_immunity ;R _> S
  if (random-float 1 < lose_immunity_rate)
  [set breed susceptibles]
end 

to advance_clock ;For recovery of infected people
  set clock (clock + 1)
end 

to runtest
  if (breed = infecteds)
  [die]
end 

to test;
  ifelse (num_to_test <= count(turtles))
    [let testsubjects n-of num_to_test turtles
    ask (testsubjects) [runtest]]
    [ask (turtles) [runtest]]
end 

;-------------------------------------------------------------------------------------------------------------------
; go to do various experiments and average
;-------------------------------------------------------------------------------------------------------------------

to go-experiments

  go

  if (count(infecteds) = 0 or count(infecteds) = count(turtles) or ticks > t_max) ;Stop if nobody is infected, or everybody is infected
  [
    ; if you enter here start a new experiment

    if (n_experiments = num_experiments)
    [
      file-close ;cierra el lammps_file
      if file-exists? output_file     
         [ file-delete output_file ]  
      file-open output_file
      file-print "# t infecteds error"
      set t 0
      while [t < t_max]
      [
        set s1 array:item acum t
        set s2 array:item acum2 t
        array:set sigma t sqrt ( s2 / num_experiments - ( (s1 / num_experiments) ^ 2 ) )
        file-write t
        file-write ( s1 / num_experiments )
        file-write ( ( array:item sigma t ) / sqrt(num_experiments ) )
        file-print ""
        set t (t + 1)
      ]
      file-close
      stop
    ]
    set n_experiments (n_experiments + 1)
    setup-experiment
  ]
end 

to lammps
  file-print "ITEM: TIMESTEP"
  file-print ticks
  file-print "ITEM: NUMBER OF ATOMS"
  file-print count turtles
  file-print "ITEM: BOX BOUNDS xx yy zz"
  file-write min-pxcor
  file-type " "
  file-print max-pxcor
  file-write min-pycor
  file-type " "
  file-print max-pycor
  file-print "0.0 0.0"
  file-print "ITEM: ATOMS id type x y z radius"
  ask (recovereds) [
    file-write who
    file-write 1
    file-type " "
    file-write xcor
    file-type " "
    file-write ycor
    file-type " "
    file-write 0.0
    file-type " "
    file-print SD_own
  ]
  ask (susceptibles) [
    file-write who
    file-write 2
    file-type " "
    file-write xcor
    file-type " "
    file-write ycor
    file-type " "
    file-write 0.0
    file-type " "
    file-print SD_own
  ]
  ask (infecteds) [
    file-write who
    file-write 3
    file-type " "
    file-write xcor
    file-type " "
    file-write ycor
    file-type " "
    file-write 0.0
    file-type " "
    file-print SD_own
  ]
end 

; go of an experiment

to go

  if (ticks mod freq_dump_file = 0)
     [lammps]

  ask (recovereds) [lose_immunity]
  ask (susceptibles) [get_infected]
  ask (infecteds) [advance_clock]
  ask (infecteds) [recover]
  ask (turtles) [move 1]
  ask (turtles) [recolor]
  if (ticks mod test_period = 0)
      [test]

  if ticks < t_max
  [
    array:set acum ticks  ( array:item acum ticks  + (100 * count(infecteds) / count(turtles) ) )
    array:set acum2 ticks ( array:item acum2 ticks + ((100 * count(infecteds) / count(turtles)) ^ 2) )

    set s1 array:item acum ticks
    set s2 array:item acum2 ticks
    array:set sigma ticks sqrt ( s2 / n_experiments - ( (s1 / n_experiments) ^ 2 ) )
  ]

  tick
end