-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------


                       Example: The spread of a disease


    In this model each agent can be in one of three states:


        Susceptible: the agent is healthy, but can get the disease.
        Infectious: the agent has the disease and can spread it to others.
        Recovered: the agent has had the disease and is now immune.


    In the first round one or more agents get infected.


    In all subsequent rounds:


        - All infectious agents have a chance to become recovered, and
        - all recovered agents get a chance to become susceptible.
        - All agents that are still infectious contact a random number of
        other agents
        - Each of these contact has a chance of resulting in an infection if
        someone is susceptible.

. set seed 1234567


. set scheme s1color


. clear mata


. drop _all


. run sir_main.do


. 
. mata:
------------------------------------------------- mata (type end to exit) -----
: 
: model = sir()


: model.N(2000)


: model.tdim(150)


: model.outbreak(5)


: model.mcontacts(4)


: model.transmissibility(0.045)


: model.mindur(10)


: model.meandur(14)


: model.pr_loss(0)


: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................        50
..................................................       100
..................................................       150


: model.export_sir()


: end
-------------------------------------------------------------------------------


. 
. twoway line S I R t,  name(sir1, replace)          ///
>    ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..) 


.    
. drop _all

    The number 0.045 for transmissibility was chosen to get a >> R0of
    2.5, which in line with COVID-19.


 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------


                          Social distancing and masks


    The purpose of social distancing and encouraging masks is to reduce the
    chance of a contact leading to an infection.


    We can change the model such that these policies reduce the
    transmissibility by a quarter

. mata:
------------------------------------------------- mata (type end to exit) -----
: model.transmissibility(0.03)


: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................        50
..................................................       100
..................................................       150


: model.export_sir()


: end
-------------------------------------------------------------------------------


. 
. twoway line S I R t,  name(sir2, replace)          ///
>    ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..) 


. drop _all

    We can similarly study the effect of a lock-down, which reduces the
    number of people we contact


    or we can adept the model such that we look at how contact tracing might
    influence the spread of the disease


 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------


                                Losing immunity


    Recently, there have an increasing number of news stories about that
    immunity against COVID-19 does not persist


    We can change to model to see how that impacts the spread of the disease

. mata:
------------------------------------------------- mata (type end to exit) -----
: model.transmissibility(0.045)


: model.pr_loss(0.02)


: model.tdim(250)


: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................        50
..................................................       100
..................................................       150
..................................................       200
..................................................       250


: model.export_sir()


: end
-------------------------------------------------------------------------------


. 
. twoway line S I R t,  name(sir3, replace)          ///
>    ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..) 


. drop _all

 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------


                               Adding a network


    The model thus far does not need a simulation; one can compute them with
    a set of differential equations.


    In Stata you can do that with epimodels by Sergiy Radyakin and Paolo
    Verme


    However, the advantage of Agent Based Models is that it is much easier to
    expand the model.


    For example, what if we don't believe that people contact each other
    randomly but instead it is a network that determines who contacts who.


    A common model that works fairly well for a lot of social networks is a
    small world network:


        We assume all agents are on a circle, and that each agent is
        connected to the 4 closes agents.


        However, a small number of these connections are replaced by a
        connection to a random agent

. clear mata


. run nw_main.do


. mata:
------------------------------------------------- mata (type end to exit) -----
: model = sir_nw()


: model.N(100)


: model.tdim(10)


: model.outbreak(2)


: model.degree(4)


: model.pr_long(.05)


: model.transmissibility(0.1)


: model.mindur(5)


: model.meandur(6)


: model.pr_loss(0.00)


: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..........
: model.export_nw()


: end
-------------------------------------------------------------------------------


