...
 
Commits (3)
......@@ -20,29 +20,29 @@ end
-- Creates a new person, assigning processess
M.create_exposed = function ( zone, exposed_from, t )
M.create_exposed = function ( zone, exposed_from )
local t = lib.esim.now
local new_exposed = {
id=#people+1,
exposed_from = exposed_from,
exposed_to = {},
in_zone = zone,
quarantine = false,
disease = { infgen = 0 }
disease = lib.disease.new_disease()
}
people[new_exposed.id] = new_exposed
if type(exposed_from)=='table' then
exposed_from.exposed_to[#exposed_from.exposed_to+1] = new_exposed
end
-- SEIR process
lib.esim.launch(lib.disease.process, new_exposed, t)
-- disease process
lib.esim.launch(lib.disease.process, new_exposed)
--assert(new_exposed.disease and new_exposed.disease.infgen==0)
-- pick and launch a mobility task
local mobility_model = mobility.pick_mobility_model(t, zone, exposed_from)
new_exposed.mobility_model = mobility_model
lib.esim.launch(mobility.process[mobility_model](new_exposed, t))
lib.esim.launch(mobility.process[mobility_model](new_exposed))
w_log(t, ' new_exposed id: ', new_exposed.id,
' mobility: ', new_exposed.mobility_model,
......
......@@ -30,10 +30,9 @@ local function tracer_process ( t, tracer_id )
t = lib.esim.sleep_to(t+conf.contact_t_find)
if i<=#case.exposed_to then
local contact = case.exposed_to[i]
if not contact.quarantined then
if not contact.disease.quarantined then
w_log(t, ' traced tid: ', tracer_id, ' case: ', case.id, ' contact: ', contact.id )
case.exposed_to[i].quarantined = true
lib.disease.Qcount = lib.disease.Qcount+1
contact.disease.quarantined = true
end
i = i + 1
end
......@@ -52,7 +51,6 @@ M.initialize = function ( c )
conf = assert(c)
lib = conf.lib
if not conf.f_eir then w_slir = function() end end
if not conf.f_log then w_log = function() end end
M.cases = queue:new()
......
reset
clear
#set title ""
#set xrange [0:200]
set yrange [0:1]
set xlabel '#tracers'
set ylabel 'p'
#set key top right
set key autotitle columnhead
set term post size 5in, 2in eps enhanced color
#set term png size 700,400
set out "contacttracing.eps"
filename = 'out/contacttracing.log'
stats filename nooutput
#plot "out/states1.log" using 1:2 w steps, \
# "out/states1.log" using 1:3 w steps, \
# "out/states1.log" using 1:4 w steps, \
plot for [i=3:STATS_columns] filename using 1:i w l
unset output
......@@ -9,9 +9,11 @@ local conf
local lib
local model
local function w_eir (t)
local function w_states ()
local t = lib.esim.now
--w_eir('t E I Q R')
conf.f_eir:write(t, ' ', M.Ecount, ' ', M.Icount, ' ', M.Qcount, ' ', M.Rcount, '\n')
conf.f_eir:write(t, ' ', M.state_count.E, ' ', M.state_count.I, ' ',
M.state_count.R, ' ', M.state_count.quarantined, '\n')
end
local function w_log (t, ...)
conf.f_log:write(t, ...)
......@@ -19,7 +21,41 @@ local function w_log (t, ...)
--conf.f_log:flush()
end
M.infector = function (person, t, to_t)
local states_mt = {
__index = function(t, k )
if k == 'state' then
return t._state
elseif k == 'quarantined' then
return t._quarantined
end
return rawget(t, k)
end,
__newindex = function(t, k, v )
if k == 'state' then
if t._state then M.state_count[t._state] = M.state_count[t._state]-1 end
t._state = v
M.state_count[v] = M.state_count[v] + 1
if v=='R' and t._quarantined then
M.state_count.quarantined = M.state_count.quarantined - 1
end
w_states()
elseif k=='quarantined' then
if v then
assert(not t._quarantined)
M.state_count.quarantined = M.state_count.quarantined + 1
t._quarantined = v
w_states()
else
error('unquarantine not supported')
end
else
rawset(t, k, v)
end
end,
}
M.infector = function (person, to_t)
local t = lib.esim.now
local zone = person.in_zone
local d = person.disease
local infgen = d.infgen + 1
......@@ -27,71 +63,72 @@ M.infector = function (person, t, to_t)
local next_infection = t+mathutil.sample_exp_distribution( model.R0/model.T_I )
while next_infection<to_t and person.disease.state~='R'
and not person.quarantined do
and not d.quarantined do
t = lib.esim.sleep_to(next_infection)
if infgen ~= d.infgen then --there was another infector process launched
return lib.esim.sleep_to(to_t)
end
if person.disease.state == 'I' and not person.quarantined then
lib.agents.create_exposed(zone, person, t)
if person.disease.state == 'I' and not d.quarantined then
lib.agents.create_exposed(zone, person)
end
next_infection = t+mathutil.sample_exp_distribution( model.R0/model.T_I )
end
return lib.esim.sleep_to(to_t)
end
local selfQuarantineProcess = function ( person, t )
local selfQuarantineProcess = function ( person )
local t = lib.esim.now
if not model.T_q then return end
local d = person.disease
local tq = t+model.T_q
lib.esim.sleep_to(tq)
if not person.quarantined then
person.quarantined = true
M.Qcount = M.Qcount+1
if not d.quarantined then
d.quarantined = true
w_log(t, ' selfquarantined id:', person.id)
lib.contacttracer.cases:pushright(person)
end
end
M.process = function ( person, t )
if M.Ecount + M.Icount > conf.outbreak_threshold then
M.process = function ( person )
local t = lib.esim.now
if M.state_count.E + M.state_count.I > conf.outbreak_threshold then
lib.esim.finish('outbreak')
return
end
local d = person.disease
d.state = 'E'
M.Ecount = M.Ecount + 1
w_eir(t)
t = lib.esim.sleep_to( t+model.T_E )
d.state = 'I'
M.Ecount, M.Icount = M.Ecount-1, M.Icount+1
w_eir(t)
lib.esim.launch(selfQuarantineProcess, person, t)
lib.esim.launch(selfQuarantineProcess, person)
local time_start_R = t+model.T_I
t = M.infector(person, t, time_start_R)
t = M.infector(person, time_start_R)
d.state = 'R'
M.Icount, M.Rcount = M.Icount-1, M.Rcount+1
if person.quarantined then
M.Qcount = M.Qcount-1
end
w_eir(t)
if M.Ecount + M.Icount == 0 then
if M.state_count.E + M.state_count.I == 0 then
lib.esim.finish('extinction')
end
end
M.new_disease = function ()
local d = {
infgen = 0,
}
d = setmetatable(d, states_mt)
return d
end
M.initialize = function ( c )
conf = assert(c)
model = c.disease
lib = conf.lib
M.Ecount, M.Icount, M.Rcount, M.Qcount = 0, 0, 0, 0
M.state_count = {E=0, I=0, R=0, quarantined=0}
if not conf.f_eir then w_eir = function() end end
if not conf.f_eir then w_states = function() end end
if not conf.f_log then w_log = function() end end
w_eir('t E I Q R')
if conf.f_eir then conf.f_eir:write('t L P I R quarantined\n') end
end
return M
\ No newline at end of file
......@@ -9,10 +9,11 @@ local conf
local lib
local model
local function w_slir (t)
local function w_states ()
local t = lib.esim.now
--w_slir('t L P I Q R')
conf.f_eir:write(t, ' ', M.Lcount, ' ', M.Pcount, ' ', M.Icount
, ' ', M.Qcount, ' ', M.Rcount, '\n')
conf.f_eir:write(t, ' ', M.state_count.L, ' ', M.state_count.P, ' ',
M.state_count.I, ' ', M.state_count.R, ' ', M.state_count.quarantined, '\n')
conf.f_eir:flush()
end
local function w_log (t, ...)
......@@ -21,6 +22,40 @@ local function w_log (t, ...)
--conf.f_log:flush()
end
local states_mt = {
__index = function(t, k )
if k == 'state' then
return t._state
elseif k == 'quarantined' then
return t._quarantined
end
return rawget(t, k)
end,
__newindex = function(t, k, v )
if k == 'state' then
if t._state then M.state_count[t._state] = M.state_count[t._state]-1 end
t._state = v
M.state_count[v] = M.state_count[v] + 1
if v=='R' and t._quarantined then
M.state_count.quarantined = M.state_count.quarantined - 1
end
w_states()
elseif k=='quarantined' then
if v then
assert(not t._quarantined)
M.state_count.quarantined = M.state_count.quarantined + 1
t._quarantined = v
w_states()
else
error('unquarantine not supported')
end
else
rawset(t, k, v)
end
end,
}
local function compute_beta_SLIR (d)
local beta
if d.state == 'I' then
......@@ -35,7 +70,8 @@ local function compute_beta_SLIR (d)
return beta
end
M.infector = function (person, t, to_t)
M.infector = function (person, to_t)
local t = lib.esim.now
local zone = person.in_zone
local d = person.disease
local infgen = d.infgen + 1
......@@ -47,14 +83,14 @@ M.infector = function (person, t, to_t)
end
local next_infection = t + mathutil.sample_exp_distribution( beta )
while next_infection<to_t and person.disease.state~='R'
and not person.quarantined do
and not d.quarantined do
t = lib.esim.sleep_to(next_infection)
if infgen ~= d.infgen then --there was another infector process launched
return lib.esim.sleep_to(to_t)
end
if person.disease.state == 'P' or person.disease.state == 'I'
and not person.quarantined then
lib.agents.create_exposed(zone, person, t)
and not d.quarantined then
lib.agents.create_exposed(zone, person)
end
local beta = compute_beta_SLIR(d)
if not beta then
......@@ -65,81 +101,74 @@ M.infector = function (person, t, to_t)
return lib.esim.sleep_to(to_t)
end
local selfQuarantineProcess = function ( person, t )
local selfQuarantineProcess = function ( person )
local t = lib.esim.now
if not model.T_q then return end
local d = person.disease
local tq = t+model.T_q
lib.esim.sleep_to(tq)
if not person.quarantined then
person.quarantined = true
M.Qcount = M.Qcount+1
if not d.quarantined then
d.quarantined = true
w_log(t, ' selfquarantined id:', person.id)
lib.contacttracer.cases:pushright(person)
end
end
M.process = function ( person, t )
if M.Lcount + M.Pcount + M.Icount >= conf.outbreak_threshold then
M.process = function ( person )
local t = lib.esim.now
if M.state_count.L + M.state_count.P + M.state_count.I
> conf.outbreak_threshold then
lib.esim.finish('outbreak')
return
end
local d = person.disease
d.state = 'L'
M.Lcount = M.Lcount + 1
d.symptomatic = math.random()>model.p
w_slir(t)
if d.symptomatic then
t = lib.esim.sleep_to( t+model._eps )
d.state = 'P'
M.Lcount, M.Pcount = M.Lcount-1, M.Pcount+1
w_slir(t)
local time_start_I = t+model._gamma
t = M.infector(person, t, time_start_I)
t = M.infector(person, time_start_I)
d.state = 'I'
M.Pcount, M.Icount = M.Pcount-1, M.Icount+1
w_slir(t)
lib.esim.launch(selfQuarantineProcess, person, t)
lib.esim.launch(selfQuarantineProcess, person)
local time_start_R = t+model._mu
t = M.infector(person, t, time_start_R)
t = M.infector(person, time_start_R)
d.state = 'R'
M.Icount, M.Rcount = M.Icount-1, M.Rcount+1
if person.quarantined then
M.Qcount = M.Qcount-1
end
w_slir(t)
else --asymptomatic
t = lib.esim.sleep_to( t+model._epsA )
d.state = 'I'
M.Lcount, M.Icount = M.Lcount-1, M.Icount+1
w_slir(t)
local time_start_R = t+model._mu
t = M.infector(person, t, time_start_R)
t = M.infector(person, time_start_R)
d.state = 'R'
M.Icount, M.Rcount = M.Icount-1, M.Rcount+1
if person.quarantined then
M.Qcount = M.Qcount-1
end
w_slir(t)
end
if M.Lcount + M.Pcount + M.Icount == 0 then
if M.state_count.L + M.state_count.P + M.state_count.I == 0 then
lib.esim.finish('extinction')
end
end
M.new_disease = function ()
local d = {
infgen = 0,
}
d = setmetatable(d, states_mt)
return d
end
M.initialize = function ( c )
conf = assert(c)
model = c.disease
lib = conf.lib
M.Lcount, M.Pcount, M.Icount, M.Rcount, M.Qcount = 0, 0, 0, 0, 0
if not conf.f_eir then w_slir = function() end end
M.state_count = {L=0, P=0, I=0, R=0, quarantined=0}
if not conf.f_eir then w_states = function() end end
if not conf.f_log then w_log = function() end end
w_slir('t L P I Q R')
if conf.f_eir then conf.f_eir:write('t L P I R quarantined\n') end
end
return M
\ No newline at end of file
......@@ -12,7 +12,7 @@ local conf = {
stop_localidad = nil, -- any
outbreak_threshold = 100, -- Max number of actives
---[[
--[[
disease = {
model = 'SLIR',
R0 = 2.5,
......@@ -29,7 +29,7 @@ local conf = {
T_q = 3.0, -- time to self quarantine
},
--]]
--[[
---[[
disease = {
model = 'SEIR',
R0 = 2,
......@@ -63,9 +63,12 @@ conf.lib = lib
local scenario_process = function (conf)
local G = require('lib.graph')
for zoneId, Exposition in pairs(conf.startExposed) do
for i = 1, Exposition.count do
lib.agents.create_exposed (G.N[zoneId], nil, Exposition.t)
end
lib.esim.launch(function()
for i = 1, Exposition.count do
lib.esim.sleep_to(Exposition.t)
lib.agents.create_exposed (G.N[zoneId], nil)
end
end)
end
end
......@@ -104,7 +107,7 @@ for i = 1, NRUNS do
end
until sim_status
conf.f_log:write(t_now, ' finished ',sim_status,' R: '..lib.disease.Rcount..'\n')
conf.f_log:write(t_now, ' finished ',sim_status,' R: '..lib.disease.state_count.R..'\n')
conf.f_eir:close()
conf.f_log:close()
......
--------------------------------------------------------------------------
-- parameters
local conf = {
startExposed = {
[13943] = {t=0, count=20}, -- exposed at t=0 , RIVERA-MASOLLER
},
--stop conditions
stop_dpto = 'MONTEVIDEO',
stop_localidad = nil, -- any
outbreak_threshold = 100, -- Max number of actives
---[[
disease = {
model = 'SLIR',
R0 = 2.5,
-- https://www.medrxiv.org/content/10.1101/2020.05.06.20092841v1
p = 0.25, --proportion of asymptomatic
_epsA = 5, --latent period
_eps = 3, --latent period
_gamma = 2, --pre-symptomatic period
_mu = 6, --2.5, --time to removed/home stay
k = 0.15, -- proportion of presymptomatic transmission
r = 0.5, -- relative infectiousness of asymptomatic individuals
T_q = 3.0, -- time to self quarantine
},
--]]
--[[
disease = {
model = 'SEIR',
R0 = 2,
T_E = 6, -- Duration of E state
T_I = 6, -- Duration of I state
T_q = 2.5, -- time to self quarantine
},
--]]
contacttracers = 0,
contact_t_find = 1,
contacttracing_time = 5,
}
local TMAX = math.huge --12*30 -- Simulation time limit
local NRUNS = 1000 -- Number of runs
--------------------------------------------------------------------------
--require('mobdebug').coro()
local lib = {}
lib.esim = require('lib.esim')
lib.agents = require('agents')
lib.disease = require('disease_'..conf.disease.model)
lib.contacttracer = require('contacttracer')
conf.lib = lib
--------------------------------------------------------------------------
-- process that sets the scenario up.
local scenario_process = function (conf)
local G = require('lib.graph')
for zoneId, Exposition in pairs(conf.startExposed) do
lib.esim.launch(function()
for i = 1, Exposition.count do
lib.esim.sleep_to(Exposition.t)
lib.agents.create_exposed (G.N[zoneId], nil)
end
end)
end
end
--os.execute('rm out/*.log out/*.txt')
local function runsim()
local sim_status_count = {outbreak=0, extinction=0, runout=0}
for i = 1, NRUNS do
-- Log files
--conf.f_eir = io.open('out/states'..i..'.log', 'w')
--conf.f_log = io.open('out/epitrace'..i..'.log', 'w')
-- Seed the random generator independetly per run
math.randomseed(i) --os.time())
--conf.f_log:write(0, ' randomseed ', i, '\n')
lib.esim.initialize( TMAX )
lib.agents.initialize( conf )
lib.disease.initialize( conf )
lib.contacttracer.initialize( conf )
lib.esim.launch( scenario_process or function() end, conf )
local event_count = 0
local t_now
local sim_status
repeat
local t, status = lib.esim.step()
if t then
event_count = event_count + 1
t_now = t
end
-- stop conditions
if not t then
sim_status = tostring(status)
end
until sim_status
--conf.f_log:write(t_now, ' finished ',sim_status,' R: '..lib.disease.Rcount..'\n')
--conf.f_eir:close()
--conf.f_log:close()
sim_status_count[sim_status] = (sim_status_count[sim_status] or 0) + 1
print (i, 'finished at t: '..t_now..' events: '
..event_count..' ('..sim_status ..')')
end
-- Print run totals
--for status, count in pairs(sim_status_count) do
-- print('runs', status, count, count/NRUNS)
--end
-- Launch plotting (on background)
--os.execute('gnuplot states.plot &')
--os.execute('luajit parse_logs.lua '..NRUNS..' &')
return sim_status_count
end
local frunlog = io.open('out/contacttracing.log', 'w')
frunlog:write('tracers runout outbreak extinction hitstop\n')
for ntracers = 0, 20, 2 do
conf.contacttracers = ntracers
local s = runsim()
frunlog:write(
ntracers, ' ',
s.runout/NRUNS, ' ',
s.outbreak/NRUNS, ' ',
s.extinction/NRUNS, ' ',
s.hit_stop/NRUNS, '\n'
)
end
frunlog:write()
os.execute('gnuplot contacttracing.plot &')
......@@ -5,8 +5,11 @@ local TMAX
local finished_condition =nil
M.now = 0
M.initialize = function ( tmax )
TMAX = assert(tmax)
M.now = 0
finished_condition = nil
events:initialize()
end
......@@ -48,6 +51,7 @@ M.step = function ()
end
local co, t_now = events:pop()
M.now = t_now
if coroutine.status(co) ~= 'dead' then
local ok, err = coroutine.resume(co, t_now)
if not ok then
......
......@@ -75,8 +75,9 @@ local function next_hour_of_day (now, hour)
end
M.process.random_bus_day = function ( person, t )
M.process.random_bus_day = function ( person )
local f = function ()
local t = lib.esim.now
local home, work
local next_go = next_hour_of_day(t, 7)
local next_ret = next_hour_of_day(t, 17)
......@@ -84,7 +85,7 @@ M.process.random_bus_day = function ( person, t )
if next_ret<next_go then -- started at work
work = person.in_zone
lib.disease.infector(person, t, next_ret)
lib.disease.infector(person, next_ret)
if person.disease.state=='R' then return end
t = esim.sleep_to(next_ret)
......@@ -99,7 +100,7 @@ M.process.random_bus_day = function ( person, t )
while true do
local next_go = next_hour_of_day(t, 7)
t = lib.disease.infector(person, t, next_go)
t = lib.disease.infector(person, next_go)
if person.disease.state=='R' then return end
work = work or pick_random_destination( person, t )
......@@ -108,7 +109,7 @@ M.process.random_bus_day = function ( person, t )
person.in_zone = work
local next_ret = next_hour_of_day(t, 17)
t = lib.disease.infector(person, t, next_ret)
t = lib.disease.infector(person, next_ret)
if person.state=='R' then return end
home = home or pick_random_source( person, t ) -- ponderated source of workers for work
......@@ -120,9 +121,9 @@ M.process.random_bus_day = function ( person, t )
return f
end
M.process.static = function ( person, t )
M.process.static = function ( person )
local f = function ()
lib.disease.infector(person, t, math.huge)
lib.disease.infector(person, math.huge)
end
return f
end
......