Package initialization

Scientific computing with Julia

David Widmann (@devmotion )

Uppsala University, Sweden

American University in Bulgaria, 17 February 2022

About me

• 👨‍🎓 PhD student at the Department of Information Technology and the Centre for Interdisciplinary Mathematics (CIM) in Uppsala

  • BSc and MSc in Mathematics (TU Munich)

  • Human medicine (LMU and TU Munich)

• 👨‍🔬 Research topic: "Uncertainty-aware deep learning"

  • I.e., statistics, probability theory, machine learning, and computer science

• 💻 Julia programming, e.g.,

  • SciML, in particular DelayDiffEq.jl

  • Turing ecosystem

Outline

We will see examples from and I will discuss my relation to

• SciML (a lot)

• Turing (a bit)

• my research (a bit)

This is a very opinionated choice, there are many other great packages and organizations in Julia for scientific computing.

Please interrupt me at any time if you have a question 🙂

"Why we created Julia"

Bezanson, Karpinski, Shah, and Edelman (2012):

In short, because we are greedy.

We want a language that's open source, with a liberal license.

We want the speed of C with the dynamism of Ruby.

We want something as usable for general programming as Python, as easy for statistics as R, as natural for string processing as Perl, as powerful for linear algebra as Matlab, as good at gluing programs together as the shell.

Something that is dirt simple to learn, yet keeps the most serious hackers happy.

We want it interactive and we want it compiled.

Why we use Julia, 10 years later

14 February 2022

Why do I use Julia?

I was curious and excited early on:

• a, supposedly fast and cool, programming language! 🎉

• that is open source 📖

• and has a nice math-like syntax 😍

At TUM I had to learn and use MATLAB and R...

...but I used the opportunity and tried to work and become familiar with it for a course project during my Erasmus exchange in Uppsala.

Initial contributions

I had noticed some things that were missing and maybe could be improved.

After I was done with my exchange year, I had some time and started to make my first small contributions.

Barabási-Albert random graphs

First PR (archived, use Graphs.jl nowadays)

Weighted sampling without replacement

Faster algorithms

SciML

Excerpt from the NumFOCUS project page:

SciML is an open source software organization created to unify the packages for scientific machine learning. This includes the development of modular scientific simulation support software, such as differential equation solvers, along with the methodologies for inverse problems and automated model discovery. By providing a diverse set of tools with a common interface, we provide a modular, easily-extendable, and highly performant ecosystem for handling a wide variety of scientific simulations.

SciML tools are used by many organizations, such as (but not limited to!) the CliMA: Climate Modeling Alliance, the New York Federal Reserve Bank, the Julia Robotics community, Pumas-AI: Pharmaceutical Modeling and Simulation, and the Brazilian National Institute for Space Research (INPE).

Disclaimer: I am a Github admin and member of the steering council of SciML, so I'm definitely biased 🙃

Ordinary differential equation (ODE): Lotka-Volterra model

Differential equation model of predator-prey interaction:

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}🐰& =\alpha 🐰-\beta 🐰🦊\\ \frac{\mathrm{d}}{\mathrm{d}t}🦊& =-\gamma 🦊+\delta 🐰🦊\end{array}$$

Problem formulation with initial value, integration time span, and default parameters:

Solve problem:

Different parameter choices:

Why don't we always use Euler's method?

Reference: Blog by Chris Rackauckas

A chaotic model of random neural networks (reference solution, computed with a high-order algorithm and very low tolerance):

A corresponding work-precision diagram:

Whenever someone in 2019 uses more code to internally implement a quick and dirty solver instead of running standard codes, a baby sheds a tear. It's not just about error or speed (or even accuracy: libraries undergo lots of tests!), it's about how much error you get with a certain speed. Euler's method (and even RK4) just don't scale well if you need non-trivial error on most equations.

Comparison of different solvers

Lotka-Volterra (non-stiff)

ROBER (stiff)

Only ODEs?!

Supported equations include

• Discrete equations (function maps, discrete stochastic (Gillespie/Markov) simulations)

• Ordinary differential equations (ODEs)

• Split and Partitioned ODEs (Symplectic integrators, IMEX Methods)

• Stochastic ordinary differential equations (SODEs or SDEs)

• Stochastic differential-algebraic equations (SDAEs)