. 
. gen x_ego = cos((ego-1)/`n'*2*_pi)*(1+mod(ego,2)*.2)


. gen y_ego = sin((ego-1)/`n'*2*_pi)*(1+mod(ego,2)*.2)


. gen x_alter = cos((alter-1)/`n'*2*_pi)*(1+mod(alter,2)*.2)
(1 missing value generated)


. gen y_alter = sin((alter-1)/`n'*2*_pi)*(1+mod(alter,2)*.2)
(1 missing value generated)


.     
. 
. bys ego (alter) : gen first = _n == 1


. 
. forvalues t = 1/10{
  2.     twoway pcspike y_ego x_ego y_alter x_alter,                ///
>                   lcolor(gs8)                               || ///
>         scatter y_ego x_ego if first & state`t'==1 ,           ///
>                    mcolor(yellow) msymbol(O)                || ///
>         scatter y_ego x_ego if first & state`t'==2 ,           ///
>                    mcolor(red) msymbol(O)                   || ///
>         scatter y_ego x_ego if first & state`t'==3 ,           ///
>                    mcolor(black) msymbol(Oh)                   ///         
>                    aspect(1) xscale(off) yscale(off)           ///
>                    legend(order(2 "infected"                   ///
>                                 3 "susceptible"                ///
>                                 4 "recovered") rows(1))        ///
>                    name(t`t', replace) title(Time `t')     
  3. }


. drop _all

 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
What is an Agent Based Model?
-------------------------------------------------------------------------------


                          Avoiding unnecesary travel


    Close knit networks tend to slow down the spread of a disease, while the
    long distance ties tend


    So that is another way in which advise against unnecesary travel or
    "social buble" rules in countries like Belgium works


 

. mata:
------------------------------------------------- mata (type end to exit) -----
: model.N(2000)


: model.tdim(150)


: model.outbreak(10)


: model.degree(10)


: model.pr_long(.2)


: model.transmissibility(.018)


: model.mindur(10)


: model.meandur(14)


: model.pr_loss(0.00)


: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................        50
..................................................       100
..................................................       150


: model.export_sir()


: end
-------------------------------------------------------------------------------


. 
. twoway line S I R t,  name(sir4, replace)          ///
>    ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..) 


. drop _all


. 
. mata:
------------------------------------------------- mata (type end to exit) -----
: model.pr_long(.05)


: model.run()
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5
..................................................        50
..................................................       100
..................................................       150


: model.export_sir()


: end
-------------------------------------------------------------------------------


. 
. twoway line S I R t,  name(sir5, replace)          ///
>    ylabel(,angle(0)) legend(rows(1)) lwidth(*1.5 ..) 


. 
. drop _all

 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------


                                Where to start


    We have a population of agents to keep track of


    We need to repeatedly apply rules to them


    Applying rules can be done by a function, but that function needs access
    to settings and the population.


    This is a good case for a class.


    A class can be thought of as a tray with named slots for functions and
    data


    The functions are local to the class


    Each function has access to all material stored in the class


 

. clear mata


. 
. mata:
------------------------------------------------- mata (type end to exit) -----
: 
: class arithmatic {
> 
>     real scalar a
>     real scalar b
>     
>     real scalar sum()
>     real scalar prod()
>     void input()
>     void run()
> }


: 
: void arithmatic::input(real scalar aval, real scalar bval)
> {
>     a = aval
>     b = bval
> }


: 
: real scalar arithmatic::sum()
> {
>     return(a + b)
> }


: 
: real scalar arithmatic::prod()
> {
>     return(a * b)
> }


: 
: void arithmatic::run()
> {
>     sum()
>     prod()
> }


: 
: // how the user interacts with this program
: math = arithmatic()


: math.input(2,4)


: math.run()
  6
  8


: end
-------------------------------------------------------------------------------

 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------


                                Basic SIR model


    Lets take a look at the code for the basic SIR model:

version 16.1

run abm_pop.mata

mata:

mata set matastrict on

class sir {
	class abm_pop          scalar agents
	
	real                   scalar tdim
	real                   scalar outbreak
	real                   scalar removed
	real                   scalar mcontacts
	real                   scalar N
	real                   scalar transmissibility
	real                   scalar mindur
	real                   scalar meandur
	real                   scalar pr_heal
	real                   scalar pr_loss
	
	// sir_chks.do
	void                          posint()
	void                          pr()
	void                          isid()
	void                          istime()
	
	// sir_set_pars.do
	transmorphic                  N()
	transmorphic                  tdim()
	transmorphic                  outbreak()
	transmorphic                  removed()
	transmorphic                  mcontacts()
	transmorphic                  transmissibility()
	transmorphic                  mindur()
	transmorphic                  meandur()
	transmorphic                  pr_loss()

	// sir_sim.do
	void                          run()
	void                          setup()
	void                          progress()
	void                          step()
	real                   matrix meet()
	real                   scalar infect()
	void                          dots()
	
	// sir_export.do
	void                          export_sir()
	void                          export_r()
}
end 

do sir_chks.do
do sir_set_pars.do
do sir_sim.do
do sir_export.do

exit

    The population is stored is stored in another class abm_pop, which is
    available from https://github.com/maartenteaches/abm_pop.


    Whatever is stored in abm_pop is assumed to persist until you store
    something else. So you only have to store the disease state of an agent
    when it changes. This can save a lot of reading and writing and make the
    simulation run quicker.


    There are various checks perfomed and we can see those

include sir_locals.do

mata:
void sir::posint(real scalar val)
{
	if (val <= 0 | val != ceil(val)) {
		_error("argument must be a positive integer")
	} 
}

void sir::pr(real scalar val)
{
	if (val <0 | val > 1) {
		_error("argument must be between 0 and 1")
	}
}

void sir::isid(real scalar id)
{
	if (id < 0 | id > N | id!=ceil(id)) {
		_error("argument is not a valid id")
	}
	
}

void sir::istime(real scalar t)
{
	if (t<0 | t > tdim | t!=ceil(t)) {
		_error("argument is not a valid time")
	}
}
end

    The functions that set the parameters are shown

include sir_locals.do

mata:

transmorphic sir::N(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		N = val
	}
	else {
		return(N)
	}
}

transmorphic sir::tdim(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		tdim= val
	}
	else {
		return(tdim)
	}
}


transmorphic sir::outbreak(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		outbreak= val
	}
	else {
		return(outbreak)
	}
}

transmorphic sir::removed(| real scalar val)
{
	if (args() == 1) {
		posint(val)
		removed= val
	}
	else {
		return(removed)
	}
}

transmorphic sir::mcontacts(| real scalar val)
{
	if ( args()==1 ) {
		if ( val < 0) _error("arguments must be positive")
		mcontacts = val
	}
	else {
		return(mcontacts)
	}
}

transmorphic sir::transmissibility(| real scalar val)
{
	if ( args()==1 ) {
		pr(val)
		transmissibility = val
	}
	else {
		return(transmissibility)
	}
}

transmorphic sir::mindur(| real scalar val)
{
	if ( args()==1 ) {
		posint(val)
		mindur = val
	}
	else {
		return(mindur)
	}
}

transmorphic sir::meandur(| real scalar val)
{
	if ( args()==1 ) {
		if (val < 0) _error("argument must be positive")
		meandur = val
	}
	else {
		return(meandur)
	}
}

transmorphic sir::pr_loss(| real scalar val)
{
	if ( args()==1 ) {
		pr(val)
		pr_loss = val
	}
	else {
		return(pr_loss)
	}
}


end

    The main simulation functions are

include sir_locals.do

mata:

void sir::run()
{
	real scalar i, t
	setup()

	dots(0)
	dots(1)

	for (t=2; t<=tdim; t++) {
		dots(t)
		for(i=1; i<=N ; i++) {
			step(i,t)
		}
	}
}

