Skip to content

Commit

Permalink
populate parameters from global store
Browse files Browse the repository at this point in the history
  • Loading branch information
@lrdossan
lrdossan committed Oct 11, 2023
1 parent c5b52f0 commit 6735403
Show file tree
Hide file tree
Showing 2 changed files with 105 additions and 44 deletions.
139 changes: 100 additions & 39 deletions caimira/monte_carlo/data.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.)
Expand Down Expand Up @@ -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'),
),
}


Expand Down
10 changes: 5 additions & 5 deletions caimira/tests/test_monte_carlo_full_models.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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'],
Expand Down Expand Up @@ -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'],
Expand Down

0 comments on commit 6735403

Please sign in to comment.