• Random differential equations (RODEs or RDEs)

• Differential algebraic equations (DAEs)

• Delay differential equations (DDEs)

• Neutral, retarded, and algebraic delay differential equations (NDDEs, RDDEs, and DDAEs)

• Stochastic delay differential equations (SDDEs)

• Experimental support for stochastic neutral, retarded, and algebraic delay differential equations (SNDDEs, SRDDEs, and SDDAEs)

• Mixed discrete and continuous equations (Hybrid Equations, Jump Diffusions)

• (Stochastic) partial differential equations ((S)PDEs) (with both finite difference and finite element methods)

My relation to SciML

I wrote my master thesis about a mathematical model of quorum sensing of P. putida.

doi:10.3390/app6050149

Here lactonase activity is regulated by PpuR-AHL complex with a specific time delay:

$$\begin{array}{rl}{S}^{\prime }(t)=& D(S_0-S(t))-\frac{{\gamma }_S{S(t)}^{n_S}}{K_m^{n_S}+{S(t)}^{n_S}}N(t)\\ N^{\prime }(t)=& N(t)(\frac{a{S(t)}^{n_S}}{K_m^{n_S}+{S(t)}^{n_S}}-D)\\ A^{\prime }(t)=& \begin{array}{r}({\alpha }_A+\frac{{\beta }_A{C(t)}^{n_1}}{C_1^{n_1}+{C(t)}^{n_1}})N(t)-{\gamma }_{A}A(t)-DA(t)\\ +{\gamma }_{3}C(t)-{\alpha }_C(R_{tot}-C(t))A(t)-K_{E}A(t)L(t)\end{array}\\ C^{\prime }(t)=& {\alpha }_C(R_{tot}-C(t))A(t)-{\gamma }_{3}C(t)\\ L^{\prime }(t)=& \frac{{\alpha }_L{C(t-\tau )}^{n_2}}{C_2^{n_2}+{C(t-\tau )}^{n_2}}N(t)-{\gamma }_{L}L(t)-DL(t)\end{array}$$

Experimental data and numerical simulations:

doi:10.3390/app6050149 (adapted from doi:10.1007/s00216-014-8063-6)

Simulation with MATLAB (dde23)

Simulation with DelayDiffEq.jl

Similar issues with default settings:

Decreasing the step size improves the solution (also in MATLAB):

But this increases the number of steps and hence makes the computation slower! 😢

I really wanted to perform computations both efficiently and accurate.

So I decided to

• work with Julia and improve DelayDiffEq

• add support for more and also stiff methods based on the vast amount of existing methods for ODEs,

• add callbacks that preserve domain constraints

Stiff methods

Methods such as MethodOfSteps(Rosenbrock23()) can solve the DDE without restricting the step size:

A simpler example: Mackey-Glass equation

Model of circulating blood cells $x(t)$ at time $t$ by Mackey and Glass

Given by

$$\begin{array}{rl}{x}^{\prime }(t)& =\frac{\beta x(t-\tau )}{1+{x(t-\tau )}^n}-\gamma x(t),t\ge 0,\\ x(t)& =x_0(t),t\in [-\tau ,0]\end{array}$$

• Blood cells are destroyed at rate $\gamma$

• Their production rate depends on concentration $x(t-\tau )$ since

  [t]here is a significant delay $\tau$ between the initiation of cellular production in the bone marrow and the release of mature cells into the blood

Define DDE problem with $x_0=0.5$, $\beta =2$, $\gamma =1$, $n=9.65$, and $\tau =2$ (values taken from here):

Only differential equations?!

Core components are

• Differential Equation Solving

• Physics-Informed Model Discovery and Learning

• Bridges to Python and R

• Compiler-Assisted Model Analysis and Sparsity Acceleration

• ML-Assisted Tooling for Model Acceleration

• Differentiable Scientific Data Structures and Simulators

• Tools for Accelerated Algorithm Development and Research

Parameter estimation

Some noisy data:

SciML supports many different packages for local and global optimization such as Optim.jl, LeastSquaresOptim.jl, JuMP, and any package that implements the MathProgBase interface such as NLOpt.jl.

The default evoluationary optimization algorithm in the package BlackBoxOptim.jl yields the following optimal parameters:

Optimized trajectories:

Neural ODEs

Example 1

