The code format is rendering with *'s for italics. Here is hopefully a more
legible paste:
load 'viewmat'
coinsert'jgl2'
wd 'pc Covid-19-Sim closeok' NB. parent control (window) named 'w0'
wd 'minwh 500 500; cc g0 isidraw;' NB. add an 'isidraw' child control
named 'g0'
wd 'pshow; pmove 40 510 0 0' NB. show the window at the given
coordinates.
wd 'sm focus term' NB. session manager: bring terminal to
front
wd 'psel Covid-19-Sim; ptop' NB. bring our window to front
vmcc =: viewmatcc_jviewmat_ NB. viewmat to a child control
step =: render @ update NB. each step, we'll call those two in
sequence
RGB=: (#:i.8) { 0 255
DSEIR =: (0 7 6 4 2 { RGB)
update =: verb define
im =: 12 12 $ infect_display starting_pop
)
render =: verb define
DSEIR vmcc im;'g0'
glpaint''
)
sys_timer_z_ =: step_base_
wd 'timer 0'
ppl =: 144 NB. How many people are in our simulation?
Beta =: 0.106271 NB. The rate at which susceptible people become
exposed. RPH - dependent on R0 = 2.2 3/25/2020
Sigma =: 0.066967 NB. The rate at which exposed people become
infectious. RPH - (1-sigma)^10=0.5, 10 is incubation period.
Gamma =: 0.056126 NB. The rate at which infectious people recover. RPH -
adjusted to 12 days, so (1-Gamma)^12=0.5.
Xi =: 0.02284 NB. The rate at which recovered people lose immunity
and become susceptible. RPH - 30 days
NB. is partially based upon man on Diamond Princess who presented
symptoms around 30 days otherwise made up.
R =: 2.2 NB. How many people will 1 person infect? RPH - based
on latest CDC estimate 3/25/2020
CFR =: 0.015 NB. Case Fatality Rate - RPH based upon 0.25% to 3%
from CDC, so 1.5% is higher than average.
D =: 0.0011619 NB. Death rate per day while symptomatic - RPH based
upon 13 days to death as seen and per CDC.
ClosenessStdDev =: 5 % 3
ClosenessMean =: 1
rand_norm =: 3 : '(2&o. 2p1 * ? y) * %: _2 * ^. ? y'
closeness =: 4 %~ 2 %~ (+ |:) ClosenessMean + ClosenessStdDev * rand_norm
(ppl,ppl) $ 0
risk =: closeness - closeness * =i.ppl
NB. risk =: (ppl,ppl)$0
NB. 0 = Dead
NB. 1 = Susceptible
NB. 2 = Exposed
NB. 3 = Infectious
NB. 4 = Recovered
starting_pop =: 2 , (ppl - 1) $ 1
infect =: 3 : 0
can_spread =. (1&< *. 4&>) y
susceptible =. y = 1
exposed =. y = 2
infectious =. y = 3
recovered =. y = 4
infectiousness =. can_spread ,./ . * risk
susceptible_to_exposed =. (Beta * susceptible) * (+/ can_spread) %~ +/ |:
susceptible * infectiousness
exposed_to_infectious =. Sigma * exposed
infectious_to_recovered =. Gamma * infectious
infectious_to_dead =. D * infectious
recovered_to_susceptible =. Xi * recovered
1 (I.recovered_to_susceptible>?ppl$0)} 4
(I.infectious_to_recovered>?ppl$0)} 0 (I.infectious_to_dead>?ppl$0)} 3
(I.exposed_to_infectious>?ppl$0)} 2 (I.susceptible_to_exposed>?ppl$0)} y
)
infect_display =: 3 : 0
out =. infect y
smoutput out
smoutput
('Susceptible';'Exposed';'Infectious';'Recovered';'Dead'),:(+/1=out);(+/2=out);(+/3=out);(+/4=out);(+/0=out)
out
)
On Fri, 27 Mar 2020 at 14:24, Robert Herman <[email protected]> wrote:
> I found this article:
>
> Modeling the COVID-19 Outbreak with J
> <https://datakinds.github.io/2020/03/15/modeling-the-coronavirus-outbreak-with-j>
>
> I took the author's code, which is a modified SEIR model (The_SEIR_model
> <https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology#The_SEIR_model>),
> and updated some of the data and recalculated the variables. I then
> pilfered some of Michal's (Tangentstorm, you rock!) code for animating
> things in J from his talks. This is helping me to push my learning of J,
> and keep it interesting.
>
> Two disclaimers: The author of the linked article doesn't claim to be an
> epidemiologist, and I am only doing this to learn J, and cannot comment on
> any of the validity of the algorithm or simulations it produces.
>
> Now, on to my two issues (full code at end):
>
> 1. I know I have to figure out a way to have the* infect_display^:144
> starting population* call feed the *update* loop, otherwise it is trying
> to run all 144 days of simulation in one call to update at a frequency set
> by *wd 'timer 600'*. I am working on this, and realize I somehow have to
> run it outside of update, and then feed update "a day at a time". The goal
> is to see the same person (grid location) evolve over the simulation. I
> know this is simple, but i haven't gotten there yet. Any tips will be
> appreciated!
>
> 2. When I had initially set up the code to run 25 people instead of 144,
> the color index I had setup with the DSEIR variable exhibit two strange
> glitches. First, squares that were supposed to be 1=White, were rendering
> as 0=Black=Dead, and sometimes colors appeared "blended". An yellowish,
> off-white for yellow, and at times bright yellow. When I went to 144 people
> and it is trying to run in update, the grid appears to be correctly
> indexing and not moving askew from frame to frame, but the colors are
> inconsistent.
>
> BTW, thanks to Linda's article in Journal of J - Vol.4, No.1 for the RGB
> list to then pull my desired colors from in the DSEIR variable (Dead,
> Susceptible, Exposed, Infectious, Recovered):
>
> *RGB=: (#:i.8) { 0 255*
>
>
> *DSEIR =: (0 7 6 4 2 { RGB)*
>
> Here is my integration of the two sources of code to try and animate the
> outbreak per the article, and thanks!:
>
> *load 'viewmat'*
>
> *coinsert'jgl2'*
>
>
> *wd 'pc Covid-19-Sim closeok' NB. parent control (window) named 'w0'*
>
> *wd 'minwh 500 500; cc g0 isidraw;' NB. add an 'isidraw' child control
> named 'g0'*
>
> *wd 'pshow; pmove 40 510 0 0' NB. show the window at the given
> coordinates.*
>
> *wd 'sm focus term' NB. session manager: bring terminal to front*
>
> *wd 'psel Covid-19-Sim; ptop' NB. bring our window to front*
>
>
> *vmcc =: viewmatcc_jviewmat_ NB. viewmat to a child control*
>
>
> *step =: render @ update NB. each step, we'll call those two in sequence*
>
>
> *RGB=: (#:i.8) { 0 255*
>
>
> *DSEIR =: (0 7 6 4 2 { RGB)*
>
>
> *update =: verb define*
>
> * im =: 12 12 $ infect_display^:144 starting_pop*
>
> *)*
>
>
> *render =: verb define*
>
> * DSEIR vmcc im;'g0'*
>
> * glpaint''*
>
> *)*
>
>
> *sys_timer_z_ =: step_base_*
>
> *wd 'timer 600'*
>
>
> *ppl =: 144 NB. How many people are in our simulation?*
>
>
> *Beta =: 0.106271 NB. The rate at which susceptible people become exposed.
> RPH - dependent on R0 = 2.2 3/25/2020*
>
> *Sigma =: 0.066967 NB. The rate at which exposed people become infectious.
> RPH - (1-sigma)^10=0.5, 10 is incubation period.*
>
> *Gamma =: 0.056126 NB. The rate at which infectious people recover. RPH -
> adjusted to 12 days, so (1-Gamma)^12=0.5.*
>
> *Xi =: 0.02284 NB. The rate at which recovered people lose immunity and
> become susceptible. RPH - 30 days*
>
> * NB. is partially based upon man on Diamond Princess who presented
> symptoms around 30 days otherwise made up.*
>
> *R =: 2.2 NB. How many people will 1 person infect? RPH - based on latest
> CDC estimate 3/25/2020*
>
> *CFR =: 0.015 NB. Case Fatality Rate - RPH based upon 0.25% to 3% from
> CDC, so 1.5% is higher than average.*
>
> *D =: 0.0011619 NB. Death rate per day while symptomatic - RPH based upon
> 13 days to death as seen and per CDC.*
>
>
> *ClosenessStdDev =: 5 % 3*
>
> *ClosenessMean =: 1*
>
>
> *rand_norm =: 3 : '(2&o. 2p1 * ? y) * %: _2 * ^. ? y'*
>
> *closeness =: 4 %~ 2 %~ (+ |:) ClosenessMean + ClosenessStdDev * rand_norm
> (ppl,ppl) $ 0*
>
> *risk =: closeness - closeness * =i.ppl *
>
> *NB. risk =: (ppl,ppl)$0*
>
>
> *NB. 0 = Dead*
>
> *NB. 1 = Susceptible*
>
> *NB. 2 = Exposed*
>
> *NB. 3 = Infectious*
>
> *NB. 4 = Recovered*
>
> *starting_pop =: 2 , (ppl - 1) $ 1*
>
>
> *infect =: 3 : 0*
>
> * can_spread =. (1&< *. 4&>) y*
>
> * susceptible =. y = 1*
>
> * exposed =. y = 2*
>
> * infectious =. y = 3*
>
> * recovered =. y = 4*
>
>
> * infectiousness =. can_spread ,./ . * risk*
>
>
> * susceptible_to_exposed =. (Beta * susceptible) * (+/ can_spread) %~ +/
> |: susceptible * infectiousness*
>
> * exposed_to_infectious =. Sigma * exposed*
>
> * infectious_to_recovered =. Gamma * infectious*
>
> * infectious_to_dead =. D * infectious*
>
> * recovered_to_susceptible =. Xi * recovered*
>
> * 1 (I.recovered_to_susceptible>?ppl$0)} 4
> (I.infectious_to_recovered>?ppl$0)} 0 (I.infectious_to_dead>?ppl$0)} 3
> (I.exposed_to_infectious>?ppl$0)} 2 (I.susceptible_to_exposed>?ppl$0)} y*
>
> *)*
>
>
> *infect_display =: 3 : 0*
>
> * out =. infect y*
>
> * smoutput out*
>
> * smoutput
> ('Susceptible';'Exposed';'Infectious';'Recovered';'Dead'),:(+/1=out);(+/2=out);(+/3=out);(+/4=out);(+/0=out)*
>
> * out*
>
> *)*
>
>
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
