Skip to content

Instantly share code, notes, and snippets.

@stedn

stedn/.block

Last active March 26, 2020 07:01
Show Gist options
  • Star 1 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save stedn/b0acdeba81751c5f8e4ccaa74ccb09a9 to your computer and use it in GitHub Desktop.
Save stedn/b0acdeba81751c5f8e4ccaa74ccb09a9 to your computer and use it in GitHub Desktop.
Download ZIP
Simple Model of Infection Spreading
Raw
.block
license: cc-by-sa-4.0
height: 675
Raw
README.md

This visualizes a highly simplified model of the early stages of the SARS-CoV-2 virus epidemic and how likely it is to spread to the SF bay area.

You can modify the parameters while the simulation is running. Of most interest to me was what happened when you drop the rate of spread from Wuhan to SF (set mu=0). There's a time window after which stopping flights is too late. Read the model description here.

The javascript code solves the following equation using the forward euler method with a quirky little hack to load people onto planes as part of the animation.

equation

If you want to see this problem analyzed in more depth, check out this project on github

Raw
index.html
<!DOCTYPE html>
<meta charset="utf-8">
<head>
<title>Wuhan->SF Virus Simulation</title>
</head>
<style>
#wrapper {
font-family: "Helvetica Neue", Helvetica, sans-serif;
font-size: 14px;
color: #333;
background: '#333';
width:700px;
margin:auto;
margin-top:10px;
border-radius: 3px;
}
#content{
background: black;
}
#play-button {
top: 140px;
left: 50px;
background: #f08080;
padding-right: 26px;
border-radius: 3px;
border: none;
color: white;
margin: 5px;
padding: 0 12px;
width: 60px;
cursor: pointer;
height: 30px;
}
#reset-button {
top: 140px;
left: 50px;
background: #696969;
padding-right: 26px;
border-radius: 3px;
border: none;
color: white;
margin: 5px;
padding: 0 12px;
width: 60px;
cursor: pointer;
height: 30px;
}
#play-button:hover {
background-color: #696969;
}
</style>
<body>
<div id="wrapper">
<div id="content">
<canvas width="700" height="400"></canvas>
</div>
<div id="vis">
<button id="play-button">Start</button> <button id="reset-button">Reset</button>
</div>
<div id="frm">Parameters:
<label for="mu"
style="display: inline-block; width: 80; text-align: right">
flight rate = <span id="mu-value"></span>
</label>
<input type="number" min="0" max="0.1" step="0.001" id="mu">
<label for="gamma"
style="display: inline-block; width: 80; text-align: right">
infection rate = <span id="gamma-value"></span>
</label>
<input type="number" min="0" max="10" step="0.1" id="gamma">
<label for="kappa"
style="display: inline-block; width: 80; text-align: right">
quarantine = <span id="kappa-value"></span>
</label>
<input type="number" min="0" max="100" step="1" id="kappa">
<label for="n_init"
style="display: inline-block; width: 80; text-align: right">
init pop = <span id="n_init-value"></span>
</label>
<input type="number" min="0" max="200" step="10" id="n_init">
</div>
<div id="my_dataviz"></div>
</div>
<script src="https://cdnjs.cloudflare.com/ajax/libs/d3/4.2.2/d3.min.js"></script>
<script>
var margin = {top: 10, right: 30, bottom: 30, left: 60},
width = 630,
height = 100;
duration = 10;
sf_color = 'tomato';
wu_color = 'dodgerblue';
// append the svg object to the body of the page
var svg = d3.select("#my_dataviz")
.append("svg")
.attr("width", width + margin.left + margin.right)
.attr("height", height + margin.top + margin.bottom)
.append("g")
.attr("transform",
"translate(" + (margin.left/2) + "," + (margin.top/2) + ")");
// model params
gamma = 1.3;
kappa = 10;
mu = 0.01;
n_wu_start = 100;
d3.select("#mu").property('value', mu);
d3.select("#gamma").property('value', gamma);
d3.select("#kappa").property('value', kappa);
d3.select("#n_init").property('value', n_wu_start);
// when the input range changes update value
d3.select("#mu").on("input", function() {
mu = +this.value;
});
d3.select("#gamma").on("input", function() {
gamma = +this.value;
});
d3.select("#kappa").on("input", function() {
kappa = +this.value;
});
d3.select("#n_init").on("input", function() {
n_wu_start = +this.value;
});
function gaussianRand() {
var rand = 0;
var cnt = 6;
for (var i = 0; i < cnt; i += 1) {
rand += Math.random();
}
return (rand / cnt) -0.5;
}
function randomCircle(rad){
pt_angle = Math.random() * 2 * Math.PI;
pt_radius_sq = Math.sqrt(Math.abs(gaussianRand()) * rad * rad);
pt_x = pt_radius_sq * Math.cos(pt_angle);
pt_y = pt_radius_sq * Math.sin(pt_angle);
return [pt_x, pt_y];
}
var geojson = {}
var context = d3.select('#content canvas')
.node()
.getContext('2d');
var projection = d3.geoConicEquidistant()
.scale(280)
.center([-179,40])
.translate([40,-50])
.rotate([180, 15, 10]);
var geoGenerator = d3.geoPath()
.projection(projection)
.pointRadius(2)
.context(context);
var wuhanLonLat = [114.3, 30.6];
var sfLonLat = [-122.4, 37.8];
var geoInterpolator = d3.geoInterpolate(wuhanLonLat, sfLonLat);
dt = 0.01;
log_inc = 0.1
// rendering params and vars
max_n = 1000;
rad_wu = 8;
rad_sf = 4;
dl = 0.05
// function to put us back to square 1
function reset() {
t = 0;
next_t_log = 0;
log = []
n_wu = n_wu_start;
n_sf = 0;
n_tr = 0;
i_wu = [];
i_wu_n = 0;
i_sf = [];
i_sf_n = 0;
i_tr = [];
i_tr_n = 0;
x_domain_size = 10;
y_sf_domain_size = n_wu;
y_wu_domain_size = n_wu*4;
resetting = true;
}
firsttime = true;
// simulation and viz function
function update() {
// compute differential euations for isolated populations
d_wu = (gamma * n_wu - kappa) * dt;
d_sf = (gamma * n_sf - kappa) * dt;
n_wu = n_wu + d_wu;
n_sf = n_sf + d_sf;
d_tr = (mu * n_wu) * dt
n_tr = n_tr + d_tr
// transition dots out of Wuhan and into "in transit" list
while(n_tr > 1){
n_tr = n_tr - 1;
n_wu = n_wu - 1;
i_tr.push(0);
i_tr_n = i_tr_n + 1;
}
// for dots that have crossed, remove from "in transit" and add to n_sf
for (i = 0; i < i_tr_n; i++) {
i_tr[i] = i_tr[i] + dl + dl*gaussianRand();
if(i_tr[i]>=1){
i_tr.splice(i, 1);
i_tr_n = i_tr_n - 1;
n_sf = n_sf + 1;
}
}
// hard limits to non-negative
if (n_wu<0) {
n_wu = 0;
}
if (n_sf<0) {
n_sf = 0;
}
// only need to vizualize up to some max_n
n_wu_viz = n_wu;
if (n_wu_viz>max_n) {
n_wu_viz = max_n;
}
n_sf_viz = n_sf;
if (n_sf_viz>max_n) {
n_sf_viz = max_n;
}
// add/remove people to the lists to get the right number of people
while (i_wu_n > Math.ceil(n_wu_viz)){
i_wu.splice(Math.floor(Math.random()*i_wu.length), 1);
i_wu_n = i_wu_n - 1;
}
while(i_wu_n < Math.floor(n_wu_viz)){
pt = randomCircle(rad_wu)
i_wu.push([wuhanLonLat[0]+pt[0], wuhanLonLat[1]+pt[1]])
i_wu_n = i_wu_n + 1
}
while (i_sf_n > Math.ceil(n_sf_viz)){
i_sf.splice(Math.floor(Math.random()*i_sf.length), 1);
i_sf_n = i_sf_n - 1;
}
while(i_sf_n < Math.floor(n_sf_viz)){
pt = randomCircle(rad_sf)
i_sf.push([sfLonLat[0]+pt[0], sfLonLat[1]+pt[1]])
i_sf_n = i_sf_n + 1
}
//DRAW THINGS
context.clearRect(0, 0, 800, 600);
// country map
context.strokeStyle = '#666';
context.lineWidth = 1;
context.globalAlpha = 1;
context.beginPath();
geoGenerator({type: 'FeatureCollection', features: geojson.features})
context.stroke();
// graticule (lat-lon lines)
context.lineWidth = 0.5;
context.globalAlpha = 0.5;
var graticule = d3.geoGraticule();
context.beginPath();
context.strokeStyle = '#ccc';
geoGenerator(graticule());
context.stroke();
// Wuhan Bubble
context.lineWidth = 2;
context.globalAlpha = 1;
context.strokeStyle = wu_color;
context.beginPath()
var circle = d3.geoCircle()
.center(wuhanLonLat)
.radius(0.75*rad_wu);
geoGenerator(circle());
context.stroke()
// SF Bubble
context.strokeStyle = sf_color;
context.beginPath()
var circle = d3.geoCircle()
.center(sfLonLat)
.radius(0.75*rad_sf);
geoGenerator(circle());
context.stroke()
// Wuhan - SF flight line
context.lineWidth = 2;
context.globalAlpha = 0.5;
context.beginPath();
context.strokeStyle = 'grey';
geoGenerator({type: 'Feature', geometry: {type: 'LineString', coordinates: [wuhanLonLat, sfLonLat]}});
context.stroke();
// dots in Wuhan
context.globalAlpha = 0.8;
context.fillStyle = wu_color;
for (i = 0; i < i_wu_n; i++) {
context.beginPath();
geoGenerator({type: 'Feature', geometry: {type: 'Point', coordinates: i_wu[i], radius: 100}});
context.fill();
}
// dots in flight
context.globalAlpha = 1;
for (i = 0; i < i_tr_n; i++) {
context.beginPath();
geoGenerator({type: 'Feature', geometry: {type: 'Point', coordinates: geoInterpolator(i_tr[i])}});
context.fill();
}
// dots in SF
context.globalAlpha = 0.75;
context.fillStyle = sf_color;
for (i = 0; i < i_sf_n; i++) {
context.beginPath();
geoGenerator({type: 'Feature', geometry: {type: 'Point', coordinates: i_sf[i]}});
context.fill();
}
// GRAPHING
// first time set up axes and lines
if(firsttime){
x = d3.scaleLinear().domain([0, x_domain_size]).range([0, width]);
x_axis = d3.axisBottom(x);
x_axis_svg = svg.append("g")
.attr("class", "axis axis--x")
.attr("transform", "translate(0," + height + ")");
x_axis_svg.call(x_axis);
y_wu = d3.scaleLinear().domain([0, y_wu_domain_size]).range([height, 0]);
y_axis_wu = d3.axisLeft(y_wu);
y_axis_svg_wu = svg.append("g")
.attr("class", "axis axis--y")
.style("fill", wu_color);
y_axis_svg_wu.call(y_axis_wu);
y_sf = d3.scaleLinear().domain([0, y_sf_domain_size]).range([height, 0]);
y_axis_sf = d3.axisRight(y_sf);
y_axis_svg_sf = svg.append("g")
.attr("class", "axis axis--y")
.style("fill", sf_color)
.attr("transform", "translate(" + width + " ,0)");
y_axis_svg_sf.call(y_axis_sf);
wu_valueline = d3.line()
.x(function(d) { return x(d[0]); })
.y(function(d) { return y_wu(d[1]); });
wu_data_line = svg.append("path")
.attr("fill", "none")
.attr("stroke", wu_color)
.attr("stroke-width", 1.5)
.attr("class", "line")
.attr("d", wu_valueline(log));
sf_valueline = d3.line()
.x(function(d) { return x(d[0]); })
.y(function(d) { return y_sf(d[2]); });
sf_data_line = svg.append("path")
.attr("fill", "none")
.attr("stroke", sf_color)
.attr("stroke-width", 1.5)
.attr("class", "line")
.attr("d", wu_valueline(log));
firsttime = false;
}
// if reset has occurred make sure to reset axes
if (resetting) {
x.domain([0, x_domain_size])
x_axis_svg.transition().duration(duration).ease(d3.easeLinear, 2).call(x_axis);
y_wu.domain([0, y_wu_domain_size])
y_axis_svg_wu.transition().duration(duration).ease(d3.easeLinear, 2).call(y_axis_wu);
y_sf.domain([0, y_sf_domain_size])
y_axis_svg_sf.transition().duration(duration).ease(d3.easeLinear, 2).call(y_axis_sf);
resetting = false;
}
// only update graph every so often
if (t >= next_t_log) {
next_t_log = t + log_inc;
log.push([t,n_wu,n_sf]);
// might need to rescale axes
if (t > x_domain_size){
x_domain_size = t
x.domain([0, x_domain_size])
x_axis_svg.transition().duration(duration).ease(d3.easeLinear, 2).call(x_axis);
}
if (n_wu > y_wu_domain_size){
y_wu_domain_size = n_wu*1.5
y_wu.domain([0, y_wu_domain_size])
y_axis_svg_wu.transition().duration(duration).ease(d3.easeLinear, 2).call(y_axis_wu);
}
if (n_sf > y_sf_domain_size){
y_sf_domain_size = n_sf*1.5
y_sf.domain([0, y_sf_domain_size])
y_axis_svg_sf.transition().duration(duration).ease(d3.easeLinear, 2).call(y_axis_sf);
}
// update the lines
wu_data_line
.transition()
.duration(duration)
.attr("d", wu_valueline(log));
sf_data_line
.transition()
.duration(duration)
.attr("d", sf_valueline(log));
}
// increment time (only needed for logging)
t = t + dt;
// if the total infeccted passes 100k stop things
if(n_wu + n_sf > 100000){
clearInterval(interval);
d3.select("#play-button").text("Start")
alert('Reached 100k infected. Violates early phase assumptions. stopping.')
}
}
// REQUEST DATA
d3.json('world110m2.json', function(err, json) {
geojson = json;
reset();
update();
interval = null;
// add a play/pause and reset buttons
// https://bl.ocks.org/officeofjane/47d2b0bfeecfcb41d2212d06d095c763
var playButton = d3.select("#play-button");
playButton
.on("click", function() {
var button = d3.select(this);
if (button.text() == "Pause") {
clearInterval(interval);
button.text("Start");
} else {
interval = setInterval(update, 50);
button.text("Pause");
}
})
var resetButton = d3.select("#reset-button");
resetButton
.on("click", function() {
var button = d3.select(this);
clearInterval(interval);
reset();
update();
d3.select("#play-button").text("Start")
})
})
// resources that I based parts of this off of:
// - https://www.d3-graph-gallery.com/graph/line_basic.html
// - http://bl.ocks.org/benjchristensen/2579599
// - https://bl.ocks.org/d3noob/7030f35b72de721622b8
// - https://www.d3indepth.com/geographic/
</script>
</body>
</html>
Display the source blob
Display the rendered blob
Raw
world110m2.json
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment