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 bf9c148c..d043f688 100644 --- a/cara/monte_carlo/data.py +++ b/cara/monte_carlo/data.py @@ -3,9 +3,10 @@ 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 +from cara.monte_carlo.sampleable import LogCustom, LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom sqrt2pi = np.sqrt(2.*np.pi) @@ -101,47 +102,61 @@ 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_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, left=0., right=0.), + max_function=0.2) + + # 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( - 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, ), '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/monte_carlo/sampleable.py b/cara/monte_carlo/sampleable.py index 27907e49..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): @@ -82,6 +84,32 @@ class Custom(SampleableDistribution): return x +class LogCustom(SampleableDistribution): + """ + 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. + 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): + 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 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):