From 0c59dcb905d1ecd360e41b01c17cdaaae1b6c6e5 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Fri, 13 May 2022 16:50:20 +0200 Subject: [PATCH 1/5] Added custom distribution for weibull_min dist --- cara/monte_carlo/data.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index bf9c148c..e70a65de 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -3,6 +3,7 @@ import typing import numpy as np from scipy import special as sp +from scipy.stats import weibull_min import cara.monte_carlo as mc from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom @@ -101,13 +102,25 @@ symptomatic_vl_frequencies = LogCustomKernel( kernel_bandwidth=0.1 ) + +# Weibull distribution with a shape factor of 3.47 and a scale factor of 7.01. +# From https://elifesciences.org/articles/65774 and first line of the figure in +# https://iiif.elifesciences.org/lax:65774%2Felife-65774-fig4-figsupp3-v2.tif/full/1500,/0/default.jpg +viral_load = np.linspace(weibull_min.ppf(0.01, c=3.47, scale=7.01), + weibull_min.ppf(0.99, c=3.47, scale=7.01), 30) +frequencies = weibull_min.pdf(viral_load, c=3.47, scale=7.01) +covid_overal_vl_data = Custom(bounds=(2, 10), function=lambda d: np.interp(d, viral_load, frequencies, right=0., left=0.), max_function=0.16) + + # Derived from data in doi.org/10.1016/j.ijid.2020.09.025 and # https://iosh.com/media/8432/aerosol-infection-risk-hospital-patient-care-full-report.pdf (page 60) viable_to_RNA_ratio_distribution = Uniform(0.01, 0.6) + # From discussion with virologists infectious_dose_distribution = Uniform(10., 100.) + # From https://doi.org/10.1101/2021.10.14.21264988 and refererences therein virus_distributions = { 'SARS_CoV_2': mc.SARSCoV2( From 9f5311f42ec2866ba9c1bb416c0effb43fec6196 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 30 May 2022 10:38:31 +0200 Subject: [PATCH 2/5] Added LogCustom class and changed max function on the distribution --- cara/apps/calculator/__init__.py | 2 +- cara/monte_carlo/data.py | 8 +++++--- cara/monte_carlo/sampleable.py | 24 ++++++++++++++++++++++++ 3 files changed, 30 insertions(+), 4 deletions(-) diff --git a/cara/apps/calculator/__init__.py b/cara/apps/calculator/__init__.py index 8d4e5e45..d7b15698 100644 --- a/cara/apps/calculator/__init__.py +++ b/cara/apps/calculator/__init__.py @@ -33,7 +33,7 @@ from .user import AuthenticatedUser, AnonymousUser # calculator version. If the calculator needs to make breaking changes (e.g. change # form attributes) then it can also increase its MAJOR version without needing to # increase the overall CARA version (found at ``cara.__version__``). -__version__ = "4.1.2" +__version__ = "4.2" class BaseRequestHandler(RequestHandler): diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index e70a65de..886037cf 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -6,7 +6,7 @@ from scipy import special as sp from scipy.stats import weibull_min import cara.monte_carlo as mc -from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom +from cara.monte_carlo.sampleable import LogCustom, LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom sqrt2pi = np.sqrt(2.*np.pi) @@ -108,8 +108,10 @@ symptomatic_vl_frequencies = LogCustomKernel( # https://iiif.elifesciences.org/lax:65774%2Felife-65774-fig4-figsupp3-v2.tif/full/1500,/0/default.jpg viral_load = np.linspace(weibull_min.ppf(0.01, c=3.47, scale=7.01), weibull_min.ppf(0.99, c=3.47, scale=7.01), 30) -frequencies = weibull_min.pdf(viral_load, c=3.47, scale=7.01) -covid_overal_vl_data = Custom(bounds=(2, 10), function=lambda d: np.interp(d, viral_load, frequencies, right=0., left=0.), max_function=0.16) +frequencies_pdf = weibull_min.pdf(viral_load, c=3.47, scale=7.01) +covid_overal_vl_data = LogCustom(bounds=(2, 10), + function=lambda d: np.interp(d, viral_load, frequencies_pdf, right=0., left=0.), + max_function=0.2) # Derived from data in doi.org/10.1016/j.ijid.2020.09.025 and diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 27907e49..fab478e0 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -82,6 +82,30 @@ class Custom(SampleableDistribution): return x +class LogCustom(SampleableDistribution): + """ + Defines a distribution which follows a custom curve vs. the the log (in base 10) + of the random variable. Uses a simple algorithm. This is appropriate for a smooth + distribution function (one should know its maximum). + """ + def __init__(self, bounds: typing.Tuple[float, float], + function: typing.Callable, max_function: float): + self.bounds = bounds + self.function = function + self.max_function = max_function + + def generate_samples(self, size: int) -> float_array_size_n: + fvalue = np.random.uniform(0,self.max_function,size) + x = np.random.uniform(*self.bounds,size) + invalid = np.where(fvalue>self.function(x))[0] + while len(invalid)>0: + fvalue[invalid] = np.random.uniform(0,self.max_function,len(invalid)) + x[invalid] = np.random.uniform(*self.bounds,len(invalid)) + invalid = np.where(fvalue>self.function(x))[0] + + return 10 ** x + + class CustomKernel(SampleableDistribution): """ Defines a distribution which follows a custom curve vs. the From 431fa3b33d965cbf86f3b33d819427700af4bedd Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Mon, 30 May 2022 11:16:30 +0200 Subject: [PATCH 3/5] Updated Custom and LogCustom classes' docstring --- cara/monte_carlo/sampleable.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/cara/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index fab478e0..fe6e72d3 100644 --- a/cara/monte_carlo/sampleable.py +++ b/cara/monte_carlo/sampleable.py @@ -62,7 +62,9 @@ class Custom(SampleableDistribution): """ Defines a distribution which follows a custom curve vs. the random variable. Uses a simple algorithm. This is appropriate for a smooth - distribution function (one should know its maximum). + distribution function. + Note: in max_function, a value slightly above the maximum of the distribution + function should be provided. """ def __init__(self, bounds: typing.Tuple[float, float], function: typing.Callable, max_function: float): @@ -84,9 +86,11 @@ class Custom(SampleableDistribution): class LogCustom(SampleableDistribution): """ - Defines a distribution which follows a custom curve vs. the the log (in base 10) + Defines a distribution which follows a custom curve vs. the log (in base 10) of the random variable. Uses a simple algorithm. This is appropriate for a smooth - distribution function (one should know its maximum). + distribution function. + Note: in max_function, a value slightly above the maximum of the distribution + function should be provided. """ def __init__(self, bounds: typing.Tuple[float, float], function: typing.Callable, max_function: float): From 4f37b2db8d4501a06b34bf469e0be8f2f346aedc Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 31 May 2022 15:21:09 +0200 Subject: [PATCH 4/5] Added the new LogCustom distribution into the virus distributions --- cara/monte_carlo/data.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 886037cf..63c83084 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -126,31 +126,31 @@ infectious_dose_distribution = Uniform(10., 100.) # From https://doi.org/10.1101/2021.10.14.21264988 and refererences therein virus_distributions = { 'SARS_CoV_2': mc.SARSCoV2( - viral_load_in_sputum=symptomatic_vl_frequencies, + viral_load_in_sputum=covid_overal_vl_data, infectious_dose=infectious_dose_distribution, viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, transmissibility_factor=1., ), 'SARS_CoV_2_ALPHA': mc.SARSCoV2( - viral_load_in_sputum=symptomatic_vl_frequencies, + viral_load_in_sputum=covid_overal_vl_data, infectious_dose=infectious_dose_distribution, viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, transmissibility_factor=0.78, ), 'SARS_CoV_2_BETA': mc.SARSCoV2( - viral_load_in_sputum=symptomatic_vl_frequencies, + viral_load_in_sputum=covid_overal_vl_data, infectious_dose=infectious_dose_distribution, viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, transmissibility_factor=0.8, ), 'SARS_CoV_2_GAMMA': mc.SARSCoV2( - viral_load_in_sputum=symptomatic_vl_frequencies, + viral_load_in_sputum=covid_overal_vl_data, infectious_dose=infectious_dose_distribution, viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, transmissibility_factor=0.72, ), 'SARS_CoV_2_DELTA': mc.SARSCoV2( - viral_load_in_sputum=symptomatic_vl_frequencies, + viral_load_in_sputum=covid_overal_vl_data, infectious_dose=infectious_dose_distribution, viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, transmissibility_factor=0.51, From 59e24e0f47f3ea38d46941e1afe26d01fc86f301 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 31 May 2022 17:44:12 +0200 Subject: [PATCH 5/5] Added the new distribution for the viral_load_in_sputum. Tests were updated --- cara/monte_carlo/data.py | 4 ++-- cara/tests/test_monte_carlo_full_models.py | 22 ++++++++++----------- cara/tests/test_predefined_distributions.py | 6 +++--- 3 files changed, 16 insertions(+), 16 deletions(-) diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py index 63c83084..d043f688 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -110,7 +110,7 @@ viral_load = np.linspace(weibull_min.ppf(0.01, c=3.47, scale=7.01), weibull_min.ppf(0.99, c=3.47, scale=7.01), 30) frequencies_pdf = weibull_min.pdf(viral_load, c=3.47, scale=7.01) covid_overal_vl_data = LogCustom(bounds=(2, 10), - function=lambda d: np.interp(d, viral_load, frequencies_pdf, right=0., left=0.), + function=lambda d: np.interp(d, viral_load, frequencies_pdf, left=0., right=0.), max_function=0.2) @@ -156,7 +156,7 @@ virus_distributions = { transmissibility_factor=0.51, ), 'SARS_CoV_2_OMICRON': mc.SARSCoV2( - viral_load_in_sputum=symptomatic_vl_frequencies, + viral_load_in_sputum=covid_overal_vl_data, infectious_dose=infectious_dose_distribution, viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, transmissibility_factor=0.2, diff --git a/cara/tests/test_monte_carlo_full_models.py b/cara/tests/test_monte_carlo_full_models.py index 2272948d..93ebd459 100644 --- a/cara/tests/test_monte_carlo_full_models.py +++ b/cara/tests/test_monte_carlo_full_models.py @@ -308,17 +308,17 @@ def waiting_room_mc(): ) -@retry() +@retry(tries=10) @pytest.mark.parametrize( "mc_model, expected_pi, expected_new_cases, expected_dose, expected_ER", [ - ["shared_office_mc", 5.55, 0.17, 2.699, 809], - ["classroom_mc", 9.58, 1.82, 9.034, 5624], - ["ski_cabin_mc", 16.0, 0.47, 17.315, 7966], + ["shared_office_mc", 5.38, 0.16, 3.350, 1056], + ["classroom_mc", 8.21, 1.56, 11.356, 7416], + ["ski_cabin_mc", 12.92, 0.39, 21.796, 10231], ["skagit_chorale_mc",61.01, 36.53, 84.730, 190422], - ["bus_ride_mc", 10.59, 7.06, 6.65, 5419], - ["gym_mc", 0.43, 0.12, 0.197, 1145], - ["waiting_room_mc", 1.34, 0.18, 0.670, 737], + ["bus_ride_mc", 10.59, 7.06, 6.650, 5419], + ["gym_mc", 0.52, 0.14, 0.249, 1450], + ["waiting_room_mc", 1.53, 0.21, 0.844, 929], ] ) def test_report_models(mc_model, expected_pi, expected_new_cases, @@ -339,10 +339,10 @@ def test_report_models(mc_model, expected_pi, expected_new_cases, @pytest.mark.parametrize( "mask_type, month, expected_pi, expected_dose, expected_ER", [ - ["No mask", "Jul", 8.46, 8.113, 809], - ["Type I", "Jul", 1.44, 0.727, 149], - ["FFP2", "Jul", 0.43, 0.197, 149], - ["Type I", "Feb", 0.54, 0.253, 149], + ["No mask", "Jul", 7.689, 10.050, 1034.435], + ["Type I", "Jul", 1.663, 0.938, 193.52], + ["FFP2", "Jul", 0.523, 0.253, 193.52], + ["Type I", "Feb", 0.659, 0.325, 193.52], ], ) def test_small_shared_office_Geneva(mask_type, month, expected_pi, diff --git a/cara/tests/test_predefined_distributions.py b/cara/tests/test_predefined_distributions.py index b576b7e3..db8c8c4b 100644 --- a/cara/tests/test_predefined_distributions.py +++ b/cara/tests/test_predefined_distributions.py @@ -34,11 +34,11 @@ def test_activity_distributions(distribution, mean, std): # - with a refined precision on the values @pytest.mark.parametrize( "distribution, mean, std",[ - ['SARS_CoV_2', 6.59, 1.74], + ['SARS_CoV_2', 6.22, 1.80], - ['SARS_CoV_2_ALPHA', 6.59, 1.74], + ['SARS_CoV_2_ALPHA', 6.22, 1.80], - ['SARS_CoV_2_GAMMA', 6.59, 1.74], + ['SARS_CoV_2_GAMMA', 6.22, 1.80], ] ) def test_viral_load_logdistribution(distribution, mean, std):