From e23f07944dbb7646d021b09658a3a8e31633f221 Mon Sep 17 00:00:00 2001 From: Edward Wenger Date: Tue, 16 Dec 2025 16:08:27 -0800 Subject: [PATCH 1/2] Refactor to instance-based configuration for thread safety MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Replace static class configuration with instance-based MalariaConfig to enable thread-safe parallel simulations. Each IntrahostComponent now owns a shared_ptr to its config, including its own RNG and SUID generator. API changes: - Add create_config(params) to create config instances from defaults - Add config.update(params) for incremental parameter updates - IntrahostComponent.create() now requires a config argument - Remove static set_params()/update_params() methods - Add config.yaml property for YAML representation New files: - MalariaConfig.h/cpp: Consolidated configuration class - test_thread_safety.py: Tests for parallel simulation isolation - MIGRATION.md: Migration guide from v0.0.4 to v0.1.0 🤖 Generated with [Claude Code](https://claude.com/claude-code) Co-Authored-By: Claude Opus 4.5 --- CMakeLists.txt | 1 + MIGRATION.md | 149 +++++++++++++++ docs/tutorials/tut_challenge.ipynb | 33 +--- docs/tutorials/tut_intro.ipynb | 173 ++---------------- docs/tutorials/tut_multihost.ipynb | 38 +--- docs/tutorials/tut_optuna.ipynb | 60 +----- examples/naive_infection.py | 10 +- include/emodlib/malaria/InfectionMalaria.cpp | 99 +++------- include/emodlib/malaria/InfectionMalaria.h | 35 +--- .../emodlib/malaria/IntrahostComponent.cpp | 60 ++---- include/emodlib/malaria/IntrahostComponent.h | 31 +--- include/emodlib/malaria/MalariaAntibody.cpp | 47 ++--- include/emodlib/malaria/MalariaAntibody.h | 18 +- include/emodlib/malaria/MalariaConfig.cpp | 86 +++++++++ include/emodlib/malaria/MalariaConfig.h | 112 ++++++++++++ .../emodlib/malaria/SusceptibilityMalaria.cpp | 89 ++++----- .../emodlib/malaria/SusceptibilityMalaria.h | 43 +---- src/emodlib/malaria/__init__.py | 77 +++++++- src/malaria.cpp | 79 ++++++-- tests/test_host.py | 46 +++-- tests/test_infection.py | 20 +- tests/test_susceptibility.py | 14 +- tests/test_thread_safety.py | 106 +++++++++++ 23 files changed, 793 insertions(+), 633 deletions(-) create mode 100644 MIGRATION.md create mode 100644 include/emodlib/malaria/MalariaConfig.cpp create mode 100644 include/emodlib/malaria/MalariaConfig.h create mode 100644 tests/test_thread_safety.py diff --git a/CMakeLists.txt b/CMakeLists.txt index c87bc37..988153c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,6 +6,7 @@ find_package(pybind11 CONFIG REQUIRED) # emodlib src files set(EMODLIB_OBJECTS + ${CMAKE_SOURCE_DIR}/include/emodlib/malaria/MalariaConfig.cpp ${CMAKE_SOURCE_DIR}/include/emodlib/malaria/InfectionMalaria.cpp ${CMAKE_SOURCE_DIR}/include/emodlib/malaria/SusceptibilityMalaria.cpp ${CMAKE_SOURCE_DIR}/include/emodlib/malaria/IntrahostComponent.cpp diff --git a/MIGRATION.md b/MIGRATION.md new file mode 100644 index 0000000..aa60160 --- /dev/null +++ b/MIGRATION.md @@ -0,0 +1,149 @@ +# Migration Guide: v0.0.4 to v0.1.0 + +## Breaking Changes + +Version 0.1.0 introduces a thread-safe, instance-based configuration system. The previous static configuration pattern has been removed. + +### Configuration API Changes + +**v0.0.4 (Old)** +```python +from emodlib.malaria import IntrahostComponent + +# Set params globally (affected all instances) +IntrahostComponent.set_params({'Run_Number': 42}) +IntrahostComponent.update_params({'infection_params': {'Antigen_Switch_Rate': 1e-9}}) + +# Create instance (used global params) +ic = IntrahostComponent.create() +``` + +**v0.1.0 (New)** +```python +from emodlib.malaria import IntrahostComponent, create_config + +# Create config instance (thread-safe, independent) +config = create_config({'Run_Number': 42}) + +# Or with nested params +config = create_config({ + 'Run_Number': 42, + 'infection_params': {'Antigen_Switch_Rate': 1e-9} +}) + +# Create instance with config +ic = IntrahostComponent.create(config) +``` + +### Removed Functions + +| Removed | Replacement | +|---------|-------------| +| `IntrahostComponent.set_params(params)` | `config = create_config(params)` | +| `IntrahostComponent.update_params(params)` | `config.update(params)` | +| `IntrahostComponent.params` | `config.yaml` or config properties | + +### Incremental Updates with config.update() + +The `update()` method provides the same incremental merge behavior as the old `update_params()`: + +```python +# Create config with initial params +config = create_config({'Run_Number': 1, 'Max_Individual_Infections': 3}) + +# Later, update just one param (keeps Run_Number=1, Max_Individual_Infections=3) +config.update({'infection_params': {'Antigen_Switch_Rate': 1e-8}}) + +# Method chaining also works +config.update({'Run_Number': 2}).update({'infection_params': {'Antigen_Switch_Rate': 1e-9}}) +``` + +### Creating Susceptibility and Infection Objects + +**v0.0.4** +```python +from emodlib.malaria import Susceptibility, Infection + +s = Susceptibility.create() +inf = Infection.create(susceptibility=s, hepatocytes=1) +``` + +**v0.1.0** +```python +from emodlib.malaria import Susceptibility, Infection, create_config + +config = create_config() +s = Susceptibility.create(config) +inf = Infection.create(susceptibility=s, config=config, hepatocytes=1) +``` + +### Viewing Configuration + +**v0.0.4** +```python +print(IntrahostComponent.params.yaml) +``` + +**v0.1.0** +```python +config = create_config() +print(config.yaml) +``` + +### Multi-Simulation Workflows + +**v0.0.4** - Simulations with different parameters required careful sequencing: +```python +IntrahostComponent.set_params({'Run_Number': 1}) +ic1 = IntrahostComponent.create() + +IntrahostComponent.set_params({'Run_Number': 2}) # Changed global state! +ic2 = IntrahostComponent.create() +``` + +**v0.1.0** - Each config is independent (thread-safe): +```python +config1 = create_config({'Run_Number': 1}) +config2 = create_config({'Run_Number': 2}) + +ic1 = IntrahostComponent.create(config1) +ic2 = IntrahostComponent.create(config2) + +# Can run in parallel without interference +``` + +## Quick Migration Checklist + +1. Add `create_config` to imports: + ```python + from emodlib.malaria import IntrahostComponent, create_config + ``` + +2. Replace `set_params()` with `create_config()`: + ```python + # Before + IntrahostComponent.set_params({'Run_Number': 42}) + + # After + config = create_config({'Run_Number': 42}) + ``` + +3. Replace `update_params()` with `config.update()`: + ```python + # Before + IntrahostComponent.update_params({'infection_params': {'X': 1}}) + + # After + config.update({'infection_params': {'X': 1}}) + ``` + +4. Pass config to `create()`: + ```python + # Before + ic = IntrahostComponent.create() + + # After + ic = IntrahostComponent.create(config) + ``` + +5. For `Susceptibility.create()` and `Infection.create()`, add config parameter. diff --git a/docs/tutorials/tut_challenge.ipynb b/docs/tutorials/tut_challenge.ipynb index 9d93611..e528a21 100644 --- a/docs/tutorials/tut_challenge.ipynb +++ b/docs/tutorials/tut_challenge.ipynb @@ -10,40 +10,17 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "\n", - "from emodlib.malaria import IntrahostComponent" - ] + "source": "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom emodlib.malaria import IntrahostComponent, create_config" }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "duration = 300\n", - "\n", - "asexuals = np.zeros(duration)\n", - "gametocytes = np.zeros(duration)\n", - "\n", - "ic = IntrahostComponent.create()\n", - "ic.challenge()\n", - "\n", - "for t in range(duration):\n", - " ic.update(dt=1)\n", - " asexuals[t] = ic.parasite_density\n", - " gametocytes[t] = ic.gametocyte_density\n", - "\n", - "df = pd.DataFrame({'days': range(duration),\n", - " 'parasite_density': asexuals,\n", - " 'gametocyte_density': gametocytes}).set_index('days')" - ] + "source": "duration = 300\n\nasexuals = np.zeros(duration)\ngametocytes = np.zeros(duration)\n\nconfig = create_config()\nic = IntrahostComponent.create(config)\nic.challenge()\n\nfor t in range(duration):\n ic.update(dt=1)\n asexuals[t] = ic.parasite_density\n gametocytes[t] = ic.gametocyte_density\n\ndf = pd.DataFrame({'days': range(duration),\n 'parasite_density': asexuals,\n 'gametocyte_density': gametocytes}).set_index('days')" }, { "cell_type": "code", @@ -205,4 +182,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/docs/tutorials/tut_intro.ipynb b/docs/tutorials/tut_intro.ipynb index 4e12df3..1e9e3cd 100644 --- a/docs/tutorials/tut_intro.ipynb +++ b/docs/tutorials/tut_intro.ipynb @@ -34,175 +34,40 @@ "attachments": {}, "cell_type": "markdown", "metadata": {}, - "source": [ - "First let's look at the default malaria intra-host model parameters that are loaded on import..." - ] + "source": "First let's create a configuration instance with default malaria intra-host model parameters and look at them..." }, { "cell_type": "code", - "execution_count": 9, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Base_Gametocyte_Mosquito_Survival_Rate: 0.002011099\n", - "Cytokine_Gametocyte_Inactivation: 0.01667\n", - "Falciparum_MSP_Variants: 32\n", - "Falciparum_Nonspecific_Types: 76\n", - "Falciparum_PfEMP1_Variants: 1070\n", - "Run_Number: 54321\n", - "Max_Individual_Infections: 5\n", - "infection_params:\n", - " Antibody_IRBC_Kill_Rate: 1.596\n", - " Antigen_Switch_Rate: 7.645570124964182e-10\n", - " Base_Gametocyte_Fraction_Male: 0.2\n", - " Base_Gametocyte_Production_Rate: 0.06150582\n", - " Base_Incubation_Period: 7\n", - " Gametocyte_Stage_Survival_Rate: 0.588569307\n", - " MSP1_Merozoite_Kill_Fraction: 0.511735322\n", - " Merozoites_Per_Hepatocyte: 15000\n", - " Merozoites_Per_Schizont: 16\n", - " Nonspecific_Antigenicity_Factor: 0.415111634\n", - " Number_Of_Asexual_Cycles_Without_Gametocytes: 1\n", - " RBC_Destruction_Multiplier: 3.29\n", - "susceptibility_params:\n", - " Antibody_CSP_Decay_Days: 90\n", - " Antibody_Capacity_Growth_Rate: 0.09\n", - " Antibody_Memory_Level: 0.34\n", - " Antibody_Stimulation_C50: 30\n", - " Erythropoiesis_Anemia_Effect: 3.5\n", - " Fever_IRBC_Kill_Rate: 1.4\n", - " Maternal_Antibody_Decay_Rate: 0.01\n", - " Max_MSP1_Antibody_Growthrate: 0.045\n", - " Min_Adapted_Response: 0.05\n", - " Nonspecific_Antibody_Growth_Rate_Factor: 0.5\n", - " Pyrogenic_Threshold: 15000.0\n", - "\n" - ] - } - ], - "source": [ - "from emodlib.malaria import IntrahostComponent\n", - "print(IntrahostComponent.params.yaml)" - ] + "outputs": [], + "source": "from emodlib.malaria import IntrahostComponent, create_config\n\nconfig = create_config()\nprint(config.yaml)" }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, - "source": [ - "We can layer in modifications to arbitrary parameters (on top of existing defaults) by passing a dictionary to `update_params` function..." - ] + "source": "We can create a new configuration with modifications to arbitrary parameters (on top of existing defaults) by passing a dictionary to `create_config` function..." }, { "cell_type": "code", - "execution_count": 10, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Base_Gametocyte_Mosquito_Survival_Rate: 0.002011099\n", - "Cytokine_Gametocyte_Inactivation: 0.01667\n", - "Falciparum_MSP_Variants: 32\n", - "Falciparum_Nonspecific_Types: 76\n", - "Falciparum_PfEMP1_Variants: 1070\n", - "Run_Number: 6789\n", - "Max_Individual_Infections: 5\n", - "infection_params:\n", - " Antibody_IRBC_Kill_Rate: 1.596\n", - " Antigen_Switch_Rate: 7.645570124964182e-10\n", - " Base_Gametocyte_Fraction_Male: 0.2\n", - " Base_Gametocyte_Production_Rate: 0.06150582\n", - " Base_Incubation_Period: 7\n", - " Gametocyte_Stage_Survival_Rate: 0.588569307\n", - " MSP1_Merozoite_Kill_Fraction: 0.511735322\n", - " Merozoites_Per_Hepatocyte: 15000\n", - " Merozoites_Per_Schizont: 16\n", - " Nonspecific_Antigenicity_Factor: 0.415111634\n", - " Number_Of_Asexual_Cycles_Without_Gametocytes: 1\n", - " RBC_Destruction_Multiplier: 3.29\n", - "susceptibility_params:\n", - " Antibody_CSP_Decay_Days: 99\n", - " Antibody_Capacity_Growth_Rate: 0.09\n", - " Antibody_Memory_Level: 0.34\n", - " Antibody_Stimulation_C50: 30\n", - " Erythropoiesis_Anemia_Effect: 3.5\n", - " Fever_IRBC_Kill_Rate: 1.4\n", - " Maternal_Antibody_Decay_Rate: 0.01\n", - " Max_MSP1_Antibody_Growthrate: 0.045\n", - " Min_Adapted_Response: 0.05\n", - " Nonspecific_Antibody_Growth_Rate_Factor: 0.5\n", - " Pyrogenic_Threshold: 15000.0\n", - "\n" - ] - } - ], - "source": [ - "IntrahostComponent.update_params(dict(Run_Number=6789, susceptibility_params=dict(Antibody_CSP_Decay_Days=99)))\n", - "print(IntrahostComponent.params.yaml)" - ] + "outputs": [], + "source": "config2 = create_config({'Run_Number': 6789, 'susceptibility_params': {'Antibody_CSP_Decay_Days': 99}})\nprint(config2.yaml)" }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, - "source": [ - "And we can layer in modified parameters on top of the module defaults, resetting those not explicitly set, with `set_params` function call" - ] + "source": "Each config instance is independent. Note that `config` above still has its original defaults (Run_Number=12345), while `config2` has our modified values. Creating a new config always starts from defaults..." }, { "cell_type": "code", - "execution_count": 11, + "execution_count": null, "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Base_Gametocyte_Mosquito_Survival_Rate: 0.002011099\n", - "Cytokine_Gametocyte_Inactivation: 0.01667\n", - "Falciparum_MSP_Variants: 32\n", - "Falciparum_Nonspecific_Types: 76\n", - "Falciparum_PfEMP1_Variants: 1070\n", - "Run_Number: 54321\n", - "Max_Individual_Infections: 5\n", - "infection_params:\n", - " Antibody_IRBC_Kill_Rate: 1.596\n", - " Antigen_Switch_Rate: 7.645570124964182e-10\n", - " Base_Gametocyte_Fraction_Male: 0.2\n", - " Base_Gametocyte_Production_Rate: 0.06150582\n", - " Base_Incubation_Period: 7\n", - " Gametocyte_Stage_Survival_Rate: 0.588569307\n", - " MSP1_Merozoite_Kill_Fraction: 0.511735322\n", - " Merozoites_Per_Hepatocyte: 15000\n", - " Merozoites_Per_Schizont: 16\n", - " Nonspecific_Antigenicity_Factor: 0.415111634\n", - " Number_Of_Asexual_Cycles_Without_Gametocytes: 1\n", - " RBC_Destruction_Multiplier: 3.29\n", - "susceptibility_params:\n", - " Antibody_CSP_Decay_Days: 90\n", - " Antibody_Capacity_Growth_Rate: 0.09\n", - " Antibody_Memory_Level: 0.34\n", - " Antibody_Stimulation_C50: 30\n", - " Erythropoiesis_Anemia_Effect: 3.5\n", - " Fever_IRBC_Kill_Rate: 1.4\n", - " Maternal_Antibody_Decay_Rate: 0.01\n", - " Max_MSP1_Antibody_Growthrate: 0.045\n", - " Min_Adapted_Response: 0.05\n", - " Nonspecific_Antibody_Growth_Rate_Factor: 0.5\n", - " Pyrogenic_Threshold: 15000.0\n", - "\n" - ] - } - ], - "source": [ - "IntrahostComponent.set_params(dict(Run_Number=54321))\n", - "print(IntrahostComponent.params.yaml)" - ] + "outputs": [], + "source": "config3 = create_config({'Run_Number': 99999})\nprint(f\"config.random_seed = {config.random_seed}\")\nprint(f\"config2.random_seed = {config2.random_seed}\")\nprint(f\"config3.random_seed = {config3.random_seed}\")" }, { "attachments": {}, @@ -216,20 +81,14 @@ "attachments": {}, "cell_type": "markdown", "metadata": {}, - "source": [ - "Now let's import the python bindings for the `emodlib.malaria` library and create a naive host..." - ] + "source": "Now let's create a naive host using one of our configs..." }, { "cell_type": "code", - "execution_count": 12, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "from emodlib.malaria import IntrahostComponent\n", - "\n", - "ic = IntrahostComponent.create()" - ] + "source": "ic = IntrahostComponent.create(config)" }, { "attachments": {}, @@ -375,4 +234,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/docs/tutorials/tut_multihost.ipynb b/docs/tutorials/tut_multihost.ipynb index ab77fd1..dcc19e5 100644 --- a/docs/tutorials/tut_multihost.ipynb +++ b/docs/tutorials/tut_multihost.ipynb @@ -10,45 +10,17 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "import numpy as np\n", - "import pandas as pd\n", - "import matplotlib.pyplot as plt\n", - "import xarray as xr\n", - "\n", - "from emodlib.malaria import IntrahostComponent" - ] + "source": "import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport xarray as xr\n\nfrom emodlib.malaria import IntrahostComponent, create_config" }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "duration = 350\n", - "n_people = 25\n", - "\n", - "asexuals = np.zeros((n_people, duration))\n", - "gametocytes = np.zeros((n_people, duration))\n", - "\n", - "pp = [IntrahostComponent.create() for _ in range(n_people)]\n", - "_ = [p.challenge() for p in pp]\n", - "\n", - "for t in range(duration):\n", - " for i, p in enumerate(pp):\n", - " p.update(dt=1)\n", - " asexuals[i, t] = p.parasite_density\n", - " gametocytes[i, t] = p.gametocyte_density\n", - " \n", - "da = xr.DataArray(dims=('individual', 'time', 'channel'),\n", - " coords=(range(n_people), range(duration), ['parasite_density', 'gametocyte_density']))\n", - " \n", - "da.loc[dict(channel='parasite_density')] = asexuals\n", - "da.loc[dict(channel='gametocyte_density')] = gametocytes" - ] + "source": "duration = 350\nn_people = 25\n\nasexuals = np.zeros((n_people, duration))\ngametocytes = np.zeros((n_people, duration))\n\nconfig = create_config()\npp = [IntrahostComponent.create(config) for _ in range(n_people)]\n_ = [p.challenge() for p in pp]\n\nfor t in range(duration):\n for i, p in enumerate(pp):\n p.update(dt=1)\n asexuals[i, t] = p.parasite_density\n gametocytes[i, t] = p.gametocyte_density\n \nda = xr.DataArray(dims=('individual', 'time', 'channel'),\n coords=(range(n_people), range(duration), ['parasite_density', 'gametocyte_density']))\n \nda.loc[dict(channel='parasite_density')] = asexuals\nda.loc[dict(channel='gametocyte_density')] = gametocytes" }, { "cell_type": "code", @@ -657,4 +629,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/docs/tutorials/tut_optuna.ipynb b/docs/tutorials/tut_optuna.ipynb index 1e0e40e..e4b02ad 100644 --- a/docs/tutorials/tut_optuna.ipynb +++ b/docs/tutorials/tut_optuna.ipynb @@ -26,39 +26,10 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "import numpy as np\n", - "import pandas as pd\n", - "import xarray as xr\n", - "import matplotlib.pyplot as plt\n", - "\n", - "from emodlib.malaria import IntrahostComponent\n", - "\n", - "\n", - "def multiple_challenges(n_people, duration):\n", - " \n", - " asexuals = np.zeros((n_people, duration))\n", - " gametocytes = np.zeros((n_people, duration))\n", - " pp = [IntrahostComponent.create() for _ in range(n_people)]\n", - " _ = [p.challenge() for p in pp]\n", - "\n", - " for t in range(duration):\n", - " for i, p in enumerate(pp):\n", - " p.update(dt=1)\n", - " asexuals[i, t] = p.parasite_density\n", - " gametocytes[i, t] = p.gametocyte_density\n", - " \n", - " da = xr.DataArray(dims=('individual', 'time', 'channel'),\n", - " coords=(range(n_people), range(duration), ['parasite_density', 'gametocyte_density']))\n", - " \n", - " da.loc[dict(channel='parasite_density')] = asexuals\n", - " da.loc[dict(channel='gametocyte_density')] = gametocytes\n", - " \n", - " return da" - ] + "source": "import numpy as np\nimport pandas as pd\nimport xarray as xr\nimport matplotlib.pyplot as plt\n\nfrom emodlib.malaria import IntrahostComponent, create_config\n\n\ndef multiple_challenges(config, n_people, duration):\n \n asexuals = np.zeros((n_people, duration))\n gametocytes = np.zeros((n_people, duration))\n pp = [IntrahostComponent.create(config) for _ in range(n_people)]\n _ = [p.challenge() for p in pp]\n\n for t in range(duration):\n for i, p in enumerate(pp):\n p.update(dt=1)\n asexuals[i, t] = p.parasite_density\n gametocytes[i, t] = p.gametocyte_density\n \n da = xr.DataArray(dims=('individual', 'time', 'channel'),\n coords=(range(n_people), range(duration), ['parasite_density', 'gametocyte_density']))\n \n da.loc[dict(channel='parasite_density')] = asexuals\n da.loc[dict(channel='gametocyte_density')] = gametocytes\n \n return da" }, { "attachments": {}, @@ -83,33 +54,14 @@ "attachments": {}, "cell_type": "markdown", "metadata": {}, - "source": [ - "Then we'll define our objective function:\n", - "- log-uniform sampling of the antigen switching rate within a defined range\n", - "- setting the model parameters\n", - "- running the multi-individual challenge time-series\n", - "- returning the mean infection-duration value" - ] + "source": "Then we'll define our objective function:\n- log-uniform sampling of the antigen switching rate within a defined range\n- creating a model config with that parameter value\n- running the multi-individual challenge time-series\n- returning the mean infection-duration value" }, { "cell_type": "code", - "execution_count": 3, + "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "def objective(trial):\n", - " \n", - " n_people = 50\n", - " duration = 500\n", - "\n", - " antigen_switch_rate = trial.suggest_float(\"Antigen_Switch_Rate\", 5e-10, 5e-8, log=True)\n", - " IntrahostComponent.set_params(dict(infection_params=dict(Antigen_Switch_Rate=antigen_switch_rate)))\n", - " \n", - " da = multiple_challenges(duration=duration, n_people=n_people)\n", - " infection_durations = get_last_nonzero_by_row(da.sel(channel='parasite_density').values)[1]\n", - " \n", - " return infection_durations.mean()" - ] + "source": "def objective(trial):\n \n n_people = 50\n duration = 500\n\n antigen_switch_rate = trial.suggest_float(\"Antigen_Switch_Rate\", 5e-10, 5e-8, log=True)\n config = create_config({'infection_params': {'Antigen_Switch_Rate': antigen_switch_rate}})\n \n da = multiple_challenges(config, duration=duration, n_people=n_people)\n infection_durations = get_last_nonzero_by_row(da.sel(channel='parasite_density').values)[1]\n \n return infection_durations.mean()" }, { "attachments": {}, @@ -2318,4 +2270,4 @@ }, "nbformat": 4, "nbformat_minor": 2 -} +} \ No newline at end of file diff --git a/examples/naive_infection.py b/examples/naive_infection.py index 9ca7871..520a9fb 100644 --- a/examples/naive_infection.py +++ b/examples/naive_infection.py @@ -2,16 +2,16 @@ import numpy as np import pandas as pd -from emodlib.malaria import IntrahostComponent +from emodlib.malaria import IntrahostComponent, create_config -def run_challenge(duration): +def run_challenge(config, duration): asexuals = np.zeros(duration) gametocytes = np.zeros(duration) fevers = np.zeros(duration) infects = np.zeros(duration) - ic = IntrahostComponent.create() + ic = IntrahostComponent.create(config) ic.challenge() for t in range(duration): @@ -33,9 +33,9 @@ def run_challenge(duration): if __name__ == "__main__": - IntrahostComponent.set_params() # default params + config = create_config() # default params - df = run_challenge(duration=300) + df = run_challenge(config, duration=300) print(df.head(10)) fig, ax = plt.subplots(1, 1, figsize=(8, 3)) diff --git a/include/emodlib/malaria/InfectionMalaria.cpp b/include/emodlib/malaria/InfectionMalaria.cpp index 6b917cb..87e963b 100644 --- a/include/emodlib/malaria/InfectionMalaria.cpp +++ b/include/emodlib/malaria/InfectionMalaria.cpp @@ -12,7 +12,7 @@ #include "emodlib/utils/Common.h" #include "emodlib/utils/Sigmoid.h" -#include "IntrahostComponent.h" +#include "MalariaConfig.h" #include "SusceptibilityMalaria.h" @@ -22,47 +22,7 @@ namespace emodlib namespace malaria { - // TODO: emodlib#8 (boost + enums) - // ParasiteSwitchType::Enum Infection::params::parasite_switch_type = ParasiteSwitchType::RATE_PER_PARASITE_7VARS; - // MalariaStrains::Enum Infection::params::malaria_strains = MalariaStrains::FALCIPARUM_RANDOM_STRAIN; - - float Infection::params::incubation_period = 7.0f; // liver stage duration - - float Infection::params::antibody_IRBC_killrate = DEFAULT_ANTIBODY_IRBC_KILLRATE; - float Infection::params::non_specific_antigenicity = DEFAULT_NON_SPECIFIC_ANTIGENICITY; - float Infection::params::MSP1_merozoite_kill = DEFAULT_MSP1_MEROZOITE_KILL; - float Infection::params::gametocyte_stage_survival = DEFAULT_GAMETOCYTE_STAGE_SURVIVAL; - float Infection::params::base_gametocyte_sexratio = DEFAULT_BASE_GAMETOCYTE_SEX_RATIO; - float Infection::params::base_gametocyte_production = DEFAULT_BASE_GAMETOCYTE_PRODUCTION; - float Infection::params::antigen_switch_rate = DEFAULT_ANTIGEN_SWITCH_RATE; - float Infection::params::merozoites_per_hepatocyte = DEFAULT_MEROZOITES_PER_HEPATOCYTE; - float Infection::params::merozoites_per_schizont = DEFAULT_MEROZOITES_PER_SCHIZONT; - float Infection::params::RBC_destruction_multiplier = DEFAULT_RBC_DESTRUCTION_MULTIPLIER; - int Infection::params::n_asexual_cycles_wo_gametocytes = DEFAULT_ASEXUAL_CYCLES_WITHOUT_GAMETOCYTES; - - - suids::distributed_generator Infection::infectionSuidGenerator(0, 0); - - - void Infection::params::Configure(const ParamSet& pset) - { - incubation_period = pset["Base_Incubation_Period"].cast(); // TODO: emodlib#6 (gaussian distribution) - - antibody_IRBC_killrate = pset["Antibody_IRBC_Kill_Rate"].cast(); - non_specific_antigenicity = pset["Nonspecific_Antigenicity_Factor"].cast(); - MSP1_merozoite_kill = pset["MSP1_Merozoite_Kill_Fraction"].cast(); - gametocyte_stage_survival = pset["Gametocyte_Stage_Survival_Rate"].cast(); - base_gametocyte_sexratio = pset["Base_Gametocyte_Fraction_Male"].cast(); - base_gametocyte_production = pset["Base_Gametocyte_Production_Rate"].cast(); - antigen_switch_rate = pset["Antigen_Switch_Rate"].cast(); - merozoites_per_hepatocyte = pset["Merozoites_Per_Hepatocyte"].cast(); - merozoites_per_schizont = pset["Merozoites_Per_Schizont"].cast(); - RBC_destruction_multiplier = pset["RBC_Destruction_Multiplier"].cast(); - n_asexual_cycles_wo_gametocytes = pset["Number_Of_Asexual_Cycles_Without_Gametocytes"].cast(); - } - - - Infection::Infection() + Infection::Infection(MalariaConfig* config) : suid(suids::nil_suid()) , m_liver_stage_timer(0.0f) @@ -87,13 +47,14 @@ namespace emodlib , m_gametosexratio(0.0) , immunity(nullptr) + , m_config(config) { } - Infection* Infection::Create(Susceptibility* _susceptibility, int initial_hepatocytes) + Infection* Infection::Create(Susceptibility* _susceptibility, MalariaConfig* config, int initial_hepatocytes) { - Infection *newinfection = new Infection(); + Infection *newinfection = new Infection(config); newinfection->Initialize(_susceptibility, initial_hepatocytes); return newinfection; @@ -101,7 +62,7 @@ namespace emodlib void Infection::Initialize(Susceptibility* _susceptibility, int initial_hepatocytes) { - suid = infectionSuidGenerator(); // next suid from generator + suid = m_config->suidGenerator(); // next suid from generator m_hepatocytes = initial_hepatocytes; // Here we set the antigenic repertoire of the infection @@ -110,14 +71,14 @@ namespace emodlib // Recker, M., S. Nee, et al. (2004). "Transient cross-reactive immune responses can orchestrate antigenic variation in malaria." Nature 429(6991): 555-558. // In our model, not all antigens are expressed at the same time, but switching occurs. This just sets the total repertoire - auto rng = IntrahostComponent::p_rng; + auto rng = m_config->rng; - m_MSPtype = rng->uniformZeroToN16(IntrahostComponent::params::falciparumMSPVars); - m_nonspectype = rng->uniformZeroToN16(IntrahostComponent::params::falciparumNonSpecTypes); + m_MSPtype = rng->uniformZeroToN16(m_config->falciparumMSPVars); + m_nonspectype = rng->uniformZeroToN16(m_config->falciparumNonSpecTypes); for (int i = 0; i < CLONAL_PfEMP1_VARIANTS; i++) { - m_IRBCtype[i] = rng->uniformZeroToN16(IntrahostComponent::params::falciparumPfEMP1Vars); + m_IRBCtype[i] = rng->uniformZeroToN16(m_config->falciparumPfEMP1Vars); m_minor_epitope_type[i] = rng->uniformZeroToN16(MINOR_EPITOPE_VARS_PER_SET) + MINOR_EPITOPE_VARS_PER_SET * m_nonspectype; } @@ -206,7 +167,7 @@ namespace emodlib // --- process start of asexual phase if the incubation period is over and there are still hepatocytes // ---------------------------------------------------------------------------------------------------------------------- if (m_asexual_phase == AsexualCycleStatus::NoAsexualCycle && - m_liver_stage_timer >= Infection::params::incubation_period) + m_liver_stage_timer >= m_config->incubation_period) { m_IRBC_count.assign(CLONAL_PfEMP1_VARIANTS, 0); @@ -216,7 +177,7 @@ namespace emodlib #pragma loop(hint_parallel(8)) for ( int i=0; imerozoites_per_hepatocyte / INITIAL_PFEMP1_VARIANTS); immunity->UpdateActiveAntibody( m_PfEMP1_antibodies[i], m_minor_epitope_type[i], m_IRBCtype[i] ); // insert into set of antigens the immune system has ever "seen" } @@ -236,7 +197,7 @@ namespace emodlib double RBCavailability = immunity->get_RBC_availability(); // Merozoite survival limited at very low density according to density-dependent probability-of-success formula - double merozoitesurvival = std::max(0.0, (1.0 - Infection::params::MSP1_merozoite_kill * m_MSP_antibody->GetAntibodyConcentration() ) * EXPCDF(-RBCavailability / MEROZOITE_LIMITING_RBC_THRESHOLD)); + double merozoitesurvival = std::max(0.0, (1.0 - m_config->MSP1_merozoite_kill * m_MSP_antibody->GetAntibodyConcentration() ) * EXPCDF(-RBCavailability / MEROZOITE_LIMITING_RBC_THRESHOLD)); // How many rupture for this infection handed to suscept object for total stimulation calculations int64_t totalIRBC = 0; @@ -264,7 +225,7 @@ namespace emodlib } // Uninfected RBC killing diminishing in proportion to RBC availability - double destruction_factor_ = std::max(1.0, Infection::params::RBC_destruction_multiplier * EXPCDF(-RBCavailability / MEROZOITE_LIMITING_RBC_THRESHOLD) ); + double destruction_factor_ = std::max(1.0, m_config->RBC_destruction_multiplier * EXPCDF(-RBCavailability / MEROZOITE_LIMITING_RBC_THRESHOLD) ); immunity->remove_RBCs( totalIRBC, m_malegametocytes[0] + m_femalegametocytes[0], destruction_factor_ ); // reset timer for next asexual cycle @@ -298,12 +259,12 @@ namespace emodlib if ( m_IRBC_count[j] <= 0 ) continue; // no IRBC means no contribution to next time step int64_t temp_sum_IRBC = 0; - if (Infection::params::antigen_switch_rate > 0) + if (m_config->antigen_switch_rate > 0) { #pragma loop(hint_parallel(8)) for ( int iswitch = 0; iswitch < SWITCHING_IRBC_VARIANT_COUNT; iswitch++ ) { - switchingIRBC[iswitch] = (iswitch < 7) ? IntrahostComponent::p_rng->Poisson(Infection::params::antigen_switch_rate * m_IRBC_count[j]) : 0; + switchingIRBC[iswitch] = (iswitch < 7) ? m_config->rng->Poisson(m_config->antigen_switch_rate * m_IRBC_count[j]) : 0; } // now test to see if these add up to more than 100 percent @@ -321,13 +282,13 @@ namespace emodlib } // Now switch to next stages based on predetermined number of switching IRBC's - tmpIRBCcount[j] = int64_t(tmpIRBCcount[j] + ((1.0 - m_gametorate) * m_IRBC_count[j] - temp_sum_IRBC) * Infection::params::merozoites_per_schizont * merozoitesurvival); - if (Infection::params::antigen_switch_rate > 0) + tmpIRBCcount[j] = int64_t(tmpIRBCcount[j] + ((1.0 - m_gametorate) * m_IRBC_count[j] - temp_sum_IRBC) * m_config->merozoites_per_schizont * merozoitesurvival); + if (m_config->antigen_switch_rate > 0) { #pragma loop(hint_parallel(8)) for ( int iswitch = 0; iswitch < SWITCHING_IRBC_VARIANT_COUNT; iswitch++) { - tmpIRBCcount[(j + iswitch + 1) % CLONAL_PfEMP1_VARIANTS] = int64_t(tmpIRBCcount[(j + iswitch + 1) % CLONAL_PfEMP1_VARIANTS] + switchingIRBC[iswitch] * Infection::params::merozoites_per_schizont * merozoitesurvival); + tmpIRBCcount[(j + iswitch + 1) % CLONAL_PfEMP1_VARIANTS] = int64_t(tmpIRBCcount[(j + iswitch + 1) % CLONAL_PfEMP1_VARIANTS] + switchingIRBC[iswitch] * m_config->merozoites_per_schizont * merozoitesurvival); } } } @@ -339,10 +300,10 @@ namespace emodlib void Infection::malariaCycleGametocytes(double merozoitesurvival) { // set gametocyte production rate for next cycle - if ( m_asexual_cycle_count >= Infection::params::n_asexual_cycles_wo_gametocytes ) + if ( m_asexual_cycle_count >= m_config->n_asexual_cycles_wo_gametocytes ) { - m_gametorate = double(Infection::params::base_gametocyte_production); // gametocyte production used by all switching calculations, here is where factors modifying production would go - m_gametosexratio = double(Infection::params::base_gametocyte_sexratio); + m_gametorate = double(m_config->base_gametocyte_production); // gametocyte production used by all switching calculations, here is where factors modifying production would go + m_gametosexratio = double(m_config->base_gametocyte_sexratio); } // check for valid range of input, and only create next cycle if valid @@ -352,13 +313,13 @@ namespace emodlib //process gametocytes--5 stages--Sinden, R. E., G. A. Butcher, et al. (1996). "Regulation of Infectivity of Plasmodium to the Mosquito Vector." Advances in Parasitology 38: 53-117. for (int j = GametocyteStages::Mature; j > 0; j--) // move developing gametocytes forward a class, moving backwards through stages to not override next stage's values { - m_malegametocytes[j] = int64_t(m_malegametocytes[j] + m_malegametocytes[j - 1] * Infection::params::gametocyte_stage_survival); + m_malegametocytes[j] = int64_t(m_malegametocytes[j] + m_malegametocytes[j - 1] * m_config->gametocyte_stage_survival); m_malegametocytes[j - 1] = 0; if (m_malegametocytes[j] < 1) m_malegametocytes[j] = 0; - m_femalegametocytes[j] = int64_t(m_femalegametocytes[j] + (m_femalegametocytes[j - 1] * Infection::params::gametocyte_stage_survival)); + m_femalegametocytes[j] = int64_t(m_femalegametocytes[j] + (m_femalegametocytes[j - 1] * m_config->gametocyte_stage_survival)); m_femalegametocytes[j - 1] = 0; if (m_femalegametocytes[j] < 1) @@ -371,8 +332,8 @@ namespace emodlib { // review of production rates and sex ratios in Sinden, R. E., G. A. Butcher, et al. (1996). "Regulation of Infectivity of Plasmodium to the Mosquito Vector." Advances in Parasitology 38: 53-117. // each factor may be variable, but here we leave it constant at the moment, conservatively not including the possible senescence of transmission in late infection - m_malegametocytes[GametocyteStages::Stage0] = int64_t(m_malegametocytes[GametocyteStages::Stage0] + m_IRBC_count[j] * m_gametorate * m_gametosexratio * merozoitesurvival * Infection::params::merozoites_per_schizont); - m_femalegametocytes[GametocyteStages::Stage0] = int64_t(m_femalegametocytes[GametocyteStages::Stage0] + m_IRBC_count[j] * m_gametorate * (1.0 - m_gametosexratio) * merozoitesurvival * Infection::params::merozoites_per_schizont); + m_malegametocytes[GametocyteStages::Stage0] = int64_t(m_malegametocytes[GametocyteStages::Stage0] + m_IRBC_count[j] * m_gametorate * m_gametosexratio * merozoitesurvival * m_config->merozoites_per_schizont); + m_femalegametocytes[GametocyteStages::Stage0] = int64_t(m_femalegametocytes[GametocyteStages::Stage0] + m_IRBC_count[j] * m_gametorate * (1.0 - m_gametosexratio) * merozoitesurvival * m_config->merozoites_per_schizont); } } } @@ -435,7 +396,7 @@ namespace emodlib if ( m_IRBC_count[i] == 0 ) continue; // don't need to estimate killing if there are no IRBC of this variant to kill! // total = antibodies (major, minor, maternal) + fever + drug - double pkill = EXPCDF(-dt * ( (m_PfEMP1_antibodies[i].major->GetAntibodyConcentration() + Infection::params::non_specific_antigenicity * m_PfEMP1_antibodies[i].minor->GetAntibodyConcentration() + immunity->get_maternal_antibodies() ) * Infection::params::antibody_IRBC_killrate + fever_cytokine_killrate + drug_killrate)); + double pkill = EXPCDF(-dt * ( (m_PfEMP1_antibodies[i].major->GetAntibodyConcentration() + m_config->non_specific_antigenicity * m_PfEMP1_antibodies[i].minor->GetAntibodyConcentration() + immunity->get_maternal_antibodies() ) * m_config->antibody_IRBC_killrate + fever_cytokine_killrate + drug_killrate)); // Now here there is an interesting issue: to save massive amounts of computational time, can use a Gaussian approximation for the true binomial, but this returns a float // This is fine for large numbers of killed IRBC's, but an issue arises for small numbers @@ -443,7 +404,7 @@ namespace emodlib double tempval1 = m_IRBC_count[i] * pkill; if ( tempval1 > 0 ) // don't need to smear the killing by a random number if it is going to be zero - tempval1 = IntrahostComponent::p_rng->eGauss() * sqrt(tempval1 * (1.0 - pkill)) + tempval1; + tempval1 = m_config->rng->eGauss() * sqrt(tempval1 * (1.0 - pkill)) + tempval1; if (tempval1 < 0.5) tempval1 = 0; @@ -508,8 +469,8 @@ namespace emodlib void Infection::apply_MatureGametocyteKillProbability(float pkill) { // Gaussian approximation of binomial errors for male and female mature gametocytes - m_femalegametocytes[ GametocyteStages::Mature ] = ApplyKillProbability( pkill, m_femalegametocytes[ GametocyteStages::Mature ], IntrahostComponent::p_rng->eGauss() ); - m_malegametocytes[ GametocyteStages::Mature ] = ApplyKillProbability( pkill, m_malegametocytes[ GametocyteStages::Mature ], IntrahostComponent::p_rng->eGauss() ); + m_femalegametocytes[ GametocyteStages::Mature ] = ApplyKillProbability( pkill, m_femalegametocytes[ GametocyteStages::Mature ], m_config->rng->eGauss() ); + m_malegametocytes[ GametocyteStages::Mature ] = ApplyKillProbability( pkill, m_malegametocytes[ GametocyteStages::Mature ], m_config->rng->eGauss() ); } void Infection::malariaCheckInfectionStatus(float dt) diff --git a/include/emodlib/malaria/InfectionMalaria.h b/include/emodlib/malaria/InfectionMalaria.h index 4796c93..eb46464 100644 --- a/include/emodlib/malaria/InfectionMalaria.h +++ b/include/emodlib/malaria/InfectionMalaria.h @@ -6,7 +6,8 @@ #pragma once -#include "emodlib/ParamSet.h" +#include + #include "emodlib/utils/suids.hpp" #include "Malaria.h" @@ -20,6 +21,7 @@ namespace emodlib namespace malaria { + struct MalariaConfig; // forward declaration class Susceptibility; class Infection @@ -27,33 +29,7 @@ namespace emodlib public: - struct params - { - // TODO: emodlib#8 (boost + enums) - // static ParasiteSwitchType::Enum parasite_switch_type; - // static MalariaStrains::Enum malaria_strains; - - static float incubation_period; - static float antibody_IRBC_killrate; - static float non_specific_antigenicity; - static float MSP1_merozoite_kill; - static float gametocyte_stage_survival; - static float base_gametocyte_sexratio; - static float base_gametocyte_production; - static float antigen_switch_rate; - static float merozoites_per_hepatocyte; - static float merozoites_per_schizont; - static float RBC_destruction_multiplier; - static int n_asexual_cycles_wo_gametocytes; - - static void Configure(const ParamSet& pset); - }; - - - static suids::distributed_generator infectionSuidGenerator; - - - static Infection *Create(Susceptibility* _susceptibility, int initial_hepatocytes=1); + static Infection *Create(Susceptibility* _susceptibility, MalariaConfig* config, int initial_hepatocytes=1); void Update(float dt); @@ -96,9 +72,10 @@ namespace emodlib double m_gametosexratio; Susceptibility* immunity; + MalariaConfig* m_config; // non-owning pointer to configuration - Infection(); + Infection(MalariaConfig* config); void Initialize(Susceptibility* _susceptibility, int initial_hepatocytes); void malariaProcessHepatocytes(float dt); diff --git a/include/emodlib/malaria/IntrahostComponent.cpp b/include/emodlib/malaria/IntrahostComponent.cpp index df2057d..d76724a 100644 --- a/include/emodlib/malaria/IntrahostComponent.cpp +++ b/include/emodlib/malaria/IntrahostComponent.cpp @@ -16,53 +16,18 @@ namespace emodlib namespace malaria { - int IntrahostComponent::params::randomSeed = 0; - - int IntrahostComponent::params::max_ind_inf = 1; - - - int IntrahostComponent::params::falciparumMSPVars = DEFAULT_MSP_VARIANTS; - int IntrahostComponent::params::falciparumNonSpecTypes = DEFAULT_NONSPECIFIC_TYPES; - int IntrahostComponent::params::falciparumPfEMP1Vars = DEFAULT_PFEMP1_VARIANTS; - - // TODO: emodlib#7 (infectiousness calculations) - float IntrahostComponent::params::base_gametocyte_mosquito_survival = DEFAULT_BASE_GAMETOCYTE_MOSQUITO_SURVIVAL; - float IntrahostComponent::params::cytokine_gametocyte_inactivation = DEFAULT_CYTOKINE_GAMETOCYTE_INACTIVATION; - - std::shared_ptr IntrahostComponent::p_rng = nullptr; - - - void IntrahostComponent::params::Configure(const ParamSet& pset) - { - randomSeed = pset["Run_Number"].cast(); - IntrahostComponent::p_rng = std::shared_ptr(new PSEUDO_DES(randomSeed, 256)); - - max_ind_inf = pset["Max_Individual_Infections"].cast(); - - falciparumMSPVars = pset["Falciparum_MSP_Variants"].cast(); - falciparumNonSpecTypes = pset["Falciparum_Nonspecific_Types"].cast(); - falciparumPfEMP1Vars = pset["Falciparum_PfEMP1_Variants"].cast(); - - // TODO: emodlib#7 (infectiousness calculations) - base_gametocyte_mosquito_survival = pset["Base_Gametocyte_Mosquito_Survival_Rate"].cast(); - cytokine_gametocyte_inactivation = pset["Cytokine_Gametocyte_Inactivation"].cast(); - - Infection::params::Configure(pset["infection_params"]); - Susceptibility::params::Configure(pset["susceptibility_params"]); - } - - - IntrahostComponent::IntrahostComponent() - : susceptibility(nullptr) + IntrahostComponent::IntrahostComponent(std::shared_ptr config) + : m_config(config) + , susceptibility(nullptr) , infections() { } - IntrahostComponent* IntrahostComponent::Create() + IntrahostComponent* IntrahostComponent::Create(std::shared_ptr config) { - IntrahostComponent* ic = new IntrahostComponent(); - ic->susceptibility = Susceptibility::Create(); + IntrahostComponent* ic = new IntrahostComponent(config); + ic->susceptibility = Susceptibility::Create(config.get()); return ic; } @@ -92,8 +57,8 @@ namespace emodlib void IntrahostComponent::Challenge() { - if (infections.size() < params::max_ind_inf) { - Infection* inf = Infection::Create(susceptibility); + if (infections.size() < m_config->max_ind_inf) { + Infection* inf = Infection::Create(susceptibility, m_config.get()); infections.push_back(inf); // TODO: emodlib#2 (Max_Individual_Infections) } } @@ -134,8 +99,8 @@ namespace emodlib float IntrahostComponent::GetInfectiousness() const { float cytokines = susceptibility->get_cytokines(); - double fever_effect = Sigmoid::basic_sigmoid(params::cytokine_gametocyte_inactivation, cytokines); - return EXPCDF(-GetGametocyteDensity() * MICROLITERS_PER_BLOODMEAL * params::base_gametocyte_mosquito_survival * (1.0 - fever_effect)); + double fever_effect = Sigmoid::basic_sigmoid(m_config->cytokine_gametocyte_inactivation, cytokines); + return EXPCDF(-GetGametocyteDensity() * MICROLITERS_PER_BLOODMEAL * m_config->base_gametocyte_mosquito_survival * (1.0 - fever_effect)); } Susceptibility* IntrahostComponent::GetSusceptibility() const @@ -148,6 +113,11 @@ namespace emodlib return infections; } + std::shared_ptr IntrahostComponent::GetConfig() const + { + return m_config; + } + } } diff --git a/include/emodlib/malaria/IntrahostComponent.h b/include/emodlib/malaria/IntrahostComponent.h index 5919dec..215139e 100644 --- a/include/emodlib/malaria/IntrahostComponent.h +++ b/include/emodlib/malaria/IntrahostComponent.h @@ -9,9 +9,7 @@ #include #include -#include "emodlib/ParamSet.h" -#include "emodlib/utils/RANDOM.h" - +#include "MalariaConfig.h" #include "InfectionMalaria.h" #include "SusceptibilityMalaria.h" @@ -27,27 +25,7 @@ namespace emodlib public: - struct params - { - static int randomSeed; - - static int max_ind_inf; - - static int falciparumMSPVars; - static int falciparumNonSpecTypes; - static int falciparumPfEMP1Vars; - - // ... infectiousness calculations - static float base_gametocyte_mosquito_survival; // TODO: emodlib#7 (infectiousness calculations) - static float cytokine_gametocyte_inactivation; - - static void Configure(const ParamSet& pset); - }; - - - static std::shared_ptr p_rng; - - static IntrahostComponent* Create(); + static IntrahostComponent* Create(std::shared_ptr config); void Update(float dt); @@ -65,13 +43,16 @@ namespace emodlib Susceptibility* GetSusceptibility() const; std::list GetInfections() const; + std::shared_ptr GetConfig() const; + private: + std::shared_ptr m_config; Susceptibility* susceptibility; std::list infections; - IntrahostComponent(); + IntrahostComponent(std::shared_ptr config); }; diff --git a/include/emodlib/malaria/MalariaAntibody.cpp b/include/emodlib/malaria/MalariaAntibody.cpp index c35b7ad..b2c00b8 100644 --- a/include/emodlib/malaria/MalariaAntibody.cpp +++ b/include/emodlib/malaria/MalariaAntibody.cpp @@ -1,7 +1,7 @@ #include "MalariaAntibody.h" #include "emodlib/utils/Sigmoid.h" -#include "SusceptibilityMalaria.h" +#include "MalariaConfig.h" #define NON_TRIVIAL_ANTIBODY_THRESHOLD (0.0000001) @@ -18,8 +18,9 @@ namespace emodlib namespace malaria { - MalariaAntibody::MalariaAntibody() - : m_antigen_count(0) + MalariaAntibody::MalariaAntibody(MalariaConfig* config) + : m_config(config) + , m_antigen_count(0) , m_antigen_present(false) { } @@ -41,9 +42,9 @@ namespace emodlib } // antibody capacity decays to a medium value (.3) dropping below .4 in ~120 days from 1.0 - if ( m_antibody_capacity > Susceptibility::params::memory_level ) + if ( m_antibody_capacity > m_config->memory_level ) { - m_antibody_capacity -= ( m_antibody_capacity - Susceptibility::params::memory_level) * Susceptibility::params::hyperimmune_decay_rate * dt; + m_antibody_capacity -= ( m_antibody_capacity - m_config->memory_level) * m_config->hyperimmune_decay_rate * dt; } } @@ -52,7 +53,7 @@ namespace emodlib // allow the decay of anti-CSP concentrations greater than unity (e.g. after boosting by vaccine) if ( m_antibody_concentration > m_antibody_capacity ) { - m_antibody_concentration -= m_antibody_concentration * dt / Susceptibility::params::antibody_csp_decay_days; + m_antibody_concentration -= m_antibody_concentration * dt / m_config->antibody_csp_decay_days; } else { @@ -70,8 +71,8 @@ namespace emodlib // Let's use the MSP version of antibody growth in the base class ... void MalariaAntibody::UpdateAntibodyCapacity( float dt, float inv_uL_blood ) { - float growth_rate = Susceptibility::params::MSP1_antibody_growthrate; - float threshold = Susceptibility::params::antibody_stimulation_c50; + float growth_rate = m_config->MSP1_antibody_growthrate; + float threshold = m_config->antibody_stimulation_c50; m_antibody_capacity += growth_rate * (1.0f - m_antibody_capacity) * float(Sigmoid::basic_sigmoid( threshold, float(m_antigen_count) * inv_uL_blood)); @@ -102,9 +103,9 @@ namespace emodlib // The minor PfEMP1 version is similar but not exactly the same... void MalariaAntibodyPfEMP1Minor::UpdateAntibodyCapacity( float dt, float inv_uL_blood ) { - float min_stimulation = Susceptibility::params::antibody_stimulation_c50 * Susceptibility::params::minimum_adapted_response; - float growth_rate = Susceptibility::params::antibody_capacity_growthrate * Susceptibility::params::non_specific_growth; - float threshold = Susceptibility::params::antibody_stimulation_c50; + float min_stimulation = m_config->antibody_stimulation_c50 * m_config->minimum_adapted_response; + float growth_rate = m_config->antibody_capacity_growthrate * m_config->non_specific_growth; + float threshold = m_config->antibody_stimulation_c50; if (m_antibody_capacity <= B_CELL_PROLIFERATION_THRESHOLD) { @@ -125,9 +126,9 @@ namespace emodlib // The major PfEMP1 version is slightly different again... void MalariaAntibodyPfEMP1Major::UpdateAntibodyCapacity( float dt, float inv_uL_blood ) { - float min_stimulation = Susceptibility::params::antibody_stimulation_c50 * Susceptibility::params::minimum_adapted_response; - float growth_rate = Susceptibility::params::antibody_capacity_growthrate; - float threshold = Susceptibility::params::antibody_stimulation_c50; + float min_stimulation = m_config->antibody_stimulation_c50 * m_config->minimum_adapted_response; + float growth_rate = m_config->antibody_capacity_growthrate; + float threshold = m_config->antibody_stimulation_c50; if (m_antibody_capacity <= B_CELL_PROLIFERATION_THRESHOLD) { @@ -169,7 +170,7 @@ namespace emodlib // allow the decay of anti-CSP concentrations greater than unity (e.g. after boosting by vaccine) if ( m_antibody_concentration > m_antibody_capacity ) { - m_antibody_concentration -= m_antibody_concentration * dt / Susceptibility::params::antibody_csp_decay_days; + m_antibody_concentration -= m_antibody_concentration * dt / m_config->antibody_csp_decay_days; } else { @@ -240,33 +241,33 @@ namespace emodlib //------------------------------------------------------------------ - IMalariaAntibody* MalariaAntibodyCSP::CreateAntibody( int variant, float capacity ) + IMalariaAntibody* MalariaAntibodyCSP::CreateAntibody( MalariaConfig* config, int variant, float capacity ) { - MalariaAntibodyCSP * antibody = new MalariaAntibodyCSP(); + MalariaAntibodyCSP * antibody = new MalariaAntibodyCSP(config); antibody->Initialize( MalariaAntibodyType::CSP, variant, capacity ); return antibody; } - IMalariaAntibody* MalariaAntibodyMSP::CreateAntibody( int variant, float capacity ) + IMalariaAntibody* MalariaAntibodyMSP::CreateAntibody( MalariaConfig* config, int variant, float capacity ) { - MalariaAntibodyMSP * antibody = new MalariaAntibodyMSP(); + MalariaAntibodyMSP * antibody = new MalariaAntibodyMSP(config); antibody->Initialize( MalariaAntibodyType::MSP1, variant, capacity ); return antibody; } - IMalariaAntibody* MalariaAntibodyPfEMP1Minor::CreateAntibody( int variant, float capacity ) + IMalariaAntibody* MalariaAntibodyPfEMP1Minor::CreateAntibody( MalariaConfig* config, int variant, float capacity ) { - MalariaAntibodyPfEMP1Minor * antibody = new MalariaAntibodyPfEMP1Minor(); + MalariaAntibodyPfEMP1Minor * antibody = new MalariaAntibodyPfEMP1Minor(config); antibody->Initialize( MalariaAntibodyType::PfEMP1_minor, variant, capacity ); return antibody; } - IMalariaAntibody* MalariaAntibodyPfEMP1Major::CreateAntibody( int variant, float capacity ) + IMalariaAntibody* MalariaAntibodyPfEMP1Major::CreateAntibody( MalariaConfig* config, int variant, float capacity ) { - MalariaAntibodyPfEMP1Major * antibody = new MalariaAntibodyPfEMP1Major(); + MalariaAntibodyPfEMP1Major * antibody = new MalariaAntibodyPfEMP1Major(config); antibody->Initialize( MalariaAntibodyType::PfEMP1_major, variant, capacity ); return antibody; diff --git a/include/emodlib/malaria/MalariaAntibody.h b/include/emodlib/malaria/MalariaAntibody.h index fd705fe..b678915 100644 --- a/include/emodlib/malaria/MalariaAntibody.h +++ b/include/emodlib/malaria/MalariaAntibody.h @@ -8,6 +8,8 @@ namespace emodlib namespace malaria { + struct MalariaConfig; // forward declaration + class MalariaAntibody : public IMalariaAntibody { public: @@ -38,6 +40,8 @@ namespace emodlib virtual int GetAntibodyVariant() const override; protected: + MalariaConfig* m_config; // non-owning pointer to configuration + float m_antibody_capacity; float m_antibody_concentration; int64_t m_antigen_count; @@ -46,7 +50,7 @@ namespace emodlib MalariaAntibodyType::Enum m_antibody_type; int m_antibody_variant; - MalariaAntibody(); + MalariaAntibody(MalariaConfig* config); void Initialize( MalariaAntibodyType::Enum type, int variant, float capacity = 0, float concentration = 0 ); }; @@ -55,7 +59,8 @@ namespace emodlib class MalariaAntibodyCSP : public MalariaAntibody { public: - static IMalariaAntibody* CreateAntibody( int variant, float capacity=0.0f ); + MalariaAntibodyCSP(MalariaConfig* config) : MalariaAntibody(config) {} + static IMalariaAntibody* CreateAntibody( MalariaConfig* config, int variant, float capacity=0.0f ); virtual void UpdateAntibodyConcentration( float dt ) override; virtual void Decay( float dt ) override; }; @@ -63,20 +68,23 @@ namespace emodlib class MalariaAntibodyMSP : public MalariaAntibody { public: - static IMalariaAntibody* CreateAntibody( int variant, float capacity=0.0f ); + MalariaAntibodyMSP(MalariaConfig* config) : MalariaAntibody(config) {} + static IMalariaAntibody* CreateAntibody( MalariaConfig* config, int variant, float capacity=0.0f ); }; class MalariaAntibodyPfEMP1Minor : public MalariaAntibody { public: - static IMalariaAntibody* CreateAntibody( int variant, float capacity=0.0f ); + MalariaAntibodyPfEMP1Minor(MalariaConfig* config) : MalariaAntibody(config) {} + static IMalariaAntibody* CreateAntibody( MalariaConfig* config, int variant, float capacity=0.0f ); virtual void UpdateAntibodyCapacity( float dt, float inv_uL_blood ) override; }; class MalariaAntibodyPfEMP1Major : public MalariaAntibody { public: - static IMalariaAntibody* CreateAntibody( int variant, float capacity=0.0f ); + MalariaAntibodyPfEMP1Major(MalariaConfig* config) : MalariaAntibody(config) {} + static IMalariaAntibody* CreateAntibody( MalariaConfig* config, int variant, float capacity=0.0f ); virtual void UpdateAntibodyCapacity( float dt, float inv_uL_blood ) override; }; } diff --git a/include/emodlib/malaria/MalariaConfig.cpp b/include/emodlib/malaria/MalariaConfig.cpp new file mode 100644 index 0000000..f2b6207 --- /dev/null +++ b/include/emodlib/malaria/MalariaConfig.cpp @@ -0,0 +1,86 @@ +/** + * @file MalariaConfig.cpp + * + * @brief Malaria configuration implementation + */ + +#include "MalariaConfig.h" + +#include + + +namespace emodlib +{ + + namespace malaria + { + + MalariaConfig::MalariaConfig() + : rng(nullptr) + , suidGenerator(0, 0) + { + } + + + void MalariaConfig::Configure(const ParamSet& pset) + { + // IntrahostComponent-level params + randomSeed = pset["Run_Number"].cast(); + rng = std::make_shared(randomSeed, 256); + + max_ind_inf = pset["Max_Individual_Infections"].cast(); + + falciparumMSPVars = pset["Falciparum_MSP_Variants"].cast(); + falciparumNonSpecTypes = pset["Falciparum_Nonspecific_Types"].cast(); + falciparumPfEMP1Vars = pset["Falciparum_PfEMP1_Variants"].cast(); + + base_gametocyte_mosquito_survival = pset["Base_Gametocyte_Mosquito_Survival_Rate"].cast(); + cytokine_gametocyte_inactivation = pset["Cytokine_Gametocyte_Inactivation"].cast(); + + // Susceptibility params (from nested susceptibility_params) + const ParamSet& susc_pset = pset["susceptibility_params"]; + + memory_level = susc_pset["Antibody_Memory_Level"].cast(); + hyperimmune_decay_rate = -log((0.4f - memory_level) / (1.0f - memory_level)) / 120.0f; + MSP1_antibody_growthrate = susc_pset["Max_MSP1_Antibody_Growthrate"].cast(); + antibody_stimulation_c50 = susc_pset["Antibody_Stimulation_C50"].cast(); + antibody_capacity_growthrate = susc_pset["Antibody_Capacity_Growth_Rate"].cast(); + minimum_adapted_response = susc_pset["Min_Adapted_Response"].cast(); + non_specific_growth = susc_pset["Nonspecific_Antibody_Growth_Rate_Factor"].cast(); + antibody_csp_decay_days = susc_pset["Antibody_CSP_Decay_Days"].cast(); + + maternal_antibody_decay_rate = susc_pset["Maternal_Antibody_Decay_Rate"].cast(); + + pyrogenic_threshold = susc_pset["Pyrogenic_Threshold"].cast(); + fever_IRBC_killrate = susc_pset["Fever_IRBC_Kill_Rate"].cast(); + + erythropoiesis_anemia_effect = susc_pset["Erythropoiesis_Anemia_Effect"].cast(); + + // Infection params (from nested infection_params) + const ParamSet& inf_pset = pset["infection_params"]; + + incubation_period = inf_pset["Base_Incubation_Period"].cast(); + antibody_IRBC_killrate = inf_pset["Antibody_IRBC_Kill_Rate"].cast(); + non_specific_antigenicity = inf_pset["Nonspecific_Antigenicity_Factor"].cast(); + MSP1_merozoite_kill = inf_pset["MSP1_Merozoite_Kill_Fraction"].cast(); + gametocyte_stage_survival = inf_pset["Gametocyte_Stage_Survival_Rate"].cast(); + base_gametocyte_sexratio = inf_pset["Base_Gametocyte_Fraction_Male"].cast(); + base_gametocyte_production = inf_pset["Base_Gametocyte_Production_Rate"].cast(); + antigen_switch_rate = inf_pset["Antigen_Switch_Rate"].cast(); + merozoites_per_hepatocyte = inf_pset["Merozoites_Per_Hepatocyte"].cast(); + merozoites_per_schizont = inf_pset["Merozoites_Per_Schizont"].cast(); + RBC_destruction_multiplier = inf_pset["RBC_Destruction_Multiplier"].cast(); + n_asexual_cycles_wo_gametocytes = inf_pset["Number_Of_Asexual_Cycles_Without_Gametocytes"].cast(); + } + + + std::shared_ptr MalariaConfig::FromParamSet(const ParamSet& pset) + { + auto config = std::make_shared(); + config->Configure(pset); + return config; + } + + } + +} diff --git a/include/emodlib/malaria/MalariaConfig.h b/include/emodlib/malaria/MalariaConfig.h new file mode 100644 index 0000000..85f60da --- /dev/null +++ b/include/emodlib/malaria/MalariaConfig.h @@ -0,0 +1,112 @@ +/** + * @file MalariaConfig.h + * + * @brief Malaria configuration - instance-based parameters for thread-safe operation + */ + +#pragma once + +#include + +#include "emodlib/ParamSet.h" +#include "emodlib/utils/RANDOM.h" +#include "emodlib/utils/suids.hpp" + +#include "Malaria.h" + + +namespace emodlib +{ + + namespace malaria + { + + /** + * @brief Configuration object that owns all malaria model parameters. + * + * This replaces the static params structs in IntrahostComponent, Susceptibility, + * and Infection classes. Each IntrahostComponent instance holds a shared_ptr + * to a MalariaConfig, enabling thread-safe operation with independent configs. + */ + struct MalariaConfig + { + // ===================================================== + // IntrahostComponent-level parameters + // ===================================================== + + int randomSeed = 0; + int max_ind_inf = 1; + + int falciparumMSPVars = DEFAULT_MSP_VARIANTS; + int falciparumNonSpecTypes = DEFAULT_NONSPECIFIC_TYPES; + int falciparumPfEMP1Vars = DEFAULT_PFEMP1_VARIANTS; + + float base_gametocyte_mosquito_survival = DEFAULT_BASE_GAMETOCYTE_MOSQUITO_SURVIVAL; + float cytokine_gametocyte_inactivation = DEFAULT_CYTOKINE_GAMETOCYTE_INACTIVATION; + + // ===================================================== + // Susceptibility parameters + // ===================================================== + + float memory_level = 0.2f; + float hyperimmune_decay_rate = 0.0f; // derived from memory_level + float MSP1_antibody_growthrate = 0.02f; + float antibody_stimulation_c50 = 10.0f; + float antibody_capacity_growthrate = 0.1f; + float minimum_adapted_response = 0.02f; + float non_specific_growth = 0.5f; + float antibody_csp_decay_days = DEFAULT_ANTIBODY_CSP_DECAY_DAYS; + + float maternal_antibody_decay_rate = 0.01f; + + float pyrogenic_threshold = 1000.0f; + float fever_IRBC_killrate = DEFAULT_FEVER_IRBC_KILL_RATE; + + float erythropoiesis_anemia_effect = 3.5f; + + // ===================================================== + // Infection parameters + // ===================================================== + + float incubation_period = 7.0f; + float antibody_IRBC_killrate = DEFAULT_ANTIBODY_IRBC_KILLRATE; + float non_specific_antigenicity = DEFAULT_NON_SPECIFIC_ANTIGENICITY; + float MSP1_merozoite_kill = DEFAULT_MSP1_MEROZOITE_KILL; + float gametocyte_stage_survival = DEFAULT_GAMETOCYTE_STAGE_SURVIVAL; + float base_gametocyte_sexratio = DEFAULT_BASE_GAMETOCYTE_SEX_RATIO; + float base_gametocyte_production = DEFAULT_BASE_GAMETOCYTE_PRODUCTION; + float antigen_switch_rate = DEFAULT_ANTIGEN_SWITCH_RATE; + float merozoites_per_hepatocyte = DEFAULT_MEROZOITES_PER_HEPATOCYTE; + float merozoites_per_schizont = DEFAULT_MEROZOITES_PER_SCHIZONT; + float RBC_destruction_multiplier = DEFAULT_RBC_DESTRUCTION_MULTIPLIER; + int n_asexual_cycles_wo_gametocytes = DEFAULT_ASEXUAL_CYCLES_WITHOUT_GAMETOCYTES; + + // ===================================================== + // Instance state (owned by config) + // ===================================================== + + std::shared_ptr rng; + suids::distributed_generator suidGenerator; + + // ===================================================== + // Methods + // ===================================================== + + MalariaConfig(); + + /** + * @brief Configure from a ParamSet (Python dict/YAML). + * + * This replaces the static Configure() methods in the params structs. + */ + void Configure(const ParamSet& pset); + + /** + * @brief Create a configured MalariaConfig from a ParamSet. + */ + static std::shared_ptr FromParamSet(const ParamSet& pset); + }; + + } + +} diff --git a/include/emodlib/malaria/SusceptibilityMalaria.cpp b/include/emodlib/malaria/SusceptibilityMalaria.cpp index 8a06256..f15b827 100644 --- a/include/emodlib/malaria/SusceptibilityMalaria.cpp +++ b/include/emodlib/malaria/SusceptibilityMalaria.cpp @@ -12,6 +12,7 @@ #include "emodlib/utils/Sigmoid.h" #include "Malaria.h" +#include "MalariaConfig.h" #include "MalariaAntibody.h" @@ -21,51 +22,9 @@ namespace emodlib namespace malaria { - float Susceptibility::params::memory_level = 0.2f; - float Susceptibility::params::hyperimmune_decay_rate = 0.0f; - float Susceptibility::params::MSP1_antibody_growthrate = 0.02f; - float Susceptibility::params::antibody_stimulation_c50 = 10.0f; - float Susceptibility::params::antibody_capacity_growthrate = 0.1f; - float Susceptibility::params::minimum_adapted_response = 0.02f; - float Susceptibility::params::non_specific_growth = 0.5f; - float Susceptibility::params::antibody_csp_decay_days = DEFAULT_ANTIBODY_CSP_DECAY_DAYS; - - // TODO: emodlib#9 (maternal antibody init) + emodlib#8 (boost + enums) - // bool Susceptibility::params::enable_maternal_antibodies_transmission = false; - // MaternalAntibodiesType::Enum Susceptibility::params::maternal_antibodies_type = MaternalAntibodiesType::OFF; - // float Susceptibility::params::maternal_antibody_protection = 0.1f; - float Susceptibility::params::maternal_antibody_decay_rate = 0.01f; - - // TODO: emodlib#9 (innate heterogeneity init) + emodlib#8 (boost + enums) - // InnateImmuneVariationType::Enum Susceptibility::params::innate_immune_variation_type = InnateImmuneVariationType::NONE; - float Susceptibility::params::pyrogenic_threshold = 1000.0f; - float Susceptibility::params::fever_IRBC_killrate = DEFAULT_FEVER_IRBC_KILL_RATE; - - float Susceptibility::params::erythropoiesis_anemia_effect = 3.5f; - - - void Susceptibility::params::Configure(const ParamSet& pset) - { - memory_level = pset["Antibody_Memory_Level"].cast(); - hyperimmune_decay_rate = -log((0.4f - memory_level) / (1.0f - memory_level)) / 120.0f; // This sets the decay rate towards memory level so that the decay from antibody levels of 1 to levels of 0.4 is consistent - MSP1_antibody_growthrate = pset["Max_MSP1_Antibody_Growthrate"].cast(); - antibody_stimulation_c50 = pset["Antibody_Stimulation_C50"].cast(); - antibody_capacity_growthrate = pset["Antibody_Capacity_Growth_Rate"].cast(); - minimum_adapted_response = pset["Min_Adapted_Response"].cast(); - non_specific_growth = pset["Nonspecific_Antibody_Growth_Rate_Factor"].cast(); - antibody_csp_decay_days = pset["Antibody_CSP_Decay_Days"].cast(); - - maternal_antibody_decay_rate = pset["Maternal_Antibody_Decay_Rate"].cast(); - - pyrogenic_threshold = pset["Pyrogenic_Threshold"].cast(); - fever_IRBC_killrate = pset["Fever_IRBC_Kill_Rate"].cast(); - - erythropoiesis_anemia_effect = pset["Erythropoiesis_Anemia_Effect"].cast(); - } - - - Susceptibility::Susceptibility() - : age(0) + Susceptibility::Susceptibility(MalariaConfig* config) + : m_config(config) + , age(0) , m_antigenic_flag(0) , m_maternal_antibody_strength(0) @@ -88,9 +47,9 @@ namespace emodlib } - Susceptibility* Susceptibility::Create() + Susceptibility* Susceptibility::Create(MalariaConfig* config) { - Susceptibility *newsusceptibility = new Susceptibility(); + Susceptibility *newsusceptibility = new Susceptibility(config); newsusceptibility->Initialize(); return newsusceptibility; @@ -108,12 +67,12 @@ namespace emodlib // Track individual pyrogenic thresholds + fever killing rates as instance variables // TODO: emodlib#9 (innate heterogeneity init) - m_ind_pyrogenic_threshold = params::pyrogenic_threshold; - m_ind_fever_kill_rate = params::fever_IRBC_killrate; + m_ind_pyrogenic_threshold = m_config->pyrogenic_threshold; + m_ind_fever_kill_rate = m_config->fever_IRBC_killrate; // TODO: emodlib#9 (maternal antibody init) - m_CSP_antibody = MalariaAntibodyCSP::CreateAntibody(0); + m_CSP_antibody = MalariaAntibodyCSP::CreateAntibody(m_config, 0); // MSP + PfEMP1 antibodies are added upon infection } @@ -121,7 +80,6 @@ namespace emodlib IMalariaAntibody* Susceptibility::RegisterAntibody(MalariaAntibodyType::Enum type, int variant, float capacity) { std::vector *variant_vector; - IMalariaAntibody* (*typed_create_antibody)(int,float); switch( type ) { @@ -130,17 +88,14 @@ namespace emodlib case MalariaAntibodyType::MSP1: variant_vector = &m_active_MSP_antibodies; - typed_create_antibody = MalariaAntibodyMSP::CreateAntibody; break; case MalariaAntibodyType::PfEMP1_minor: variant_vector = &m_active_PfEMP1_minor_antibodies; - typed_create_antibody = MalariaAntibodyPfEMP1Minor::CreateAntibody; break; case MalariaAntibodyType::PfEMP1_major: variant_vector = &m_active_PfEMP1_major_antibodies; - typed_create_antibody = MalariaAntibodyPfEMP1Major::CreateAntibody; break; default: @@ -160,7 +115,20 @@ namespace emodlib if (antibody == nullptr) // make a new antibody if it hasn't been created yet { - antibody = typed_create_antibody(variant, capacity); + switch( type ) + { + case MalariaAntibodyType::MSP1: + antibody = MalariaAntibodyMSP::CreateAntibody(m_config, variant, capacity); + break; + case MalariaAntibodyType::PfEMP1_minor: + antibody = MalariaAntibodyPfEMP1Minor::CreateAntibody(m_config, variant, capacity); + break; + case MalariaAntibodyType::PfEMP1_major: + antibody = MalariaAntibodyPfEMP1Major::CreateAntibody(m_config, variant, capacity); + break; + default: + break; + } variant_vector->push_back(antibody); } @@ -192,10 +160,10 @@ namespace emodlib recalculateBloodCapacity(age); // Red blood cell dynamics - if (Susceptibility::params::erythropoiesis_anemia_effect > 0) + if (m_config->erythropoiesis_anemia_effect > 0) { // This is the amount of "erythropoietin", assume absolute amounts of erythropoietin correlate linearly with absolute increases in hemoglobin - float anemia_erythropoiesis_multiplier = exp( Susceptibility::params::erythropoiesis_anemia_effect * (1 - get_RBC_availability()) ); + float anemia_erythropoiesis_multiplier = exp( m_config->erythropoiesis_anemia_effect * (1 - get_RBC_availability()) ); m_RBC = int64_t(m_RBC - (m_RBC * .00833 - m_RBCproduction * anemia_erythropoiesis_multiplier) * dt); // *.00833 ==/120 (AVERAGE_RBC_LIFESPAN) } else @@ -211,7 +179,7 @@ namespace emodlib m_parasite_density = 0; // this is accumulated in updateImmunityPfEMP1Minor // decay maternal antibodies - m_maternal_antibody_strength -= dt * m_maternal_antibody_strength * Susceptibility::params::maternal_antibody_decay_rate; + m_maternal_antibody_strength -= dt * m_maternal_antibody_strength * m_config->maternal_antibody_decay_rate; if ( m_maternal_antibody_strength < 0 ) { m_maternal_antibody_strength = 0; } // antibody capacities increase and antibodies released if antigen present, only process if antigens are present at all @@ -477,6 +445,11 @@ namespace emodlib { m_ind_fever_kill_rate = _rate; } + + MalariaConfig* Susceptibility::GetConfig() const + { + return m_config; + } } } diff --git a/include/emodlib/malaria/SusceptibilityMalaria.h b/include/emodlib/malaria/SusceptibilityMalaria.h index 128fc1d..6923261 100644 --- a/include/emodlib/malaria/SusceptibilityMalaria.h +++ b/include/emodlib/malaria/SusceptibilityMalaria.h @@ -6,7 +6,7 @@ #pragma once -#include "emodlib/ParamSet.h" +#include #include "MalariaEnums.h" #include "IMalariaAntibody.h" @@ -18,43 +18,14 @@ namespace emodlib namespace malaria { + struct MalariaConfig; // forward declaration + class Susceptibility { public: - struct params - { - // Used in MalariaAntibody for boost-decay functions - static float memory_level; - static float hyperimmune_decay_rate; - static float MSP1_antibody_growthrate; - static float antibody_stimulation_c50; - static float antibody_capacity_growthrate; - static float minimum_adapted_response; - static float non_specific_growth; - static float antibody_csp_decay_days; - - // Used in Susceptibility for: - // ...maternal protection - // static bool enable_maternal_antibodies_transmission; - // static MaternalAntibodiesType::Enum maternal_antibodies_type; // TODO: emodlib#9 (innate init) - // static float maternal_antibody_protection; - static float maternal_antibody_decay_rate; - - // ...innate immunity effects - // static InnateImmuneVariationType::Enum innate_immune_variation_type; // TODO: emodlib#9 (innate init) - static float pyrogenic_threshold; - static float fever_IRBC_killrate; - - // ... red blood cell effects - static float erythropoiesis_anemia_effect; - - static void Configure(const ParamSet& pset); - }; - - - static Susceptibility *Create(); + static Susceptibility *Create(MalariaConfig* config); IMalariaAntibody* RegisterAntibody(MalariaAntibodyType::Enum type, int variant, float capacity=0.0f); void UpdateActiveAntibody( pfemp1_antibody_t &pfemp1_variant, int minor_variant, int major_variant ); void remove_RBCs(int64_t infectedAsexual, int64_t infectedGametocytes, double RBC_destruction_multiplier); @@ -84,8 +55,12 @@ namespace emodlib float get_fever_kill_rate() const; void set_fever_kill_rate(float _rate); + MalariaConfig* GetConfig() const; + private: + MalariaConfig* m_config; // non-owning pointer to configuration + float age; // TODO: emodlib#10 (demographic components) // containers for antibody objects @@ -110,7 +85,7 @@ namespace emodlib float m_parasite_density; - Susceptibility(); + Susceptibility(MalariaConfig* config); void Initialize(); // TODO: emodlib#9 (innate init) + emodlib#10 (demographic/transmission components) void recalculateBloodCapacity( float _age ); diff --git a/src/emodlib/malaria/__init__.py b/src/emodlib/malaria/__init__.py index bdd7931..1656b11 100644 --- a/src/emodlib/malaria/__init__.py +++ b/src/emodlib/malaria/__init__.py @@ -1,23 +1,82 @@ import os -from .._emodlib_py.malaria import Infection, IntrahostComponent, Susceptibility -from ..params import Params, set_params, update_params +import yaml + +from .._emodlib_py.malaria import ( + Infection, + IntrahostComponent, + MalariaConfig, + Susceptibility, +) +from ..params import Params, deep_update + + +def _config_yaml_property(self): + """Return YAML string representation of the config's source params.""" + return yaml.dump(dict(self._params), default_flow_style=False, sort_keys=False) + + +def _config_update(self, params): + """ + Update this config with new parameters (nested merge on top of current). + + This is the instance-based equivalent of the old update_params() method. + Modifies the config in-place and returns self for chaining. + + Args: + params: Dict of parameters to merge on top of current config. + Can be nested dict matching config.yml structure. + + Returns: + self (for method chaining) + + Example: + config = create_config({'Run_Number': 1}) + config.update({'infection_params': {'Antigen_Switch_Rate': 1e-8}}) + """ + merged = deep_update(self._params, params) + self.configure(merged) + self._params = merged + return self + + +# Add yaml property and update method to MalariaConfig +MalariaConfig.yaml = property(_config_yaml_property) +MalariaConfig.update = _config_update def params_from_default_file(): + """Load default parameters from the bundled config.yml file.""" path = os.path.join(os.path.realpath(os.path.dirname(__file__)), "config.yml") return Params.from_yaml(path) +# Store default params at class level for reference IntrahostComponent.default_params = params_from_default_file() -# monkey-patch set_params + update_params -# nested on top of defaults + current values, respectively -IntrahostComponent.set_params = set_params -IntrahostComponent.update_params = update_params -# initialize default parameters -IntrahostComponent.set_params() +def create_config(params=None): + """ + Create a new MalariaConfig instance. + + Args: + params: Optional dict of parameters to override defaults. + Can be a flat dict or nested dict matching config.yml structure. + + Returns: + A configured MalariaConfig instance. + + Example: + config = create_config({'Run_Number': 42}) + ic = IntrahostComponent.create(config) + """ + if params is None: + merged = Params(IntrahostComponent.default_params) + else: + merged = deep_update(IntrahostComponent.default_params, params) + cfg = MalariaConfig.from_params(merged) + cfg._params = merged # Store for yaml property + return cfg -__all__ = ["IntrahostComponent", "Susceptibility", "Infection"] +__all__ = ["IntrahostComponent", "Susceptibility", "Infection", "MalariaConfig", "create_config"] diff --git a/src/malaria.cpp b/src/malaria.cpp index d26a9ae..3363bf4 100644 --- a/src/malaria.cpp +++ b/src/malaria.cpp @@ -6,6 +6,7 @@ #include "pybind11/pybind11.h" #include "pybind11/stl.h" +#include "emodlib/malaria/MalariaConfig.h" #include "emodlib/malaria/IntrahostComponent.h" #include "emodlib/malaria/MalariaAntibody.h" @@ -51,16 +52,64 @@ void add_malaria_bindings(py::module& m) { using namespace py::literals; - // ==== Binding of the intrahost component ==== // - py::class_ (m, "IntrahostComponent") + // ==== Binding of the MalariaConfig ==== // + py::class_> (m, "MalariaConfig", py::dynamic_attr()) + + .def(py::init<>()) - .def_static("create", &IntrahostComponent::Create) + .def("configure", &MalariaConfig::Configure, + "Configure from a ParamSet dictionary", + "pset"_a) - .def_static("_configure_from_params", - &IntrahostComponent::params::Configure, - "Configure the IntrahostComponent params from a ParamSet dictionary", + .def_static("from_params", &MalariaConfig::FromParamSet, + "Create a configured MalariaConfig from a ParamSet dictionary", "pset"_a) + // IntrahostComponent-level params + .def_readwrite("random_seed", &MalariaConfig::randomSeed) + .def_readwrite("max_ind_inf", &MalariaConfig::max_ind_inf) + .def_readwrite("falciparum_MSP_variants", &MalariaConfig::falciparumMSPVars) + .def_readwrite("falciparum_nonspecific_types", &MalariaConfig::falciparumNonSpecTypes) + .def_readwrite("falciparum_PfEMP1_variants", &MalariaConfig::falciparumPfEMP1Vars) + .def_readwrite("base_gametocyte_mosquito_survival", &MalariaConfig::base_gametocyte_mosquito_survival) + .def_readwrite("cytokine_gametocyte_inactivation", &MalariaConfig::cytokine_gametocyte_inactivation) + + // Susceptibility params + .def_readwrite("memory_level", &MalariaConfig::memory_level) + .def_readwrite("hyperimmune_decay_rate", &MalariaConfig::hyperimmune_decay_rate) + .def_readwrite("MSP1_antibody_growthrate", &MalariaConfig::MSP1_antibody_growthrate) + .def_readwrite("antibody_stimulation_c50", &MalariaConfig::antibody_stimulation_c50) + .def_readwrite("antibody_capacity_growthrate", &MalariaConfig::antibody_capacity_growthrate) + .def_readwrite("minimum_adapted_response", &MalariaConfig::minimum_adapted_response) + .def_readwrite("non_specific_growth", &MalariaConfig::non_specific_growth) + .def_readwrite("antibody_csp_decay_days", &MalariaConfig::antibody_csp_decay_days) + .def_readwrite("maternal_antibody_decay_rate", &MalariaConfig::maternal_antibody_decay_rate) + .def_readwrite("pyrogenic_threshold", &MalariaConfig::pyrogenic_threshold) + .def_readwrite("fever_IRBC_killrate", &MalariaConfig::fever_IRBC_killrate) + .def_readwrite("erythropoiesis_anemia_effect", &MalariaConfig::erythropoiesis_anemia_effect) + + // Infection params + .def_readwrite("incubation_period", &MalariaConfig::incubation_period) + .def_readwrite("antibody_IRBC_killrate", &MalariaConfig::antibody_IRBC_killrate) + .def_readwrite("non_specific_antigenicity", &MalariaConfig::non_specific_antigenicity) + .def_readwrite("MSP1_merozoite_kill", &MalariaConfig::MSP1_merozoite_kill) + .def_readwrite("gametocyte_stage_survival", &MalariaConfig::gametocyte_stage_survival) + .def_readwrite("base_gametocyte_sexratio", &MalariaConfig::base_gametocyte_sexratio) + .def_readwrite("base_gametocyte_production", &MalariaConfig::base_gametocyte_production) + .def_readwrite("antigen_switch_rate", &MalariaConfig::antigen_switch_rate) + .def_readwrite("merozoites_per_hepatocyte", &MalariaConfig::merozoites_per_hepatocyte) + .def_readwrite("merozoites_per_schizont", &MalariaConfig::merozoites_per_schizont) + .def_readwrite("RBC_destruction_multiplier", &MalariaConfig::RBC_destruction_multiplier) + .def_readwrite("n_asexual_cycles_wo_gametocytes", &MalariaConfig::n_asexual_cycles_wo_gametocytes); + + + // ==== Binding of the intrahost component ==== // + py::class_ (m, "IntrahostComponent") + + .def_static("create", &IntrahostComponent::Create, + "Create an IntrahostComponent with the given configuration", + "config"_a) + .def("update", &IntrahostComponent::Update, "Update the intrahost model state by dt", @@ -94,12 +143,9 @@ void add_malaria_bindings(py::module& m) { py::class_ (m, "Susceptibility") - .def_static("create", &Susceptibility::Create) - - .def_static("configure", - &Susceptibility::params::Configure, - "Configure the Susceptibility params from a ParamSet dictionary", - "pset"_a) + .def_static("create", &Susceptibility::Create, + "Create a Susceptibility with the given configuration", + "config"_a) .def("update", &Susceptibility::Update, @@ -124,13 +170,8 @@ void add_malaria_bindings(py::module& m) { py::class_ (m, "Infection") .def_static("create", &Infection::Create, - "Create an Infection object with pointer to Susceptibility", - "susceptibility"_a, "hepatocytes"_a=1) - - .def_static("configure", - &Infection::params::Configure, - "Configure the Infection params from a ParamSet dictionary", - "pset"_a) + "Create an Infection object with pointer to Susceptibility and configuration", + "susceptibility"_a, "config"_a, "hepatocytes"_a=1) .def("update", &Infection::Update, diff --git a/tests/test_host.py b/tests/test_host.py index c191883..0064839 100644 --- a/tests/test_host.py +++ b/tests/test_host.py @@ -3,7 +3,7 @@ import pytest import yaml -from emodlib.malaria import IntrahostComponent +from emodlib.malaria import IntrahostComponent, create_config def describe(c, t=None): @@ -34,31 +34,29 @@ def test_config(): print("Load default model parameters...\n") params = params_from_default_file() - assert IntrahostComponent.params["Run_Number"] == params["Run_Number"] + # Test that we can create a config and modify it + config = create_config() + assert config.random_seed == params["Run_Number"] run_number = 123456 - params["Run_Number"] = run_number + config2 = create_config({'Run_Number': run_number}) + assert config2.random_seed == run_number - print("Update nested parameters...") - IntrahostComponent.update_params( - dict(Run_Number=run_number, infection_params=dict(Merozoites_Per_Schizont=10)) - ) - assert IntrahostComponent.params["Run_Number"] == run_number - assert ( - IntrahostComponent.params["infection_params"]["Merozoites_Per_Schizont"] == 10 - ) - assert ( - IntrahostComponent.params["susceptibility_params"]["Antibody_CSP_Decay_Days"] - == params["susceptibility_params"]["Antibody_CSP_Decay_Days"] - ) + # Test that modifying merozoites_per_schizont works + config3 = create_config({ + 'Run_Number': run_number, + 'infection_params': {'Merozoites_Per_Schizont': 10} + }) + assert config3.random_seed == run_number + assert config3.merozoites_per_schizont == 10 def test_bindings(): - print("Set default parameters...") - IntrahostComponent.set_params() + print("Create config...") + config = create_config() print("Create...") - ic = IntrahostComponent.create() + ic = IntrahostComponent.create(config) print("Challenge...") ic.challenge() @@ -90,11 +88,11 @@ def test_bindings(): def test_infection_clearance(): - print("Set default parameters...") - IntrahostComponent.set_params() + print("Create config...") + config = create_config() print("Create...") - ic = IntrahostComponent.create() + ic = IntrahostComponent.create(config) print("Update...") n_infections_cache = 0 @@ -119,11 +117,11 @@ def test_max_infections(): print("Load default model parameters...\n") params = params_from_default_file() - print("Set default parameters...") - IntrahostComponent.set_params() + print("Create config...") + config = create_config() print("Create...") - ic = IntrahostComponent.create() + ic = IntrahostComponent.create(config) print("Challenge...") for i in range(10): diff --git a/tests/test_infection.py b/tests/test_infection.py index bee9267..bf38c2b 100644 --- a/tests/test_infection.py +++ b/tests/test_infection.py @@ -1,34 +1,34 @@ import pytest -from emodlib.malaria import Infection, IntrahostComponent, Susceptibility +from emodlib.malaria import Infection, IntrahostComponent, Susceptibility, create_config @pytest.fixture -def params(): - return IntrahostComponent.params +def config(): + return create_config() @pytest.fixture -def default_susceptibility(params): - yield Susceptibility.create() +def default_susceptibility(config): + yield Susceptibility.create(config) @pytest.fixture -def default_infection(default_susceptibility): - yield Infection.create(susceptibility=default_susceptibility, hepatocytes=1) +def default_infection(default_susceptibility, config): + yield Infection.create(susceptibility=default_susceptibility, config=config, hepatocytes=1) -def test_infection(default_infection, params): +def test_infection(default_infection, config): infs = [default_infection for _ in range(10)] msp_types = [i.msp_type for i in infs] print(msp_types) - assert all([t < params["Falciparum_MSP_Variants"] for t in msp_types]) + assert all([t < config.falciparum_MSP_variants for t in msp_types]) major_types = infs[0].pfemp1_major_types print(major_types) assert len(major_types) == 50 - assert all([t < params["Falciparum_PfEMP1_Variants"] for t in major_types]) + assert all([t < config.falciparum_PfEMP1_variants for t in major_types]) def test_antibody(default_infection): diff --git a/tests/test_susceptibility.py b/tests/test_susceptibility.py index 8075b51..10b2b91 100644 --- a/tests/test_susceptibility.py +++ b/tests/test_susceptibility.py @@ -1,10 +1,11 @@ import pytest -from emodlib.malaria import IntrahostComponent, Susceptibility +from emodlib.malaria import IntrahostComponent, Susceptibility, create_config def test_aging(): - s = Susceptibility.create() + config = create_config() + s = Susceptibility.create(config) s.age = 30 print("immune system age = %d days" % s.age) @@ -16,7 +17,8 @@ def test_aging(): def test_maternal_antibodies(): - s = Susceptibility.create() + config = create_config() + s = Susceptibility.create(config) s.maternal_antibody_strength = 0.8 @@ -31,10 +33,10 @@ def test_maternal_antibodies(): def test_immune_init(): - print("Set default parameters...") - IntrahostComponent.set_params() + print("Create config...") + config = create_config() - ic = IntrahostComponent.create() + ic = IntrahostComponent.create(config) s = ic.susceptibility print("immune system age = %d days" % s.age) diff --git a/tests/test_thread_safety.py b/tests/test_thread_safety.py new file mode 100644 index 0000000..59a717a --- /dev/null +++ b/tests/test_thread_safety.py @@ -0,0 +1,106 @@ +""" +Thread-safety tests for the malaria intrahost model. + +These tests verify that multiple simulations can run in parallel with independent +configurations without interfering with each other. +""" + +import concurrent.futures + +import pytest + +from emodlib.malaria import IntrahostComponent, create_config + + +def run_simulation(seed, n_steps=100): + """Run a simulation with the given seed and return results.""" + config = create_config({'Run_Number': seed}) + ic = IntrahostComponent.create(config) + ic.challenge() + + results = [] + for _ in range(n_steps): + ic.update(dt=1) + results.append(ic.parasite_density) + + return results + + +def test_parallel_simulations_are_deterministic(): + """ + Run same simulation twice sequentially - should get identical results + since they use the same seed. + """ + results1 = run_simulation(seed=42) + results2 = run_simulation(seed=42) + + assert results1 == results2, "Sequential runs with same seed should be deterministic" + + +def test_parallel_simulations_are_independent(): + """ + Run same simulation twice in parallel with same seed. + With instance-based config, both should get identical results. + """ + with concurrent.futures.ThreadPoolExecutor(max_workers=2) as executor: + future1 = executor.submit(run_simulation, seed=42) + future2 = executor.submit(run_simulation, seed=42) + results1 = future1.result() + results2 = future2.result() + + assert results1 == results2, "Parallel runs with same seed should produce identical results" + + +def test_different_seeds_produce_different_results(): + """Different seeds should produce different results.""" + results1 = run_simulation(seed=1) + results2 = run_simulation(seed=2) + + assert results1 != results2, "Different seeds should produce different results" + + +def test_different_configs_are_isolated(): + """ + Different configs should not affect each other. + This is the key thread-safety test. + """ + config1 = create_config({'Run_Number': 1, 'Max_Individual_Infections': 3}) + config2 = create_config({'Run_Number': 2, 'Max_Individual_Infections': 10}) + + ic1 = IntrahostComponent.create(config1) + ic2 = IntrahostComponent.create(config2) + + # Challenge both many times + for _ in range(20): + ic1.challenge() + ic2.challenge() + + # Verify each uses its own max_ind_inf config + assert ic1.n_infections <= 3, f"ic1 should have at most 3 infections, got {ic1.n_infections}" + assert ic2.n_infections <= 10, f"ic2 should have at most 10 infections, got {ic2.n_infections}" + + # Since ic2 has a higher limit and same number of challenges, it should have more infections + # (probabilistically, though not guaranteed - let's just check the bounds) + print(f"ic1 infections: {ic1.n_infections}, ic2 infections: {ic2.n_infections}") + + +def test_multiple_parallel_simulations(): + """Run many simulations in parallel.""" + num_simulations = 10 + seeds = list(range(num_simulations)) + + with concurrent.futures.ThreadPoolExecutor(max_workers=4) as executor: + futures = [executor.submit(run_simulation, seed=s, n_steps=50) for s in seeds] + results = [f.result() for f in futures] + + # All results should be different (different seeds) + for i, r1 in enumerate(results): + for j, r2 in enumerate(results): + if i != j: + assert r1 != r2, f"Results {i} and {j} should be different" + + print(f"Successfully ran {num_simulations} parallel simulations") + + +if __name__ == "__main__": + pytest.main(["-vv", "-s", __file__]) From 9b194c596a4a5e07e10c6dad93f57aabd6901af7 Mon Sep 17 00:00:00 2001 From: Edward Wenger Date: Tue, 16 Dec 2025 19:35:30 -0800 Subject: [PATCH 2/2] cleanup unused dead code --- src/emodlib/params.py | 16 ---------------- 1 file changed, 16 deletions(-) diff --git a/src/emodlib/params.py b/src/emodlib/params.py index ac7d73e..cd491ad 100644 --- a/src/emodlib/params.py +++ b/src/emodlib/params.py @@ -41,19 +41,3 @@ def deep_update( else: updated_mapping[k] = v return Params(updated_mapping) - - -@classmethod -def update_params(cls, params={}): - """Nested update on top of current parameters""" - cfg = deep_update(cls.params, params) - cls.params = cfg - cls._configure_from_params(cfg) # setting static variables in bound C++ classes - - -@classmethod -def set_params(cls, params={}): - """Nested update on top of default parameters""" - cfg = deep_update(cls.default_params, params) - cls.params = cfg - cls._configure_from_params(cfg) # setting static variables in bound C++ classes