Skip to content

Commit

Permalink
Close issue 62, change names input parameters 'draw_functions' to 'pa…
Browse files Browse the repository at this point in the history
…rameters' and 'samples' instead of 'param_dict' and 'samples_dict' (#73)
  • Loading branch information
twallema authored Jul 12, 2024
1 parent 2e69d7b commit f7d8ec4
Show file tree
Hide file tree
Showing 19 changed files with 229 additions and 133 deletions.
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -64,7 +64,7 @@ The [enzyme kinetics](enzyme_kinetics.md) and [influenza 17-18](influenza_1718.m
- Version 0.2.4 (2023-12-04, PR #62)
> Validated the use of Python 3.11. Efficiency gains in simulation of jump processes. Ommitted dependency on Numba. All changes related to publishing our software manuscript in Journal of Computational Science. Improved nomenclature in model defenition.
- IN PROGRESS: 0.2.5
> Validated the use of Python 3.12. Validated pySODM on macOS Sonoma 14.5.
> Validated the use of Python 3.12. Validated pySODM on macOS Sonoma 14.5. Input arguments of 'draw functions' must now be named 'parameters' and 'samples' instead of 'param_dict' and 'samples_dict'.
- Version 0.1 (2022-12-23, PR #14)
> Application pySODM to three use cases. Documentation website. Unit tests for ODE, JumpProcess and calibration.
- Version 0.1.1 (2023-01-09, PR #20)
Expand Down
30 changes: 15 additions & 15 deletions docs/enzyme_kinetics.md
Original file line number Diff line number Diff line change
Expand Up @@ -251,14 +251,14 @@ Finally, we can use the *draw functions* to propagate the parameter samples in o
```python
if __name__ == '__main__':

def draw_fcn(param_dict, samples_dict):
def draw_fcn(parameters, samples):
# Always draw correlated samples at the SAME INDEX!
idx, param_dict['Vf_Ks'] = random.choice(list(enumerate(samples_dict['Vf_Ks'])))
param_dict['R_AS'] = samples_dict['R_AS'][idx]
param_dict['R_AW'] = samples_dict['R_AW'][idx]
param_dict['R_Es'] = samples_dict['R_Es'][idx]
param_dict['K_eq'] = samples_dict['K_eq'][idx]
return param_dict
idx, parameters['Vf_Ks'] = random.choice(list(enumerate(samples['Vf_Ks'])))
parameters['R_AS'] = samples['R_AS'][idx]
parameters['R_AW'] = samples['R_AW'][idx]
parameters['R_Es'] = samples['R_Es'][idx]
parameters['K_eq'] = samples['K_eq'][idx]
return parameters

# Loop over datasets
for i,df in enumerate(data):
Expand Down Expand Up @@ -425,16 +425,16 @@ Next, we load our previously obtained samples of the kinetic parameters and defi
```python
# Load samples
f = open(os.path.join(os.path.dirname(__file__),'data/username_SAMPLES_2023-06-07.json'))
samples_dict = json.load(f)
samples = json.load(f)

# Define draw function
def draw_fcn(param_dict, samples_dict):
idx, param_dict['Vf_Ks'] = random.choice(list(enumerate(samples_dict['Vf_Ks'])))
param_dict['R_AS'] = samples_dict['R_AS'][idx]
param_dict['R_AW'] = samples_dict['R_AW'][idx]
param_dict['R_Es'] = samples_dict['R_Es'][idx]
param_dict['K_eq'] = samples_dict['K_eq'][idx]
return param_dict
def draw_fcn(parameters, samples):
idx, parameters['Vf_Ks'] = random.choice(list(enumerate(samples['Vf_Ks'])))
parameters['R_AS'] = samples['R_AS'][idx]
parameters['R_AW'] = samples['R_AW'][idx]
parameters['R_Es'] = samples['R_Es'][idx]
parameters['K_eq'] = samples['K_eq'][idx]
return parameters
```

A first experiment with the continuous flow reactor was performed using a reaction mixture containing 30 mM D-glucose, 60 mM lauric acid and 28 mM water. The reactants were pumped through the reactor at a flow rate of 0.2 mL/min, resulting in an average residence time of 13.5 minutes. After ten retention times, when the outlet concentration had stabilized, three samples were withdrawn at the outlet. Then, the reactor was cut short by 0.10 m, and again three samples were withdrawn at the outlet. In this way, the reactant profile acrosss the reactor length was obtained. I omit the code to replicate the following figures from this documentation as no new concepts are introduced beyond this point. Our packed-bed model does an adequate job at describing the data.
Expand Down
15 changes: 10 additions & 5 deletions docs/guidelines.md
Original file line number Diff line number Diff line change
Expand Up @@ -51,23 +51,28 @@ Always keep in mind additional dependencies lower `pySODM`'s life expectancy. Co

### Tests

Three test scripts are defined in `~/src/tests/`: 1) `test_ODE.py` to test the initializiation and simulation of ODE models, 2) `test_JumpProcess.py`, similar but for stochastic jump process models and 3) `test_calibration.py`, to test the common calibration workflow. Testing these routines is usefull to verify that modifications made to the `pySODM` code do not break code. Or alternatively, if the code does break, to see where it breaks. The test routines must pass when performing a PR to Github. To run a test locally,
Three unit tests are defined in `~/src/tests/`: 1) `test_ODE.py` to test the initializiation and simulation of ODE models, 2) `test_JumpProcess.py`, to test the initialization and simulation of stochastic jump process models and 3) `test_calibration.py`, to test the calibration workflow. Passing these routines is mandatory to perform a pull request to pySODM's GitHub repository. It verifies that modifications made to the `pySODM` code do not break any of its features. Or alternatively, unit tests can be exploited to see where code breaks. To run the tests locally,
```
conda activate pySODM
conda activate PYSODM
pytest test_calibration.py
```
Make sure you have the development dependencies installed to use `pytest` (to do this install pySODM using the command `pip install -e ".[develop]"`)

### Documentation

#### Website

Documentation consists of both the technical matter about the code as well as background information on the models. To keep these up to date and centralized, we use [Sphinx](https://www.sphinx-doc.org/en/master/) which enables us to keep the documentation together on a website.

The Sphinx setup provides the usage of both `.rst` file, i.e. [restructuredtext](https://docutils.sourceforge.io/docs/ref/rst/restructuredtext.html) as well as `.md` files, i.e. [Markdown](https://www.markdownguide.org/basic-syntax/). The latter is generally experienced as easier to write, while the former provides more advanced functionalities. Existing pages of the documentation reside in `~/docs/` and can be adjusted directly. In case you want to build the documentation locally, make sure you have the development dependencies installed (`pip install -e ".[develop]"`) to run the sphinx build script. To build locally, use,
```
The Sphinx setup provides the usage of both `.rst` file, i.e. [restructuredtext](https://docutils.sourceforge.io/docs/ref/rst/restructuredtext.html) as well as `.md` files, i.e. [Markdown](https://www.markdownguide.org/basic-syntax/). The latter is generally experienced as easier to write, while the former provides more advanced functionalities. Existing pages of the documentation reside in `~/docs/` and can be adjusted directly. In case you want to build the documentation locally, make sure you have the development dependencies installed (`pip install -e ".[develop]"`) to run the sphinx build script. To build locally, run the following command in the home directory,
```bash
python setup.py build_sphinx
```
When you want to create a new page, make sure to add the page to the `index.rst` in order to make the page part of the website. The website is build automatically using [github actions](https://github.com/twallema/pySODM/blob/master/.github/workflows/deploy.yml#L22-L24) and the output is deployed to [https://twallema.github.io/pySODM/](https://twallema.github.io/pySODM/). The resulting html-website is created in the directory `build/html`. Double click any of the `html` file in the folder to open the website in your browser (no server required).
which compiles the html-website in the directory `build/html`. Double click any of the `html` file in the folder to open the website in your browser (no server required). Or alternatively, run the following command in the 'docs' directory,
```bash
make html
```
which compiles the website in the directory `docs/_build`. When you want to create a new page, make sure to add the page to the `index.rst` in order to make the page part of the website. The website is build automatically using [github actions](https://github.com/twallema/pySODM/blob/master/.github/workflows/deploy.yml#L22-L24) and the output is deployed to [https://twallema.github.io/pySODM/](https://twallema.github.io/pySODM/).

#### Docstrings

Expand Down
10 changes: 5 additions & 5 deletions docs/models.md
Original file line number Diff line number Diff line change
Expand Up @@ -32,14 +32,14 @@ Upon intialization of the model class, the following arguments must be provided.

* **sim(time, N=1, draw_function=None, samples=None, processes=None, method='RK23', rtol=1e-3, tau=None, output_timestep=1, warmup=0)**

Simulate a model over a given time period using `scipy.integrate.solve_ivp()`. Can optionally perform `N` repeated simulations. Can change the values of model parameters at every repeated simulation by drawing samples from a dictionary `samples` using a function `draw_function`.
Simulate a model over a given time period using `scipy.integrate.solve_ivp()`. Can optionally perform `N` repeated simulations. Can change the values of model parameters at every repeated simulation by drawing parameter samples from a dictionary `samples` using a function `draw_function`.

**Parameters**

* **time** - (int/float or list) - The start and stop "time" for the simulation run. Three possible inputs: 1) int/float, 2) list of int/float of type `[start_time, stop_time]`, 3) list of pd.Timestamp or str of type `[start_date, stop_date]`.
* **N** - (int) - optional - Number of repeated simulations to perform.
* **samples** - (dict) - optional - Dictionary containing samples of model parameters. Obligatory input to a *draw function*.
* **draw_function** - (function) - optional - A function consisting of two inputs: the model parameters dictionary `param_dict`, the previously documented samples dictionary `samples_dict`. Function must return the model parameters dictionary `param_dict`. Function can be used to update a model parameter's value during every repeated simulation. Usefull to propagate MCMC samples of model parameters or perform sensitivity analyses.
* **samples** - (dict) - optional - Dictionary containing samples of the distribution of model parameters. Obligatory input to a *draw function*. Advanced users: dictionary can contain any arbitrary input to a *draw function*.
* **draw_function** - (function) - optional - A function consisting of two inputs: the model parameters dictionary `parameters`, the previously documented samples dictionary `samples`. Function must return the model parameters dictionary `parameters` with its keys unaltered. Function can be used to update a model parameter's value during every repeated simulation. Usefull to propagate MCMC samples of model parameters or perform sensitivity analyses.
* **processes** - (int) - optional - Number of cores to use when {math}`N > 1`.
* **method** - (str) - optional - Integration methods of `scipy.integrate.solve_ivp()` (read the [docs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)).
* **rtol** - (int/float) - optional - Relative tolerance of `scipy.integrate.solve_ivp()` (read the [docs](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)).
Expand Down Expand Up @@ -86,8 +86,8 @@ Upon intialization of the model class, the following arguments must be provided.

* **time** - (int/float or list) - The start and stop "time" for the simulation run. Three possible inputs: 1) int/float, 2) list of int/float of type `[start_time, stop_time]`, 3) list of pd.Timestamp or str of type `[start_date, stop_date]`.
* **N** - (int) - optional - Number of repeated simulations to perform.
* **samples** - (dict) - optional - Dictionary containing samples of model parameters. Obligatory input to a *draw function*.
* **draw_function** - (function) - optional - A function consisting of two inputs: the model parameters dictionary `param_dict`, the previously documented samples dictionary `samples_dict`. Function must return the model parameters dictionary `param_dict`. Function can be used to update a model parameter's value during every repeated simulation. Usefull to propagate MCMC samples of model parameters or perform sensitivity analyses.
* **samples** - (dict) - optional - Dictionary containing samples of the distribution of model parameters. Obligatory input to a *draw function*. Advanced users: dictionary can contain any arbitrary input to a *draw function*.
* **draw_function** - (function) - optional - A function consisting of two inputs: the model parameters dictionary `parameters`, the previously documented samples dictionary `samples`. Function must return the model parameters dictionary `parameters`. Function can be used to update a model parameter's value during every repeated simulation. Usefull to propagate MCMC samples of model parameters or perform sensitivity analyses.
* **processes** - (int) - optional - Number of cores to use when {math}`N > 1`.
* **method** - (str) - optional - 'SSA': Stochastic Simulation Algorithm. 'tau-leap': Tau-leaping algorithm. Consult the following [blog](https://lewiscoleblog.com/gillespie-algorithm) for more background information.
* **tau** - (int/float) - optional - Leap value of the tau-leaping method. Consult the following [blog](https://lewiscoleblog.com/gillespie-algorithm) for more background information.
Expand Down
2 changes: 1 addition & 1 deletion docs/optimization.md
Original file line number Diff line number Diff line change
Expand Up @@ -257,7 +257,7 @@
Samples path, identifier and run_date are combined to find the right .hdf5 `emcee` backend and the `.json` containing the settings.
> **Returns:**
> * **samples_dict** (dict) - Dictionary containing the discarded and thinned MCMC samples and settings.
> * **samples** (dict) - Dictionary containing the discarded and thinned MCMC samples and settings.
## utils.py
Expand Down
45 changes: 35 additions & 10 deletions docs/quickstart.md
Original file line number Diff line number Diff line change
Expand Up @@ -236,20 +236,30 @@ class SIR(JumpProcess):

## Advanced simulation features

### Draw function
### Draw functions

The `sim()` method of `ODE` and `JumpProcess` can be used to perform {math}`N` repeated simulations (keyword `N`). Additionally a *draw function* can be used to update model parameters during every run. A draw function to randomly draw `gamma` from a uniform simulation before every run is implemented as follows,
The `sim()` method of the `ODE` and `JumpProcess` classes can be used to perform {math}`N` repeated simulations by using the optional argument `N`. A *draw function* can be used to alter the value of a model parameters during consecutive model runs by drawing them rrandomly from a list of samples, thereby offering a powerful tool for sensitivity analysis and uncertainty propagation.

Draw functions **always** take two dictionaries as input: 1) The dictionary of model parameters, 2) A dictionary containing samples of a model parameter. Assuming we have a list containing 100 samples of the parameter `gamma`, drawn randomly from 0 to 5,

```python
def draw_function(param_dict, samples_dict):
param_dict['gamma'] = np.random.uniform(low=1, high=5)
return param_dict
# make an example of a dictionary containing parameter samples
samples_dict = {'gamma': np.random.uniform(low=0, high=5, n=100)}

out = model.sim(121, N=10, draw_function=draw_function, samples={})
# define a 'draw function'
def draw_function(parameters, samples):
""" randomly selects a sample of 'gamma' from the provided dictionary of samples and
assigns it to the dictionary of model parameters
"""
parameters['gamma'] = np.random.choice(samples['gamma'])
return parameters

# simulate 10 trajectories, exploit 10 cores to speed up the computation
out = model.sim(121, N=10, draw_function=draw_function, samples=samples_dict, processes=10)
print(out)
```

An additional dimension `'draws'` has been added to the `xarray.Dataset` to accomodate the results of the repeated simulations.
An additional dimension `'draws'` will be added to the `xarray.Dataset` to accomodate the results of the repeated simulations.

```bash
<xarray.Dataset>
Expand All @@ -264,9 +274,24 @@ Data variables:
R (draws, age_groups, time) float64 0.0 0.3439 ... 684.0 684.0
```

### Time-dependent parameter function
This example is more easily coded up by drawing the random values within the *draw function*,

```python
# define a 'draw function'
def draw_function(parameters, samples):
parameters['gamma'] = np.random.uniform(low=0, high=5)
return parameters

# simulate the model
out = model.sim(121, N=10, draw_function=draw_function, samples={})
```
**_NOTE_** Internally, a call to `draw_function` is made within the `sim()` function, where it is given the dictionary of model parameters and the input argument `samples` of the `sim()` function.

**_NOTE:_** Technically, you can supply any input you'd like to a *draw function* by exploiting its `samples` dictionary input. Aside from its name and type, `samples` is not subject to any input checks.

### Time-dependent parameter functions

Parameters can also be varied through the course of one simulation using a *time-dependent parameter function*, which can be an arbitrarily complex function. A generic time-dependent parameter function has the timestep, the dictionary of model states and the value of the varied parameter as its inputs. Additionally, the function can take any number of arguments (including other model parameters).
Parameters can also be varied through the course of one simulation using a *time-dependent parameter function* (TDPF), which can be an arbitrarily complex function. A time-dependent parameter function **always** has the current timestep (`t`), the dictionary containing the model states (`states`) at time `t`, and the original value of the parameter it is applied to (`param`) as its inputs. Additionally, the function can take any number of additional user-defined arguments (including other model parameters).

```python
def vary_my_parameter(t, states, param, an_additional_parameter):
Expand All @@ -278,7 +303,7 @@ def vary_my_parameter(t, states, param, an_additional_parameter):
return param
```

All we need to do is use the `time_dependent_parameters` keyword to declare the parameter our function should be applied to. Additionally, `an_additional_parameter` should be added to the parameters dictionary.
When initialising the model, all we need to do is use the `time_dependent_parameters` keyword to declare what parameter our TDPF should be applied to. Additional parameters introduced in TDPFs, in this example `an_additional_parameter`, should be added to the parameters dictionary.

```python
model = SIR(states={'S': 1000, 'I': 1},
Expand Down
2 changes: 1 addition & 1 deletion docs/start_to_collaborate.md
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@ First of all, thanks for considering contributing to pySODM! It's people like yo

Using pySODM for a research project and writing a paper? Consider citing it:

Tijs W. Alleman, Christian Stevens and Jan. M. Baetens. (2023). pySODM: Simulating and Optimising Dynamical Models in Python 3. arXiv preprint. DOI: https://doi.org/10.48550/arXiv.2301.10664
Tijs W. Alleman, Christian Stevens and Jan. M. Baetens. (2023). pySODM: Simulating and Optimising Dynamical Models in Python 3. Journal of Computational Science. DOI: https://doi.org/10.1016/j.jocs.2023.102148

### Ask a question

Expand Down
Loading

0 comments on commit f7d8ec4

Please sign in to comment.