atropy_core: DLR approximation for the kinetic CME

• atropy_core: DLR approximation for the kinetic CME
  • Objective
  • Requirements
  • Optional requirements
  • Installation
    • Intel MKL
    • OpenMP
    • pybind11
    • Python environment
  • Run the program
  • Input
    • Preparing input data
  • Output
  • Example problems
  • References

Objective

atropy_core solves the chemical master equation (CME), $$\partial_t{P}(t,\mathbf{x}) = \sum_{\mu = 0}^{M-1}\left(\alpha_\mu(\mathbf{x}-\bm{\nu}\mu)P(t,\mathbf{x}-\bm{\nu}\mu) - \alpha_\mu(\mathbf{x})P(t,\mathbf{x})\right)$$

according to the algorithm proposed in [1], which is based on the projector-splitting integrator for tree tensor networks [2].

$P(t,\mathbf{x}),\mathrm{d}t$ is the probability of finding a population number of $\mathbf{x} = (x_0, \dots, x_{N-1})$ molecules of species $S_0, \dots, S_{N-1}$ at time $t$. The CME describes the time evolution of this probability distribution $P(t,\mathbf{x})$ in a chemical reaction network with $N$ different species $S_0, \dots, S_{N-1}$, which can react via $M$ reaction channels $R_0, \dots, R_{M-1}$. For a given reaction $\mu$, the stoichiometric vector $\bm{\nu}\mu$ denotes the population change by that reaction and the propensity functions $\alpha\mu(\mathbf{x})$ and $\alpha_\mu(\mathbf{x}-\bm{\nu}\mu)$ are proportional to the transition probabilities $T(\mathbf{x}+\bm{\nu}\mu|\mathbf{x})$ and $T(\mathbf{x}|\mathbf{x}-\bm{\nu}_\mu)$.

atropy_core makes use of the low-rank framework Ensign [3].

Requirements

• C++20 compatible C++ compiler
• CMake (3.27 or later)
• Eigen 3.4
• netCDF4 (built together with HDF5)
• Python (>3.10)

Check via nc-config --has-hdf5, whether HDF5 was used in the netCDF4 build.

Optional requirements

• Fortran compiler (if OpenBLAS is used)
• Intel MKL
• OpenMP

See below for a description how Intel MKL and OpenMP can be enabled.

Installation

Build the program in the <build> directory by executing

cmake -B <build> -DCMAKE_BUILD_TYPE=Release
cmake --build <build>

in the project root. The generated executable atropy_core can be found in <build>.

To enable compiler options for debugging, use -DCMAKE_BUILD_TYPE=Debug instead. Unit tests for C++ files are provided in the tests folder. They can be built by activating the -DBUILD_TESTING flag. Run the tests with

ctest --test-dir <build>

Intel MKL

If you prefer to use Intel MKL as the BLAS and LAPACK backend instead of OpenBLAS set the -DMKL_INCLUDEDIR and -DMKL_LIBDIR variables and build with

cmake -B <build> -DCMAKE_BUILD_TYPE=Release -DMKL_ENABLED=ON
cmake --build <build>

Make sure to add the MKL libraries to your LD_LIBRARY_PATH, i.e.

export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/path/to/intel/mkl/lib/intel64_lin/

before running the executable.

OpenMP

OpenMP can be activated via

cmake -B <build> -DCMAKE_BUILD_TYPE=Release -DOPENMP=ON

Make sure that the OMP_NUM_THREADS environment variable is in accordance with your hardware specification and run the unit tests via

ctest --test-dir <build>

to ensure that OpenMP and atropy_core work correctly.

macOS: Note that XCode compilers do not support OpenMP. For using OpenMP on macOS, a manual installation (e.g. of gcc11) is required and the CXX, CC and FC environment variables have to be set accordingly.

pybind11

The tree tensor network integrator IntegrateTTN of atropy_core,

IntegrateTTN(std::string input, std::string output, int snapshot, double tau, double tfinal, unsigned int substeps, char method)

where parameters are the same as the command line arguments for atropy_core (cf. next section), is available in the pybind11 module atropy_core_pybind11 and thus can be used in Python. The Python wrapper for atropy_core, called atropy, makes use of this module. Build atropy_core with the -DPYBIND11_ENABLED flag to generate the pybind11 module.

macOS: pybind11 might not find the Python.h header during the CMake build process. In that case, set export CPLUS_INCLUDE_PATH=<path/to/python/include> accordingly.

Python environment

To use the Python programs included in atropy_core, a Python environment with external packages specified in pyproject.toml needs to be configured and enabled. This can be done using uv pip install or classical pip install.

For uv:

uv venv my_venv --python <python_version>
source my_venv/bin/activate
uv pip install .

For Python venv:

python -m venv path/to/my_venv
source path/to/my_venv/bin/activate
pip install .

For anaconda:

conda create -n my_venv python
conda activate my_venv
pip install .

All scripts have to be executed from the project root. When using a IDE, make sure to adjust the settings accordingly. Unit tests for Python files are located in the atropy_core/tests folder. They can be run in the Python environment via

pytest atropy_core/tests

Run the program

atropy_core has to be run with

  ./build/atropy_core [OPTION...]

and expects the following command line arguments:

• -i, --input: Name of the input .nc file (default: input/input.nc)
• -o, --output: Name of the output folder, stored in output/
• -s, --snapshot: Number of steps between two snapshots
• -t, --tau: Time step size
• -f, --tfinal: Final integration time
• -n, --substeps: Number of integration substeps (default: 1)
• -m, --method: Integration method (e (explicit Euler), r (explicit RK4), i (implicit Euler), c (Crank-Nicolson)) (default: i)
• -h, --help: Print usage

Input

Input netCDF files have to be stored as input/input.nc (the directory can be changed using the -i flag) and can be generated with the input scripts provided in atropy_core/examples.

Caution: As Ensign stores arrays in column-major (Fortran) order, it is assumed that input arrays also follow this convention.

Preparing input data

Let us consider the input script set_lambda_phage.py located in the atropy_core/examples folder, which generates input data for the lambda phage model. It gives an example on how the initial conditions have to be set up. A short documentation for this script is provided by

python3 atropy_core/examples/set_lambda_phage.py --help

Note that atropy_core assumes that the propensity function is factorizable for the species in different partitions. However, the input scripts rely on the ReactionSystem class (cf. atropy_core/reaction.py), which assumes that the propensity function is factorizable in all species. This is a valid assumption for most scenarios. For problems where species in a partition are not factorizable, the propensity function can be adjusted manually after initializing the Tree with the method initialize.

Output

atropy_core automatically creates a folder in output/ with a name set by the -o/--output parameter. The low-rank factors and coupling coefficients as well as the chosen model parameters are stored in this folder as output_t<ts>.nc (<ts> denotes the time step) in intervals according to the -s/--snapshot parameter.

Example problems

Input generation scripts for the example problems (toggle switch, lambda phage and reaction cascade) are provided in atropy_core/examples and the corresponding model files can be found in atropy_core/examples/models.

The Python file output_helper.py contains helper functions such as readTree and the TimeSeries class for plotting the NetCDF files.

References

[1] Einkemmer, L., Mangott, J., and Prugger, M.: "A hierarchical dynamical low-rank algorithm for the stochastic description of large reaction networks", arXiv, https://arxiv.org/abs/2407.11792 (2024)

[2] Ceruti, G., Lubich, C., and Walach, H.: "Time integration of Tree Tensor networks", SIAM J. Numer. Anal. 59 (2021)

[3] Cassini, F., and Einkemmer, L.: "Efficient 6D Vlasov simulation using the dynamical low-rank framework Ensign", Comp. Phys. Comm. 280 (2022)