Skip to content

Instantly share code, notes, and snippets.

@anttilipp
Created December 19, 2021 09:00
Show Gist options
  • Save anttilipp/137e02c6c862a4d3e203824beecf083b to your computer and use it in GitHub Desktop.
Save anttilipp/137e02c6c862a4d3e203824beecf083b to your computer and use it in GitHub Desktop.
Plots COVID risk graph, probability of infected person present given infection rate in population and number of persons present.
# Antti Lipponen, 18 December 2021
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import zoom
import seaborn as sns
from matplotlib import rcParams
from matplotlib.colors import Normalize
import matplotlib as mpl
from matplotlib.colors import ListedColormap
import matplotlib.patheffects as PathEffects
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Lato']
rcParams['hatch.linewidth'] = 8
rcParams['hatch.color'] = '#00000020'
def compute_prob(incidence, N):
pHasCovid = incidence / 100_000
pNoCovid = 1 - pHasCovid
pNoCovidForGroup = pNoCovid**N
pCovidForGroup = 1 - pNoCovidForGroup
return pCovidForGroup * 100
def computeisoline(p, N):
pCovidForGroup = p / 100
pNoCovidForGroup = 1 - pCovidForGroup
pNoCovid = pNoCovidForGroup**(1 / N)
pHasCovid = 1 - pNoCovid
incidence = pHasCovid * 100_000
return incidence
xlim = [1, 250]
ylim = [50, 10000]
Xlocs = np.linspace(1, 250, 100).astype(int)
Ylocs = np.logspace(np.log10(40), np.log10(10000), 100).astype(int)
Xgrid, Ygrid = np.meshgrid(Xlocs, Ylocs)
Zgrid = np.nan * np.ones_like(Xgrid)
for ii in range(len(Xlocs)):
for jj in range(len(Ylocs)):
Zgrid[jj, ii] = compute_prob(Ylocs[jj], Xlocs[ii])
Xgrid = zoom(Xgrid, 2)
Ygrid = zoom(Ygrid, 2)
Zgrid = zoom(Zgrid, 2)
colors = sns.color_palette("rocket_r", n_colors=20)
for ii in range(4, 20, 2):
meancolor = np.array(colors[ii:ii + 2]).mean(axis=0)
colors[ii] = meancolor
colors[ii + 1] = meancolor
cmap = ListedColormap(colors.as_hex())
norm = Normalize(vmin=0, vmax=100)
levels = np.linspace(0, 100, 500) + 0.01
fig = plt.figure(figsize=(3, 4), dpi=100)
ax = fig.add_axes([0.2, 0.26, 0.75, 0.71])
ax.contourf(Xgrid, Ygrid, Zgrid, levels=levels, cmap=cmap)
for p in [5, 10, 15, 20, 30, 40, 50, 60, 70, 80, 90]:
Nxx = np.linspace(xlim[0], xlim[1], 250)
yy = computeisoline(p, Nxx)
ax.plot(Nxx, yy, 'w-', linewidth=0.3, zorder=99)
ax.grid(True, alpha=0.9, c='#a0a0a0', linestyle='--', linewidth=0.2)
p = 50
xx = np.linspace(6.5, 255, 260)
ax.fill_between(xx, y1=computeisoline(p, xx), y2=11000, hatch='/', zorder=99, fc='#dedede', clip_on=False, alpha=0.25)
ax.plot([6.5 + 0.25, 6.5, 255, 255, 250], [10000, 11000, 11000, computeisoline(p, xx)[-1], computeisoline(p, xx)[-1] + 5], linewidth=0.5, c='#000000', alpha=0.2, clip_on=False)
txt = ax.text(264, 80, 'More likely to have a person with COVID-19 present than none with COVID-19.', fontsize=5, c='#303030', ha='right', va='bottom', alpha=1.0, rotation=-90, zorder=100, clip_on=False)
ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(ylim[0], ylim[1])
ax.set_yscale('log')
ax.set_xlabel('Number of persons present')
ax.set_ylabel('Infected persons / 100,000 persons in the region', labelpad=-1)
ax.set_yticks([100, 1000, 10000])
ax.set_yticklabels(['100', '1,000', '10,000'], rotation=65)
ax.set_xticks([2, 50, 100, 150, 200, 250])
ax.set_xticklabels(['2', '50', '100', '150', '200', '250'], rotation=0)
axcb = fig.add_axes([0.05, 0.11, 0.9, 0.02])
cb = fig.colorbar(mpl.cm.ScalarMappable(norm=norm, cmap=cmap), cax=axcb, orientation='horizontal')
cb.set_label(label='Probability for having at least 1 infected person present (%)', fontsize=7.5, labelpad=1.5)
cb.ax.tick_params(labelsize=8)
ax.text(0.01, 0.005, 'Computed with basic probability theory. Assuming each person is equally likely to be infected.', fontsize=4.5, transform=fig.transFigure, ha='left', va='bottom', alpha=0.35)
ax.text(0.99, 0.005, '@anttilip', fontsize=4.5, transform=fig.transFigure, ha='right', va='bottom', alpha=0.25)
txt = ax.text(35, 53.5, '<5%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(110, 51.5, '5-10%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100, rotation=-30)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(220, 49.5, '10-15%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100, rotation=-13.5)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(225, 69, '15-20%', fontsize=8, c='#303030', ha='center', va='bottom', alpha=0.9, zorder=100, rotation=-13.9)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 123, '20-30%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-13.9)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 188, '30-40%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 260, '40-50%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14.4)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 350, '50-60%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14.5)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 460, '60-70%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-14.5)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 600, '70-80%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-15)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 835, '80-90%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100, rotation=-16)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
txt = ax.text(247.5, 2850, '>90%', fontsize=8, c='#303030', ha='right', va='center', alpha=0.9, zorder=100)
txt.set_path_effects([PathEffects.withStroke(linewidth=1, foreground='w')])
plt.savefig("graph.png", dpi=900)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment