From 6735403b19fab3941286493e90e623eb1eeec38c Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Wed, 11 Oct 2023 16:23:33 +0200 Subject: [PATCH] populate parameters from global store --- caimira/monte_carlo/data.py | 139 +++++++++++++----- caimira/tests/test_monte_carlo_full_models.py | 10 +- 2 files changed, 105 insertions(+), 44 deletions(-) diff --git a/caimira/monte_carlo/data.py b/caimira/monte_carlo/data.py index e0eaa09a..6accd236 100644 --- a/caimira/monte_carlo/data.py +++ b/caimira/monte_carlo/data.py @@ -6,7 +6,8 @@ from scipy.stats import weibull_min import caimira.monte_carlo as mc -from caimira.monte_carlo.sampleable import LogCustom, LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom +from caimira.monte_carlo.sampleable import LogCustom, LogNormal, LogCustomKernel, CustomKernel, Uniform, Custom +import caimira.apps.calculator.global_store.constants as constants sqrt2pi = np.sqrt(2.*np.pi) sqrt2 = np.sqrt(2.) @@ -105,61 +106,121 @@ def integrate(self, dmin, dmax): # 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) +viral_load = np.linspace( + weibull_min.ppf( + constants.covid_overal_vl_data['start'], + c=constants.covid_overal_vl_data['shape_factor'], + scale=constants.covid_overal_vl_data['scale_factor'] + ), + weibull_min.ppf( + constants.covid_overal_vl_data['stop'], + c=constants.covid_overal_vl_data['shape_factor'], + scale=constants.covid_overal_vl_data['scale_factor'] + ), + int(constants.covid_overal_vl_data['num']) +) +frequencies_pdf = weibull_min.pdf( + viral_load, + c=constants.covid_overal_vl_data['shape_factor'], + scale=constants.covid_overal_vl_data['scale_factor'] +) +covid_overal_vl_data = LogCustom(bounds=(constants.covid_overal_vl_data['min_bound'], constants.covid_overal_vl_data['max_bound']), + function=lambda d: np.interp(d, viral_load, frequencies_pdf, constants.covid_overal_vl_data[ + 'interpolation_fp_left'], constants.covid_overal_vl_data['interpolation_fp_right']), + max_function=constants.covid_overal_vl_data['max_function']) + + +def custom_distribution_lookup(dict: dict, key_part: str): + for key, value in dict.items(): + if (key_part in key): + return value['associated_distribution'] + + +def evaluate_reference(reference_variable: str) -> typing.Any: + try: + return eval(reference_variable) + except NameError: + return f"Variable '{reference_variable}' is not defined." + + +def evaluate_custom_distribution(dist: str, params: typing.Dict) -> typing.Any: + if dist == 'Numpy Linear Space (linspace)': + return np.linspace(params['start'], params['stop'], params['num']) + elif dist == 'Numpy Normal Distribution (random.normal)': + return np.random.normal(params['mean_gaussian'], params['standard_deviation_gaussian']) + elif dist == 'Numpy Log-normal Distribution (random.lognormal)': + return np.random.lognormal(params['mean_gaussian'], params['standard_deviation_gaussian']) + elif dist == 'Numpy Uniform Distribution (random.uniform)': + return np.random.uniform(params['low'], params['high']) + else: + raise ValueError('Bad request - variables not found.') + + +def param_evaluation(root: typing.Dict, param: typing.Union[str, typing.Any]) -> typing.Any: + value = root.get(param) + + if isinstance(value, str): + if value.startswith('Ref:'): + reference_variable = value.split(' - ')[1].strip() + return evaluate_reference(reference_variable) + elif value == 'Custom': + distribution: typing.Dict = custom_distribution_lookup( + root, 'custom distribution') + for dist, params in distribution.items(): + return evaluate_custom_distribution(dist, params) + elif isinstance(value, float) or isinstance(value, int): + return value + else: + raise TypeError('Bad request - type not allowed.') # 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) +viable_to_RNA_ratio = Uniform(constants.viable_to_RNA_ratio['low'], constants.viable_to_RNA_ratio['high']) # From discussion with virologists -infectious_dose_distribution = Uniform(10., 100.) +infectious_dose = Uniform(constants.infectious_dose['low'], constants.infectious_dose['high']) # From https://doi.org/10.1101/2021.10.14.21264988 and refererences therein virus_distributions = { 'SARS_CoV_2': mc.SARSCoV2( - 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., - ), + viral_load_in_sputum=param_evaluation(constants.virus_distributions['SARS_CoV_2'], 'viral_load_in_sputum'), + infectious_dose=param_evaluation(constants.virus_distributions['SARS_CoV_2'], 'infectious_dose'), + viable_to_RNA_ratio=param_evaluation(constants.virus_distributions['SARS_CoV_2'], 'viable_to_RNA_ratio'), + transmissibility_factor=param_evaluation(constants.virus_distributions['SARS_CoV_2'], 'transmissibility_factor'), + ), 'SARS_CoV_2_ALPHA': mc.SARSCoV2( - 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, - ), + viral_load_in_sputum=param_evaluation(constants.virus_distributions['SARS_CoV_2_ALPHA'], 'viral_load_in_sputum'), + infectious_dose=param_evaluation(constants.virus_distributions['SARS_CoV_2_ALPHA'], 'infectious_dose'), + viable_to_RNA_ratio=param_evaluation(constants.virus_distributions['SARS_CoV_2_ALPHA'], 'viable_to_RNA_ratio'), + transmissibility_factor=param_evaluation(constants.virus_distributions['SARS_CoV_2_ALPHA'], 'transmissibility_factor'), + ), 'SARS_CoV_2_BETA': mc.SARSCoV2( - 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, - ), + viral_load_in_sputum=param_evaluation(constants.virus_distributions['SARS_CoV_2_BETA'], 'viral_load_in_sputum'), + infectious_dose=param_evaluation(constants.virus_distributions['SARS_CoV_2_BETA'], 'infectious_dose'), + viable_to_RNA_ratio=param_evaluation(constants.virus_distributions['SARS_CoV_2_BETA'], 'viable_to_RNA_ratio'), + transmissibility_factor=param_evaluation(constants.virus_distributions['SARS_CoV_2_BETA'], 'transmissibility_factor'), + ), 'SARS_CoV_2_GAMMA': mc.SARSCoV2( - 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, - ), + viral_load_in_sputum=param_evaluation(constants.virus_distributions['SARS_CoV_2_GAMMA'], 'viral_load_in_sputum'), + infectious_dose=param_evaluation(constants.virus_distributions['SARS_CoV_2_GAMMA'], 'infectious_dose'), + viable_to_RNA_ratio=param_evaluation(constants.virus_distributions['SARS_CoV_2_GAMMA'], 'viable_to_RNA_ratio'), + transmissibility_factor=param_evaluation(constants.virus_distributions['SARS_CoV_2_GAMMA'], 'transmissibility_factor'), + ), 'SARS_CoV_2_DELTA': mc.SARSCoV2( - 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, - ), + viral_load_in_sputum=param_evaluation(constants.virus_distributions['SARS_CoV_2_DELTA'], 'viral_load_in_sputum'), + infectious_dose=param_evaluation(constants.virus_distributions['SARS_CoV_2_DELTA'], 'infectious_dose'), + viable_to_RNA_ratio=param_evaluation(constants.virus_distributions['SARS_CoV_2_DELTA'], 'viable_to_RNA_ratio'), + transmissibility_factor=param_evaluation(constants.virus_distributions['SARS_CoV_2_DELTA'], 'transmissibility_factor'), + ), 'SARS_CoV_2_OMICRON': mc.SARSCoV2( - 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, - ), + viral_load_in_sputum=param_evaluation(constants.virus_distributions['SARS_CoV_2_OMICRON'], 'viral_load_in_sputum'), + infectious_dose=param_evaluation(constants.virus_distributions['SARS_CoV_2_OMICRON'], 'infectious_dose'), + viable_to_RNA_ratio=param_evaluation(constants.virus_distributions['SARS_CoV_2_OMICRON'], 'viable_to_RNA_ratio'), + transmissibility_factor=param_evaluation(constants.virus_distributions['SARS_CoV_2_OMICRON'], 'transmissibility_factor'), + ), } diff --git a/caimira/tests/test_monte_carlo_full_models.py b/caimira/tests/test_monte_carlo_full_models.py index b6e3948c..4fb4c1fd 100644 --- a/caimira/tests/test_monte_carlo_full_models.py +++ b/caimira/tests/test_monte_carlo_full_models.py @@ -5,7 +5,7 @@ import caimira.monte_carlo as mc from caimira import models,data -from caimira.monte_carlo.data import activity_distributions, virus_distributions, expiration_distributions, infectious_dose_distribution, viable_to_RNA_ratio_distribution +from caimira.monte_carlo.data import activity_distributions, virus_distributions, expiration_distributions, infectious_dose, viable_to_RNA_ratio from caimira.apps.calculator.model_generator import build_expiration SAMPLE_SIZE = 500_000 @@ -173,8 +173,8 @@ def skagit_chorale_mc(): presence=models.SpecificInterval(((0, 2.5), )), virus=mc.SARSCoV2( viral_load_in_sputum=10**9, - infectious_dose=infectious_dose_distribution, - viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, + infectious_dose=infectious_dose, + viable_to_RNA_ratio=viable_to_RNA_ratio, transmissibility_factor=1., ), mask=models.Mask.types['No mask'], @@ -214,8 +214,8 @@ def bus_ride_mc(): presence=models.SpecificInterval(((0, 1.67), )), virus=mc.SARSCoV2( viral_load_in_sputum=5*10**8, - infectious_dose=infectious_dose_distribution, - viable_to_RNA_ratio=viable_to_RNA_ratio_distribution, + infectious_dose=infectious_dose, + viable_to_RNA_ratio=viable_to_RNA_ratio, transmissibility_factor=1., ), mask=models.Mask.types['No mask'],