void sir::setup()
{
	real vector id, key
	real scalar i
	
	if (N==.)                   _error("N has not been set")
	if (tdim==.)                _error("tdim has not been set")
	if (outbreak==.)            _error("outbreak has not been set")
	if (outbreak > N)           _error("the initial outbreak is larger than population")
	if (removed==.)             removed = 0
	if (removed + outbreak > N) _error("intial outbreak + removed agents is larger than the population")
	if (mcontacts==.)           _error("mcontacts has not been set")
	if (mcontacts > N)          _error("average number of mcontacts is larger than population")
	if (transmissibility==.)    _error("transmissibility has not been set")
	if (mindur==.)              _error("mindur has not been specified")
	if (meandur==.)             _error("meandur has not been specified")
	if (meandur <= mindur)      _error("meandur must be larger than mindur")
	if (pr_loss==.)             pr_loss = 0
	
	agents.N(N)
	agents.k(3)
	agents.setup()
	
	id = jumble((1..N)')
	
	for (i = 1; i<= N; i++) {
		key = id[i],`state',1
		agents.put(key,(i<=outbreak ? `infectious' : (i <= outbreak+removed ? `removed' : `susceptible')))
		key[2] = `exposes'
		if (i <= outbreak) agents.put(key, J(2,0,.) )
		
		key[2] = `dur'
		agents.put(key,1)
	}

	
// if meandur - mindur < 1 then pr_heal > 1, so there is no randomness to the duration
// and the fixed duration of the disease is mindur	
	pr_heal = 1/(meandur-mindur)
}

void sir::step(real scalar id, real scalar t)
{
	real matrix exposed
	real scalar i, k
	real rowvector key
	
	isid(id)
	istime(t)

	key = id, `dur',t

	agents.put(key, agents.get(key, "last")+1)
	progress(id,t)
	key = id, `state',t
	if (agents.get(key, "last") == `infectious' ) {
		exposed = meet(id, t)
		for(i=1; i<= cols(exposed); i++) {
			exposed[2,i] = infect(exposed[1,i],t) 
		}
		key[2] = `exposes'
		agents.put(key,exposed)
	}
}

void sir::progress(real scalar id, real scalar t)
{
    real rowvector key1, key2, key3
	
	isid(id)
	istime(t)
	
	key1 = id, `state',t
	key2 = id, `dur',t
	if (agents.get(key1, "last") == `infectious'  & agents.get(key2, "last") > mindur) {
		if (runiform(1,1) < pr_heal) {
			agents.put(key1, `removed')
			agents.put(key2,1)
			key3 = id, `exposes',t
			agents.put(key3,NULL)
		}
	}
    if (agents.get(key1, "last") == `removed' & runiform(1,1) < pr_loss) {
		agents.put(key1, `susceptible')
		agents.put(key2,1)
	}
}

real matrix sir::meet(real scalar id, real scalar t)
{
	real scalar k
	real colvector res
	real rowvector key

	isid(id)
	
	k = rpoisson(1,1,mcontacts)
	if (k == 0) {
		return(J(2,0,.))
	}
	else{
		res =  ceil(runiform(1,k):*(N-1))
		res = res + (res:>=id)
		res = res \ J(1,k,.)
		return (res)
	}
}

real scalar sir::infect(real scalar id, real scalar t)
{
	real vector key
	real scalar infected

	isid(id)
	istime(t)

	infected = 0
	key = id,`state', t
	if (agents.get(key, "last")==`susceptible') {
		if (runiform(1,1) < transmissibility) {
			agents.put(key, `infectious')
			key[2] = `dur'
			agents.put(key, 1)
			key[2] = `exposes'
			agents.put(key,J(2,0,.))
			infected = 1
		}
	}
	return(infected)
}



void sir::dots(real scalar t) 
{
	if (t==0) {
		printf("----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5\n")
	}
	else if (mod(t, 50) == 0) {
		printf(". %9.0f\n", t)
	}
	else {
		printf(".")
	}
	displayflush()
}
end

    The "variables" in agents are actually numbered and not named. To make it
    easier to read the code I use local macros defined in

// address in abm_pop class agents
local state       = 1
local dur         = 2
local exposes     = 3

// states
local susceptible = 1
local infectious  = 2 
local removed     = 3

    The functions that allow you to export the results to Stata are

include sir_locals.do

mata:

void sir::export_sir()
{
    real matrix res
	pointer matrix p
    string rowvector varnames
	real scalar xid, n, i, j

	res = J(N,tdim,.)
	p = agents.extract(`state', tdim)
	for (i=1; i<=N; i++) {
		for(j=1; j<=tdim; j++) {
			res[i,j] = *p[i,j]
		}
	}
	res =  (1..tdim)', mean(res:==`susceptible')', 
	                   mean(res:==`infectious')', 
					   mean(res:==`removed')'
    varnames = "t", "S", "I", "R"
    xid = st_addvar("float", varnames)
    n = st_nobs()
    n = rows(res) - n
    if (n > 0) {
        st_addobs(n)
    }
    st_store(.,xid, res)
}

void sir::export_r()
{
	real matrix res
	pointer matrix p_exposes, p_state
	real scalar xid, n, i, j, t
	string rowvector varnames
	real colvector n_ill

	t = .
	
	res = (1..tdim)',J(tdim,1,0)
	n_ill = J(tdim,1,0)
	p_exposes = agents.extract(`exposes',tdim)
	p_state = agents.extract(`state',tdim)
	for (i=1; i<=N; i++) {
		if (*p_state[i,1]==`infectious') {
			t=1  
		} 
		else {
		    t=.
		}
		for(j=2;j<=tdim;j++) {
			if (*p_state[i,j-1]==`susceptible' & 
			    *p_state[i,j]  ==`infectious') {
					t=j
					n_ill[t] = n_ill[t] + 1
				}
			if (*p_state[i,j-1]==`infectious' & 
			    (*p_state[i,j]  ==`removed' | *p_state[i,j]  ==`susceptible')) t = .
			if (t!=.) {
			    res[t,2] = res[t,2] + rowsum((*p_exposes[i,j])[2,.])
			}
		}
	}
	res[.,2] = res[.,2]:/n_ill
	varnames = "t", "reprod"
    xid = st_addvar("float", varnames)
    n = st_nobs()
    n = rows(res) - n
    if (n > 0) {
        st_addobs(n)
    }
    st_store(.,xid, res)
	
}


end

 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------

-------------------------------------------------------------------------------
How to implement an ABM in Mata
-------------------------------------------------------------------------------


                            SIR model in a network


    To add the network we make use the abm_nw class, which is available from 
    https://github.com/maartenteaches/abm_nw.


    This takes care of creating the network and keeps track of who can
    contact who.


    With that the changes required are fairly small


    The class definition is

version 16.1

run abm_pop.mata
run abm_nw.mata                                          // <-- new

mata:

mata set matastrict on

class sir_nw {
	class abm_pop          scalar agents
	class abm_nw           scalar nw                     // <-- new
	
	real                   scalar tdim
	real                   scalar outbreak
	real                   scalar removed
//	real                   scalar mcontacts              // <-- no longer necessary
	real                   scalar N
	real                   scalar transmissibility
	real                   scalar mindur
	real                   scalar meandur
	real                   scalar pr_heal
	real                   scalar pr_loss
	real                   scalar degree                 // <-- new
	real                   scalar pr_long                // <-- new
	
	// sir_chks.do
	void                          posint()
	void                          pr()
	void                          isid()
	void                          istime()
	
	// sir_set_pars.do
	transmorphic                  N()
	transmorphic                  tdim()
	transmorphic                  outbreak()
	transmorphic                  removed()
//	transmorphic                  mcontacts()           // <-- no longer necessary
	transmorphic                  transmissibility()
	transmorphic                  mindur()
	transmorphic                  meandur()
	transmorphic                  pr_loss()
	transmorphic                  degree()              // <-- new
	transmorphic                  pr_long()             // <-- new

	// sir_sim.do
	void                          setup()
//	real                   matrix meet()                // <-- no longer necessary
	real                   scalar infect()
	void                          progress()
	void                          step()
	void                          run()
	void                          dots()
	
	// sir_export.do
	void                          export_sir()
	void                          export_r()
	void                          export_nw()          // <-- new
}
end 

do nw_chks.do
do nw_set_pars.do
do nw_sim.do
do nw_export.do

exit

    Checks and setting parameters are pretty much the same


    The main simulation functions are

include nw_locals.do

mata:

void sir_nw::setup()
{
	real vector id, key
	real scalar i
	
	if (N==.)                   _error("N has not been set")
	if (tdim==.)                _error("tdim has not been set")
	if (outbreak==.)            _error("outbreak has not been set")
	if (outbreak > N)           _error("the initial outbreak is larger than population")
	if (removed==.)             removed = 0
	if (removed + outbreak > N) _error("intial outbreak + removed agents is larger than the population")
	if (transmissibility==.)    _error("transmissibility has not been set")
	if (mindur==.)              _error("mindur has not been specified")
	if (meandur==.)             _error("meandur has not been specified")
	if (meandur <= mindur)      _error("meandur must be larger than mindur")
	if (pr_loss==.)             pr_loss = 0
	if (degree==.)              _error("degree has not been specified")  // <-- new
	if (pr_long==.)             _error("pr_long has not been specified") // <-- new
	
	nw.clear()                                                           // <-- new   
	nw.N_nodes(0,N)                                                      // <-- new
	nw.directed(0)                                                       // <-- new
	nw.tdim(0)                                                           // <-- new
	nw.sw(degree,pr_long)                                                // <-- new    
	nw.setup()                                                           // <-- new
	
	agents.N(N)
	agents.k(3)
	agents.setup()
	
	id = jumble((1..N)')
	
	for (i = 1; i<= N; i++) {
		key = id[i],`state',1
		agents.put(key,(i<=outbreak ? `infectious' : (i <= outbreak+removed ? `removed' : `susceptible')))
		key[2] = `exposes'
		if (i <= outbreak) agents.put(key, J(2,0,.) )
		
		key[2] = `dur'
		agents.put(key,1)
	}

	
// if meandur - mindur < 1 then pr_heal > 1, so there is no randomness to the duration
// and the fixed duration of the disease is mindur	
	pr_heal = 1/(meandur-mindur)
}


real scalar sir_nw::infect(real scalar id, real scalar t)
{
	real vector key
	real scalar infected

	isid(id)
	istime(t)

	infected = 0
	key = id,`state', t
	if (agents.get(key, "last")==`susceptible') {
		if (runiform(1,1) < transmissibility) {
			agents.put(key, `infectious')
			key[2] = `dur'
			agents.put(key, 1)
			key[2] = `exposes'
			agents.put(key,J(2,0,.))
			infected = 1
		}
	}
	return(infected)
}

void sir_nw::step(real scalar id, real scalar t)
{
	real matrix exposed
	real scalar i, k
	real rowvector key
	
	isid(id)
	istime(t)

	key = id, `dur',t

	agents.put(key, agents.get(key, "last")+1)
	progress(id,t)
	key = id, `state',t
	if (agents.get(key, "last") == `infectious' ) {
		exposed = nw.neighbours(id)                        // <-- new
		exposed = exposed \ J(1,cols(exposed),.)           // <-- new
		for(i=1; i<= cols(exposed); i++) {
			exposed[2,i] = infect(exposed[1,i],t) 
		}
		key[2] = `exposes'
		agents.put(key,exposed)
	}
}

void sir_nw::progress(real scalar id, real scalar t)
{
    real rowvector key1, key2, key3
	
	isid(id)
	istime(t)
	
	key1 = id, `state',t
	key2 = id, `dur',t
	if (agents.get(key1, "last") == `infectious'  & agents.get(key2, "last") > mindur) {
		if (runiform(1,1) < pr_heal) {
			agents.put(key1, `removed')
			agents.put(key2,1)
			key3 = id, `exposes',t
			//agents.put(key3,NULL)
		}
	}
    if (agents.get(key1, "last") == `removed' & runiform(1,1) < pr_loss) {
		agents.put(key1, `susceptible')
		agents.put(key2,1)
	}
}

void sir_nw::run()
{
	real scalar i, t
	setup()

	dots(0)
	dots(1)

	for (t=2; t<=tdim; t++) {
		dots(t)
		for(i=1; i<=N ; i++) {
			step(i,t)
		}
	}
}

void sir_nw::dots(real scalar t) 
{
	if (t==0) {
		printf("----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5\n")
	}
	else if (mod(t, 50) == 0) {
		printf(". %9.0f\n", t)
	}
	else {
		printf(".")
	}
	displayflush()
}
end

 
 
-------------------------------------------------------------------------------

<< index >>

-------------------------------------------------------------------------------