Neural ODE:

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}🐰\\ 🦊\end{array}\right]& =\mathrm{N}\mathrm{e}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{k}\left(\left[\begin{array}{c}🐰\\ 🦊\end{array}\right]\right)\end{array}$$

Example 2

Physics-informed neural ODE model (or rather biology-informed in this case):

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}🐰\\ 🦊\end{array}\right]& =\left[\begin{array}{c}🐰\\ 🦊\end{array}\right]\odot (\left[\begin{array}{c}\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}(c_1)\\ -\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}(c_2)\end{array}\right]+\mathrm{N}\mathrm{e}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{l}\mathrm{N}\mathrm{e}\mathrm{t}\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{k}\left(\left[\begin{array}{c}🐰\\ 🦊\end{array}\right]\right))\end{array}$$

Probabilistic modelling

We only know noisy (i.e., slightly incorrect) number of predator and prey for some days but not parameters.

Idea: model uncertainty of unknown parameters with probability distributions

Bayesian approach

Model uncertainty with conditional probability

$$p_{\mathrm{\Theta }|Y}(\theta |y)$$

of unknown parameters $\theta$ given observations $y$.

E.g.,

• $\theta$: unknown parameters $\alpha$, $\beta$, $\gamma$, $\delta$, and $u(0)$ that we want to infer

• $y$: number of observed predator and prey

Bayes' theorem

Conditional probability $p_{\mathrm{\Theta }|Y}(\theta |y)$ can be calculated as

$$p_{\mathrm{\Theta }|Y}(\theta |y)=\frac{p_{Y|\mathrm{\Theta }}(y|\theta )p_{\mathrm{\Theta }}(\theta )}{p_Y(y)}$$

Workflow

• Choose model $p_{Y|\mathrm{\Theta }}(y|\theta )$

  • e.g., describes how daily cases of newly infected individuals depend on parameters

• Choose prior $p_{\mathrm{\Theta }}(\theta )$

  • should incorporate initial beliefs and knowledge about $\theta$

  • e.g., $\alpha$, $\beta$, $\gamma$, $\delta$, and $u$ are non-negative

• Compute

  $$p_{\mathrm{\Theta }|Y}(\theta |y)=\frac{p_{Y|\mathrm{\Theta }}(y|\theta )p_{\mathrm{\Theta }}(\theta )}{p_Y(y)}=\frac{p_{Y|\mathrm{\Theta }}(y|\theta )p_{\mathrm{\Theta }}(\theta )}{\int p_{Y|\mathrm{\Theta }}(y|{\theta }^{\prime })p_{\mathrm{\Theta }}({\theta }^{\prime })\mathrm{d}{\theta }^{\prime }}$$

Discrete approximation

Idea: Approximate $p_{\mathrm{\Theta }|Y}(\theta |y)$ with a weighted mixture of point measures

$$p_{\mathrm{\Theta }|Y}(\cdot |y)\approx \sum _i{w}_i{\delta }_{{\theta }_i}(\cdot )$$

where $w_i>0$ and $\sum _i{w}_i=1$.

This implies that

$$\mathbb{E}(\mathrm{\Theta }|Y=y)\approx \sum _i{w}_i{\theta }_i$$

and more generally

$$\mathbb{E}(\varphi (\mathrm{\Theta })|Y=y)\approx \sum _i{w}_i\varphi ({\theta }_i)$$

number of samples:

500

Probabilistic programming

Source: Lecture slides

Developing probabilistic models and inference algorithms is a time-consuming and error-prone process.

• Probabilistic model written as a computer program

• Automatic inference (integral part of the programming language)

Advantages:

• Fast development of models

• Expressive models

• Widely applicable inference algorithms

Turing

Turing is a probabilistic programming language (PPL).

General design

• Probabilistic models are implemented as a Julia function

• One may use any Julia code inside of the model

• Random variables and observations are declared with the ~ operator:

  @model function mymodel(x, y)
      ...
      # random variable `a` with prior distribution `dist_a`
      a ~ dist_a
  
      ...
  
      # observation `y` with data distribution `dist_y`
      y ~ dist_y
      ...
  end

• PPL is implemented in DynamicPPL.jl, including @model macro

• Turing.jl integrates and reexports different parts of the ecosystem such as the PPL, inference algorithms, and tools for automatic differentiation