API tools faq
paste
creamygoat

Herd Immunity Threshold Graph Plotter

creamygoat
Sep 16th, 2021 (edited)
1,398
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 10.39 KB | None | 0 0
raw download clone embed print report
  1. #!/usr/bin/env python3
  2.  
  3.  
  4. import numpy as np
  5. import matplotlib
  6. import matplotlib.pyplot as plt
  7. import matplotlib.ticker as ticker
  8. from matplotlib.path import Path
  9. import matplotlib.lines as lines
  10. import matplotlib.patches as patches
  11. import matplotlib.colors as colors
  12. import matplotlib.hatch as mhatch
  13.  
  14.  
  15. fathatch_path_vertices = [
  16.   [-0.5, -0.45],
  17.   [-0.5, -0.5],
  18.   [-0.375, -0.5],
  19.   [-0.5, -0.375],
  20.   [-0.5, -0.45],
  21.  
  22.   [-0.5, 0],
  23.   [-0.5, -0.125],
  24.   [-0.125, -0.5],
  25.   [0.125, -0.5],
  26.   [-0.5, 0.125],
  27.   [-0.5, 0],
  28.  
  29.   [-0.5, 0.45],
  30.   [-0.5, 0.375],
  31.   [0.375, -0.5],
  32.   [0.5, -0.5],
  33.   [0.5, -0.375],
  34.   [-0.375, 0.5],
  35.   [-0.5, 0.5],
  36.   [-0.5, 0.45],
  37.  
  38.   [0, 0.5],
  39.   [-0.125, 0.5],
  40.   [0.5, -0.125],
  41.   [0.5, 0.125],
  42.   [0.125, 0.5],
  43.   [0, 0.5],
  44.  
  45.   [0.45, 0.5],
  46.   [0.375, 0.5],
  47.   [0.5, 0.375],
  48.   [0.5, 0.5],
  49.   [0.45, 0.5],
  50. ]
  51.  
  52. fathatch_path_cmds = [
  53.   Path.MOVETO,
  54.   Path.LINETO,
  55.   Path.LINETO,
  56.   Path.LINETO,
  57.   Path.CLOSEPOLY,
  58.  
  59.   Path.MOVETO,
  60.   Path.LINETO,
  61.   Path.LINETO,
  62.   Path.LINETO,
  63.   Path.LINETO,
  64.   Path.CLOSEPOLY,
  65.  
  66.   Path.MOVETO,
  67.   Path.LINETO,
  68.   Path.LINETO,
  69.   Path.LINETO,
  70.   Path.LINETO,
  71.   Path.LINETO,
  72.   Path.LINETO,
  73.   Path.CLOSEPOLY,
  74.  
  75.   Path.MOVETO,
  76.   Path.LINETO,
  77.   Path.LINETO,
  78.   Path.LINETO,
  79.   Path.LINETO,
  80.   Path.CLOSEPOLY,
  81.  
  82.   Path.MOVETO,
  83.   Path.LINETO,
  84.   Path.LINETO,
  85.   Path.LINETO,
  86.   Path.CLOSEPOLY,
  87. ]
  88.  
  89. fat_left_hatch_path = patches.Path(
  90.   [[x, y] for x, y in fathatch_path_vertices],
  91.   fathatch_path_cmds,
  92. )
  93.  
  94. fat_right_hatch_path = patches.Path(
  95.   [[-x, y] for x, y in fathatch_path_vertices],
  96.   fathatch_path_cmds,
  97. )
  98.  
  99.  
  100. class FatLHatch(mhatch.Shapes):
  101.     filled = True
  102.     size = 1.0
  103.     path = fat_left_hatch_path
  104.     def __init__(self, hatch, density):
  105.       self.num_rows = 2 * (hatch.count('C')) # * density
  106.       self.shape_vertices = self.path.vertices
  107.       self.shape_codes = self.path.codes
  108.       mhatch.Shapes.__init__(self, hatch, density)
  109.  
  110.  
  111. class FatRHatch(mhatch.Shapes):
  112.     filled = True
  113.     size = 1.0
  114.     path = fat_right_hatch_path
  115.     def __init__(self, hatch, density):
  116.       self.num_rows = 2 * (hatch.count('c')) # * density
  117.       self.shape_vertices = self.path.vertices
  118.       self.shape_codes = self.path.codes
  119.       mhatch.Shapes.__init__(self, hatch, density)
  120.  
  121.  
  122. def lin_fn(x):
  123.   # Linear reconstruction filter
  124.   return (x + 1 if x < 0 else 1 - x) if -1 < x < 1 else 0
  125.  
  126.  
  127. def lanczos(x, a):
  128.   # Lanczos reconstruction filter with size parameter a, a positive integer
  129.   #
  130.   # a = 1 is useless. Larger values increasingly emulate the normalised sinc
  131.   # function. (Lanczos is just the normalised sinc function windowed with the
  132.   # horizontally scaled central lobe of another normalised sinc function.)
  133.   return (np.sinc(x) * np.sinc(x / a)) if -a < x < a else 0
  134.  
  135.  
  136. def lanczos_fn(a):
  137.   # Return a 1-argt lanczos function tailored with a specific size parameter.
  138.   return lambda x: lanczos(x, a)
  139.  
  140.  
  141. def resample(values, new_x_posns, filter_fn, support):
  142.   result = []
  143.   n = len(values)
  144.   a = max(0, -support[0])
  145.   b = max(0, support[1])
  146.   for x1 in new_x_posns:
  147.     x = n * x1
  148.     rng = range(
  149.       int(np.floor(x)) + support[0] + 1,
  150.       int(np.floor(x)) + support[1] + 1
  151.     )
  152.     s = sum(values[np.clip(i, 0, n - 1)] * filter_fn(x - i) for i in rng)
  153.     result.append(s)
  154.   return result
  155.  
  156.  
  157. def resample_lanczos2(values, new_x_posns):
  158.   result = resample(values, new_x_posns, lanczos_fn(2), (-2, 2))
  159.   return result
  160.  
  161.  
  162. def hit_fn(r0, ver):
  163.   # Re = R_0 (1 - X.E)
  164.   return (1 - 1 / r0) / ver
  165.  
  166.  
  167. def soft_edge_grad_image(colour):
  168.   c = colors.to_rgba(colour)
  169.   ramp = [0.2, 0.6, 0.7, 0.8, 0.9, 0.95]
  170.   x = ramp + [1.0] + list(reversed(ramp))
  171.   n = len(x)
  172.   c3 = np.tile(np.atleast_2d(c[:3]), (n, 1))
  173.   alpha = np.atleast_2d(np.array(x)).T * c[3]
  174.   im = np.expand_dims(np.hstack([c3, alpha]), 0)
  175.   return im
  176.  
  177.  
  178. def main():
  179.  
  180.   r0_max = 20
  181.  
  182.   mhatch._hatch_types.append(FatLHatch)
  183.   mhatch._hatch_types.append(FatRHatch)
  184.  
  185.   fig = plt.figure(figsize=(5.5, 7.0))
  186.   ax = plt.axes()
  187.   ax.xaxis.set_major_locator(ticker.MultipleLocator(1))
  188.   ax.yaxis.set_major_locator(ticker.MultipleLocator(10))
  189.   ax.yaxis.set_minor_locator(ticker.MultipleLocator(5))
  190.   plt.xlim(0, r0_max)
  191.   plt.ylim(0, 105)
  192.  
  193.   # ~ P = [
  194.     # ~ (1, 0.0),
  195.     # ~ (1.25, 0.20),
  196.     # ~ (1.6, 0.375),
  197.     # ~ (2.0, 0.50),
  198.     # ~ (2.5, 0.60),
  199.     # ~ (4.0, 0.75),
  200.     # ~ (5.0, 0.80),
  201.     # ~ (8.0, 0.875),
  202.     # ~ (10.0, 0.90),
  203.     # ~ (12.5, 0.92),
  204.     # ~ (16.0, 0.9375),
  205.     # ~ (20.0, 0.95),
  206.   # ~ ]
  207.  
  208.   legend_patches = []
  209.   labels = []
  210.  
  211.   num_subdivs = 10
  212.   x = np.linspace(0, r0_max, r0_max * num_subdivs, endpoint=True)[1:]
  213.  
  214.   lp100 = dict(
  215.     color = 'black',
  216.     linewidth = 1.5,
  217.   )
  218.   lp10 = dict(
  219.     color = 'black',
  220.     linewidth = 0.6,
  221.   )
  222.   lp5 = dict(
  223.     color = 'black',
  224.     linewidth = 0.3,
  225.     linestyle = "dashed"
  226.   )
  227.  
  228.   for ver_pc in [65, 70, 75, 80, 85, 90, 95, 100]:
  229.     if ver_pc == 100:
  230.       lparams = lp100
  231.     else:
  232.       lparams = lp10 if ver_pc % 10 == 0 else lp5
  233.     label = "{}%".format(ver_pc)
  234.     y = [100 * hit_fn(x_i, ver_pc/100) for x_i in x]
  235.     ax.plot(x, y, label=label, **lparams)
  236.  
  237.   legend_patches.append(lines.Line2D([], [], **lp100))
  238.   labels.append("VE = 100%")
  239.   legend_patches.append(lines.Line2D([], [], **lp10))
  240.   labels.append("VE 10% multiple")
  241.   legend_patches.append(lines.Line2D([], [], **lp5))
  242.   labels.append("VE 5% multiple")
  243.  
  244.   y = [100 * hit_fn(x_i, 1) for x_i in x]
  245.   curve_points = np.vstack([x, y]).T
  246.   upper_clip_path = patches.Polygon(
  247.     np.vstack([curve_points, [[r0_max, 100], [0, 100]]]),
  248.     closed=True, transform=ax.transData,
  249.   )
  250.   lower_clip_path = patches.Polygon(
  251.     np.vstack([curve_points, [[r0_max, 0], [0, 0]]]),
  252.     closed=True, transform=ax.transData,
  253.   )
  254.  
  255.   D = {
  256.     'influenza': dict(
  257.       label = "Influenza",
  258.       r0_span = [1.0, 2.0],
  259.       colour = "red",
  260.     ),
  261.     'sars': dict(
  262.       label = "SARS",
  263.       r0_span = [2.0, 4.0],
  264.       colour = "#00cc55",
  265.     ),
  266.     'polio': dict(
  267.       label = "Polio",
  268.       r0_span = [5.0, 7.0],
  269.       colour = "#0066ee",
  270.     ),
  271.     'chickenpox': dict(
  272.       label = "Chickenpox",
  273.       r0_span = [10.0, 12.0],
  274.       colour = "orange",
  275.     ),
  276.     'measles': dict(
  277.       label = "Measles",
  278.       r0_span = [12.0, 18.0],
  279.       colour = "#ff00cc",
  280.     ),
  281.     'urcovid19': dict(
  282.       label = "Covid-19 (Ancestral)",
  283.       r0_span = [2.3, 3.5],
  284.       colour = "#ffdd00",
  285.       special_hatching = True,
  286.       #hatch = r"//",
  287.       hatch = r"C",
  288.     ),
  289.     'covid19alpha': dict(
  290.       #label = "Covid-19-\N{GREEK SMALL LETTER ALPHA}",
  291.       label = "Covid-19 (Alpha)",
  292.       r0_span = [4.0, 5.0],
  293.       colour = "#ffdd00",
  294.       special_hatching = True,
  295.       #hatch = r"\\",
  296.       hatch = r"CC",
  297.     ),
  298.     'covid19delta': dict(
  299.       label = "Covid-19 (Delta)",
  300.       r0_span = [5.0, 8.0],
  301.       colour = "#ffdd00",
  302.       special_hatching = True,
  303.       #hatch = r"///",
  304.       hatch = r"c",
  305.     ),
  306.   }
  307.  
  308.   # ~ D = {
  309.     # ~ 'lurgy': dict(
  310.       # ~ label = "Lurgy",
  311.       # ~ r0_span = [1.5, 8.0],
  312.       # ~ colour = "red",
  313.     # ~ ),
  314.   # ~ }
  315.  
  316.   for key, drec in D.items():
  317.  
  318.     ver = 1
  319.     band_alpha = 0.33
  320.  
  321.     r0_span = drec['r0_span']
  322.     hit_span = [100 * hit_fn(x_i, ver) for x_i in r0_span]
  323.     colour = drec.get('colour', None)
  324.     label = drec.get("label", None)
  325.     hatch=drec.get("hatch", None)
  326.  
  327.     if drec.get('special_hatching', False):
  328.       # Special Covid-19 hatching
  329.  
  330.       v_xy = [
  331.         (r0_span[0], 0),
  332.         (r0_span[1], 0),
  333.         (r0_span[1], 200),
  334.         (r0_span[0], 200),
  335.       ]
  336.       vpatch = patches.Polygon(
  337.         v_xy,
  338.         facecolor=colour,
  339.         closed=True,
  340.         fill=False,
  341.         edgecolor=colour,
  342.         hatch=hatch,
  343.         label=label,
  344.       )
  345.       vp = ax.add_patch(vpatch)
  346.       vp.set_clip_path(lower_clip_path)
  347.  
  348.       h_xy = [
  349.         (0, hit_span[0]),
  350.         (r0_max, hit_span[0]),
  351.         (r0_max, hit_span[1]),
  352.         (0, hit_span[1]),
  353.       ]
  354.  
  355.       hpatch = patches.Polygon(
  356.         h_xy,
  357.         facecolor=colour,
  358.         closed=True,
  359.         fill=False,
  360.         edgecolor=colour,
  361.         hatch=hatch,
  362.       )
  363.       hp = ax.add_patch(hpatch)
  364.       hp.set_clip_path(upper_clip_path)
  365.  
  366.       lpatch = patches.Patch(
  367.         facecolor=colour,
  368.         fill=False,
  369.         edgecolor=colour,
  370.         hatch=hatch,
  371.         label=label
  372.       )
  373.       legend_patches.append(lpatch)
  374.       labels.append(label)
  375.  
  376.     else:
  377.       # soft-edge bands
  378.  
  379.       image = soft_edge_grad_image(colour)
  380.       vim = ax.imshow(
  381.         image,
  382.         aspect='auto',
  383.         extent=(*r0_span, 0, 100),
  384.         origin='lower',
  385.         interpolation='bicubic',
  386.         alpha=band_alpha,
  387.       )
  388.       vim.set_clip_path(lower_clip_path)
  389.  
  390.       gdivs = 40
  391.       x = np.linspace(*r0_span, gdivs, endpoint=True)
  392.       y = np.array([100 * hit_fn(x_i, ver) for x_i in x])
  393.       xg = np.linspace(0, 1, gdivs, endpoint=True)
  394.       yg = (y - hit_span[0]) / (hit_span[1] - hit_span[0])
  395.       imx = np.interp(xg, yg, xg)
  396.       row = image[0]
  397.       row1 = np.array(resample_lanczos2(row, imx))
  398.       row1 = np.clip(row1, 0, 1)
  399.  
  400.       image1 = np.swapaxes(np.expand_dims(row1, 0), 0, 1)
  401.  
  402.       him = ax.imshow(
  403.         image1,
  404.         aspect='auto',
  405.         extent=(0, r0_max, *hit_span),
  406.         origin='lower',
  407.         interpolation='bicubic',
  408.         alpha=band_alpha,
  409.       )
  410.       him.set_clip_path(upper_clip_path)
  411.  
  412.       lpatch = patches.Patch(
  413.         facecolor=colour,
  414.         edgecolor=None,
  415.         label=label,
  416.         alpha=band_alpha,
  417.       )
  418.       legend_patches.append(lpatch)
  419.       labels.append(label)
  420.  
  421.   legend = ax.legend(
  422.     legend_patches,
  423.     labels,
  424.   )
  425.  
  426.   plt.xlabel("Basic Reproduction Number ($R_0$)")
  427.   plt.ylabel("Herd Immunity Threshold (%)")
  428.   plt.title("Herd Immunity Threshold (of Total Population)\n"
  429.             "for Selected Vaccine Efficacies", fontsize=15)
  430.   major_gl_props = dict(
  431.     color = 'black',
  432.     linestyle = "solid",
  433.     linewidth = 0.2,
  434.     alpha = 0.5
  435.   )
  436.   minor_gl_props = dict(
  437.     color = 'black',
  438.     linestyle = "solid",
  439.     linewidth = 0.2,
  440.     alpha = 0.25
  441.   )
  442.   plt.grid(True, 'major', 'both', **major_gl_props)
  443.   plt.grid(True, 'minor', 'both', **minor_gl_props)
  444.   ax.set_axisbelow(True)
  445.  
  446.   plt.savefig("hit_covid19.png", dpi=300)
  447.   plt.show()
  448.  
  449.  
  450. if __name__ == '__main__':
  451.   main()
  452.  
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!