diff --git a/.github/workflows/draft-pdf.yml b/.github/workflows/draft-pdf.yml new file mode 100644 index 00000000..d6f17740 --- /dev/null +++ b/.github/workflows/draft-pdf.yml @@ -0,0 +1,24 @@ +name: Draft PDF +on: [push] + +jobs: + paper: + runs-on: ubuntu-latest + name: Paper Draft + steps: + - name: Checkout + uses: actions/checkout@v4 + - name: Build draft PDF + uses: openjournals/openjournals-draft-action@master + with: + journal: joss + # This should be the path to the paper within your repo. + paper-path: paper/paper.md + - name: Upload + uses: actions/upload-artifact@v4 + with: + name: paper + # This is the output path where Pandoc will write the compiled + # PDF. Note, this should be the same directory as the input + # paper.md + path: paper/paper.pdf \ No newline at end of file diff --git a/README.md b/README.md index 5182a3e8..2bb9fc98 100644 --- a/README.md +++ b/README.md @@ -88,7 +88,7 @@ The following ODE solvers are available in diffsol 1. A variable order Backwards Difference Formulae (BDF) solver, suitable for stiff problems and singular mass matrices. The basic algorithm is derived in [(Byrne & Hindmarsh, 1975)](#1), however this particular implementation follows that implemented in the Matlab routine ode15s [(Shampine & Reichelt, 1997)](#4) and the SciPy implementation [(Virtanen et al., 2020)](#5), which features the NDF formulas for improved stability 2. A Singly Diagonally Implicit Runge-Kutta (SDIRK or ESDIRK) solver, suitable for moderately stiff problems and singular mass matrices. Two different butcher tableau are provided, TR-BDF2 [(Hosea & Shampine, 1996)](#2) and ESDIRK34 [(Jørgensen et al., 2018)](#3), or users can supply their own. -3. A variable order Explict Runge-Kutta (ERK) solver, suitable for non-stiff problems. One butcher tableau is provided, the 4th order TSIT45 [(Tsitouras, 2011)](#5), or users can supply their own. +3. A variable order Explicit Runge-Kutta (ERK) solver, suitable for non-stiff problems. One butcher tableau is provided, the 4th order TSIT45 [(Tsitouras, 2011)](#5), or users can supply their own. All solvers feature: diff --git a/paper/.gitignore b/paper/.gitignore new file mode 100644 index 00000000..33fe13a5 --- /dev/null +++ b/paper/.gitignore @@ -0,0 +1,2 @@ +*.pdf +jats \ No newline at end of file diff --git a/paper/paper.bib b/paper/paper.bib new file mode 100644 index 00000000..c3d8f95e --- /dev/null +++ b/paper/paper.bib @@ -0,0 +1,182 @@ +@software{cranelift, + author = {{Bytecode Alliance}}, + title = {cranelift}, + url = {https://wasmtime.dev}, + version = {0.115.1}, + date = {2025-03-13}, +} + +@software{nalgebra, + author = {{Dimforge}}, + title = {nalgebra}, + url = {https://github.com/dimforge/nalgebra}, + version = {0.33.2}, + date = {2025-03-13}, +} + +@software{faer, + author = {Sarah El Kazdadi}, + title = {faer-rs}, + url = {https://faer-rs.github.io}, + version = {0.1.0}, + date = {2025-03-13}, +} + +@article{virtanen2020scipy, + title={SciPy 1.0: fundamental algorithms for scientific computing in {P}ython}, + author={Virtanen, Pauli and Gommers, Ralf and Oliphant, Travis E and Haberland, Matt and Reddy, Tyler and Cournapeau, David and Burovski, Evgeni and Peterson, Pearu and Weckesser, Warren and Bright, Jonathan and others}, + journal={Nature methods}, + volume={17}, + number={3}, + pages={261--272}, + year={2020}, + publisher={Nature Publishing Group US New York}, + doi={10.1038/s41592-019-0686-2} +} + +@article{shampine1997matlab, + title={The {MATLAB} {ODE} suite}, + author={Shampine, Lawrence F and Reichelt, Mark W}, + journal={SIAM journal on scientific computing}, + volume={18}, + number={1}, + pages={1--22}, + year={1997}, + publisher={SIAM}, + doi = {10.1137/s1064827594276424} +} + +@article{gardner2022sundials, + title = {Enabling new flexibility in the {SUNDIALS} suite of nonlinear and differential/algebraic equation solvers}, + author = {Gardner, David J and Reynolds, Daniel R and Woodward, Carol S and Balos, Cody J}, + journal = {ACM Transactions on Mathematical Software (TOMS)}, + publisher = {ACM}, + volume = {48}, + number = {3}, + pages = {1--24}, + year = {2022}, + doi = {10.1145/3539801} +} + +@article{DifferentialEquations.jl-2017, + author = {Rackauckas, Christopher and Nie, Qing}, + doi = {10.5334/jors.151}, + journal = {The Journal of Open Research Software}, + keywords = {Applied Mathematics}, + note = {Exported from https://app.dimensions.ai on 2019/05/05}, + number = {1}, + pages = {}, + title = {DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in {J}ulia}, + url = {https://app.dimensions.ai/details/publication/pub.1085583166 and http://openresearchsoftware.metajnl.com/articles/10.5334/jors.151/galley/245/download/}, + volume = {5}, + year = {2017} +} + +@phdthesis{kidger2021on, + title={On Neural Differential Equations}, + author={Patrick Kidger}, + year={2021}, + school={University of Oxford}, + doi={10.48550/arXiv.2202.02435} +} + +@article{byrne1975polyalgorithm, + title={A polyalgorithm for the numerical solution of ordinary differential equations}, + author={Byrne, George D. and Hindmarsh, Alan C.}, + journal={ACM Transactions on Mathematical Software (TOMS)}, + volume={1}, + number={1}, + pages={71--96}, + year={1975}, + publisher={ACM New York, NY, USA}, + doi={10.1145/355626.355636} +} + +@article{bank1985transient, + title={Transient simulation of silicon devices and circuits}, + author={Bank, Randolph E and Coughran, William M and Fichtner, Wolfgang and Grosse, Eric H and Rose, Donald J and Smith, R Kent}, + journal={IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems}, + volume={4}, + number={4}, + pages={436--451}, + year={1985}, + publisher={IEEE}, + doi={10.1109/T-ED.1985.22232} +} + +@article{hosea1996analysis, + title={Analysis and implementation of TR-BDF2}, + author={Hosea, ME and Shampine, LF}, + journal={Applied Numerical Mathematics}, + volume={20}, + number={1-2}, + pages={21--37}, + year={1996}, + publisher={Elsevier}, + doi={10.1016/0168-9274(95)00115-8} +} + +@article{jorgensen2018family, + title={A family of ESDIRK integration methods}, + author={J{\o}rgensen, John Bagterp and Kristensen, Morten Rode and Thomsen, Per Grove}, + journal={arXiv preprint arXiv:1803.01613}, + year={2018}, + doi= {10.48550/arXiv.1803.01613} +} + +@inproceedings{moses22enzyme, +author = {Moses, William S. and Narayanan, Sri Hari Krishna and Paehler, Ludger and Churavy, Valentin and Schanen, Michel and H\"{u}ckelheim, Jan and Doerfert, Johannes and Hovland, Paul}, +title = {Scalable Automatic Differentiation of Multiple Parallel Paradigms through Compiler Augmentation}, +year = {2022}, +isbn = {9784665454445}, +publisher = {IEEE Press}, +booktitle = {Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis}, +articleno = {60}, +numpages = {18}, +keywords = {automatic differentiation, tasks, OpenMP, compiler, Julia, parallel, Enzyme, C++, RAJA, hybrid parallelization, MPI, distributed, LLVM}, +location = {Dallas, Texas}, +series = {SC '22}, +doi = {10.1109/SC41404.2022.00065} +} + +@inproceedings{lattner2004llvm, + title={LLVM: A compilation framework for lifelong program analysis \& transformation}, + author={Lattner, Chris and Adve, Vikram}, + booktitle={International symposium on code generation and optimization, 2004. CGO 2004.}, + pages={75--86}, + year={2004}, + organization={IEEE}, + doi={10.1109/CGO.2004.1281665} +} + +@article{tsitouras2011runge, + title={Runge--{K}utta pairs of order 5 (4) satisfying only the first column simplifying assumption}, + author={Tsitouras, Ch}, + journal={Computers \& mathematics with applications}, + volume={62}, + number={2}, + pages={770--775}, + year={2011}, + publisher={Elsevier}, + doi={10.1016/j.camwa.2011.06.002} +} + +@software{PyO3_Project_and_Contributors, +author = {{PyO3 Project and Contributors}}, +license = {["Apache-2.0", "MIT"]}, +url = {https://github.com/PyO3/pyo3}, +year = {2025}, +version = {v0.27.2}, +title = {{PyO3}} +} + +@article{fritzson2020OpenModelica, +author = {Fritzson, Peter and Pop, Adrian and Abdelhak, Karim and Ashgar, Adeel and Bachmann, Bernhard and Braun, Willi and Bouskela, Daniel and Braun, Robert and Buffoni, Lena and Casella, Francesco and Castro, Rodrigo and Franke, Rüdiger and Fritzson, Dag and Gebremedhin, Mahder and Heuermann, Andreas and Lie, Bernt and Mengist, Alachew and Mikelsons, Lars and Moudgalya, Kannan and Ochel, Lennart and Palanisamy, Arunkumar and Ruge, Vitalij and Schamai, Wladimir and Sjölund, Martin and Thiele, Bernhard and Tinnerholm, John and Östlund, Per}, +doi = {10.4173/mic.2020.4.1}, +journal = {Modeling, Identification and Control}, +number = {4}, +pages = {241--295}, +title = {The {OpenModelica} Integrated Environment for Modeling, Simulation, and Model-Based Development}, +volume = {41}, +year = {2020} +} diff --git a/paper/paper.md b/paper/paper.md new file mode 100644 index 00000000..67018b8f --- /dev/null +++ b/paper/paper.md @@ -0,0 +1,79 @@ +--- +title: '`diffsol`: Rust crate for solving differential equations' +tags: + - rust + - scientific computing + - solver + - ordinary differential equations + - differential algebraic equations + - runge-kutta + - backward differentiation formula +authors: + - name: Martin Robinson + orcid: 0000-0002-1572-6782 + corresponding: true + affiliation: 1 + - name: Alex Allmont + orcid: 0009-0001-4162-0180 + affiliation: 1 + + +affiliations: + - name: Oxford Research Software Engineering Group, Doctoral Training Centre, University of Oxford, Oxford, UK + index: 1 + ror: 052gg0110 +date: 6 Nov 2025 +bibliography: paper.bib +--- + +# Summary + +Ordinary Differential Equations (ODEs) are powerful tools for modelling a wide range of physical systems. Unlike purely data-driven models, ODEs can be based on the underlying physics, biology, or chemistry of the system being modelled, making them particularly useful for predicting the behaviour of a system under conditions that have not been observed. ODEs can be used to model everything from the motion of planets to the spread of infectious diseases. + +[`diffsol`](https://github.com/martinjrobins/diffsol) is a Rust crate for solving ordinary differential equations (ODEs) or semi-explicit differential algebraic equations (DAEs). It can solve equations in the following form: + +$$ +M \frac{dy}{dt} = f(t, y, p) +$$ + +where $y$ is the state of the system, $p$ are a set of parameters, $t$ is time, $f(t, y, p)$ is a function that describes how the state of the system changes over time, and $M$ is an optional and possibly singular mass matrix. The solution to an ODE is a function $y(t)$ that satisfies the ODE and any initial conditions. + +The equations (e.g., $f(t, y, p)$) can be provided by the user either using Rust code or a custom Domain Specific Language (DSL) called `DiffSL`. `DiffSL` uses automatic differentiation using Enzyme [@moses22enzyme] to calculate the necessary gradients, and JIT compilation (using either `LLVM` [@lattner2004llvm] or `Cranelift` [@cranelift]) to generate efficient native code at runtime. The DSL is ideal for using `diffsol` from a higher-level language like Python or R while still maintaining similar performance to pure rust. Diffsol currently provides Python bindings through the [`pydiffsol`](https://github.com/alexallmont/pydiffsol) package, with further language bindings planned. + +ODE solvers require linear algebra containers (e.g., vectors, matrices), operators, and linear solvers. `diffsol` allows users to choose both dense and sparse matrices and solvers from the `nalgebra` [@nalgebra] or `faer` [@faer] crates, and uses a trait-based approach to allow other linear algebra libraries to be added at a later date. + +# Statement of need + +ODE solvers have a long history in scientific computing, and many libraries currently exist. Some notable examples include `scipy.integrate.odeint` [@virtanen2020scipy] in Python, `ode45` [@shampine1997matlab] in MATLAB, and the `Sundials` suite of solvers [@gardner2022sundials] in C. Rust is a systems programming language that is gaining popularity in the scientific computing community due to its performance, safety, and ease of use. There is currently no ODE solver library written in Rust that provides the same level of functionality as these other libraries, and this is the gap that `diffsol` aims to fill. + +ODE solvers written in lower-level languages like C, Fortran, or Rust offer significant performance benefits. However, these solvers are often more difficult to wrap and use in higher-level languages like Python or MATLAB, primarily because users must supply their equations in the language of the solver. `diffsol` solves this issue by providing its own custom `DiffSL` DSL which is JIT compiled to efficient native code at run-time, meaning that users from a higher-level language like Python or R can specify their equations using a simple string-based format while still maintaining performance similar to that of pure Rust. Two other popular ODE solvers that take advantage of JIT compilation are `DifferentialEquations.jl` [@DifferentialEquations.jl-2017] in Julia and `diffrax` [@kidger2021on] in Python. However, both these packages compile the entire solver as well as the equations, which is a significant amount of code. `diffSol` only compiles the equations, meaning that it has a significantly smaller "time-to-first-plot" for users. Another popular differential equations solver package utilising a DSL is OpenModelica [@fritzson2020OpenModelica]. Wrappers to this package in higher-level languages like Python rely on messaging to a separate OpenModelica server, which can be slow and more complicated to set up. In contrast, `diffsol` can be integrated directly into higher-level languages using language bindings and linking to a single shared library, see for example, the `pydiffsol` Python bindings discussed below. + +# Features + +The following solvers are available in `diffsol`: + +1. A variable order Backwards Difference Formulae (BDF) solver, suitable for stiff problems and singular mass matrices. The basic algorithm is derived in [@byrne1975polyalgorithm], however this particular implementation follows that implemented in the MATLAB routine ode15s [@shampine1997matlab] and the SciPy implementation [@virtanen2020scipy], which features the NDF formulas for improved stability. +2. A Singly Diagonally Implicit Runge-Kutta (SDIRK or ESDIRK) solver, suitable for moderately stiff problems and singular mass matrices. Two different butcher tableau are provided, TR-BDF2 [@bank1985transient; @hosea1996analysis] and ESDIRK34 [@jorgensen2018family], or users can supply their own. +3. A variable order Explicit Runge-Kutta (ERK) solver, suitable for non-stiff problems. One butcher tableau is provided, the 4th order TSIT45 [@tsitouras2011runge], or users can supply their own. + +All solvers feature: + +- Linear algebra containers and linear solvers from the `nalgebra` or `faer` crates, including both dense and sparse matrix support. +- Adaptive step-size control to given relative and absolute tolerances. Tolerances can be set separately for the main equations, quadrature of the output function, and sensitivity analysis. +- Dense output, interpolating to times provided by the user. +- Event handling, stopping when a given condition $g_e(t, y, p)$ is met or at a specific time. +- Numerical quadrature of an optional output $g_o(t, y, p)$ function over time. +- Forward sensitivity analysis, calculating the gradient of an output function or the solver states $y$ with respect to the parameters $p$. +- Adjoint sensitivity analysis, calculating the gradient of cost function $G(p)$ with respect to the parameters $p$. The cost function can be the integral of a continuous output function $g(t, y, p)$ or a sum of a set of discrete functions $g_i(t_i, y, p)$ at time points $t_i$. + +# Bindings in higher-level languages + +[`pydiffsol`](https://github.com/alexallmont/pydiffsol) provides Python bindings to `diffsol` using the `PyO3` [@PyO3_Project_and_Contributors] crate. It allows users to define ODEs in Python using the `DiffSL` DSL, and solve them using the Rust `diffsol` library. `pydiffsol` aims to provide a simple and easy-to-use interface for solving ODEs in Python, while still maintaining the performance benifits of using Rust under the hood. + +The goal is to develop further [bindings to other higher-level languages](https://github.com/martinjrobins/diffsol/issues/131), including R, JAX, and JavaScript, exploiting the `DiffSL` DSL to provide high performance while maintaining ease of use and positioning the core `diffsol` library as a widely used cross-language, cross-platform, high-performance ODE solver. + +# Acknowledgements + +We gratefully acknowledge the support of all the past and future contributors to the `diffsol` project, for their advice, enthusiasm, bug reports and code. In particular, we would like to thank the authors of the `pharmsol` crate, Julian Otalvaro and Markus Hovd. + +# References