Using Julia via pyjulia

This notebook demonstrates the application of pyABC to models defined in Julia.

pyABC’s Julia interface is the class pyabc.external.julia.Julia.

[2]:
import julia

julia.install()
# For further information, see https://pyjulia.readthedocs.io/en/latest/installation.html.
# There are some known problems, e.g. with statically linked Python interpreters, see
# https://pyjulia.readthedocs.io/en/latest/troubleshooting.html for details.
[ Info: Julia version info
Julia Version 1.12.2
Commit ca9b6662be4 (2025-11-20 16:25 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  uname: Darwin 25.3.0 Darwin Kernel Version 25.3.0: Wed Jan 28 20:53:15 PST 2026; root:xnu-12377.81.4~5/RELEASE_ARM64_T6000 arm64 arm
  CPU: Apple M1 Pro:
                 speed         user         nice          sys         idle          irq
       #1-10  2400 MHz     207298 s          0 s     122616 s    3611702 s          0 s
  Memory: 16.0 GB (181.53125 MB free)
  Uptime: 42240.0 sec
  Load Avg:  5.88818359375  4.306640625  3.7685546875
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, apple-m1)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  HOMEBREW_PREFIX = /opt/homebrew
  INFOPATH = /opt/homebrew/share/info:
  PYTHONPATH = /Users/jonas.arruda/PyCharm Projects/pyABC
  HOME = /Users/jonas.arruda
  HOMEBREW_REPOSITORY = /opt/homebrew
  PATH = /Users/jonas.arruda/PyCharm Projects/pyABC/.venv/bin:/usr/local/bin:/System/Cryptexes/App/usr/bin:/usr/bin:/bin:/usr/sbin:/sbin:/var/run/com.apple.security.cryptexd/codex.system/bootstrap/usr/local/bin:/var/run/com.apple.security.cryptexd/codex.system/bootstrap/usr/bin:/var/run/com.apple.security.cryptexd/codex.system/bootstrap/usr/appleinternal/bin:/opt/pmk/env/global/bin:/Library/TeX/texbin:/Users/jonas.arruda/.juliaup/bin:/Users/jonas.arruda/.codeium/windsurf/bin:/Users/jonas.arruda/miniconda3/bin:/Library/Frameworks/Python.framework/Versions/3.10/bin:/Library/Frameworks/Python.framework/Versions/3.11/bin:/Library/Frameworks/Python.framework/Versions/3.12/bin:/opt/homebrew/bin:/opt/homebrew/sbin:/Users/jonas.arruda/.local/bin
  XPC_FLAGS = 0x0
  HOMEBREW_CELLAR = /opt/homebrew/Cellar
  TERM = xterm-color
[ Info: Julia executable: /Users/jonas.arruda/.julia/juliaup/julia-1.12.2+0.aarch64.apple.darwin14/bin/julia
[ Info: Trying to import PyCall...
Info: PyCall is already installed and compatible with Python executable.

PyCall:
    python: /Users/jonas.arruda/PyCharm Projects/pyABC/.venv/bin/python
    libpython: /Users/jonas.arruda/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/libpython3.12.dylib
Python:
    python: /Users/jonas.arruda/PyCharm Projects/pyABC/.venv/bin/python
    libpython: /Users/jonas.arruda/.local/share/uv/python/cpython-3.12.12-macos-aarch64-none/lib/libpython3.12.dylib
[6]:
from julia import Pkg

Pkg.add('Catalyst')
Pkg.add('JumpProcesses')
   Resolving package versions...
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
   Resolving package versions...
     Project No packages added to or removed from `~/.julia/environments/v1.12/Project.toml`
    Manifest No packages added to or removed from `~/.julia/environments/v1.12/Manifest.toml`
[ ]:
import tempfile

import matplotlib.pyplot as plt

import pyabc
from pyabc import ABCSMC, RV, Distribution, MulticoreEvalParallelSampler
from pyabc.external.julia import Julia

pyabc.settings.set_figure_params('pyabc')  # for beautified plots

As demonstration example, we use an SIR disease dynamics model. For simulation, we use an implementation of Gillespie’s direct algorithm in the DifferentialEquations.jl package.

The code consists of multiple functions in the file model_julia/SIR.jl, wrapped in the namespace of a module SIR. Importing the module and dependencies can take some time due to pre-processing.

[ ]:
%%time
jl = Julia(module_name='SIR', source_file='model_julia/SIR.jl')
[ ]:
jl.display_source_ipython()

The Julia code defines functions for model and distance. The model returns a dictionary, whose entries are internally converted to numpy arrays, while the distance returns a single floating value. Further, also the observed data can be defined in and imported from Julia. Note that also only a subset of model, distance and observation can be defined in Julia, and combined with Python objects.

[ ]:
model = jl.model()
distance = jl.distance()
obs = jl.observation()

_ = plt.plot(obs['t'], obs['u'])

We use a prior on log-scale, the Julia model applies the corresponding transformation internally.

[ ]:
gt_par = {'p1': -4.0, 'p2': -2.0}

# parameter limits and prior
par_limits = {
    'p1': (-5, -3),
    'p2': (-3, -1),
}
prior = Distribution(
    **{key: RV('uniform', lb, ub - lb) for key, (lb, ub) in par_limits.items()}
)

If we use parallelization, we need to call all julia functions once before, to perform necessary pre-compiling:

[ ]:
distance(model(gt_par), model(gt_par))

We are all set to run an analysis, even with parallelization via multi-processing:

[ ]:
abc = ABCSMC(
    model,
    prior,
    distance,
    sampler=MulticoreEvalParallelSampler(),
)
db = tempfile.mkstemp(suffix='.db')[1]
abc.new('sqlite:///' + db, obs)
h = abc.run(max_nr_populations=10)

The posterior marginals show that the data constrain the parameters well:

[ ]:
for t in [0, h.max_t]:
    pyabc.visualization.plot_kde_matrix_highlevel(
        h,
        t=t,
        limits=par_limits,
        refval=gt_par,
        refval_color='grey',
    )
    plt.gcf().suptitle(f'Posterior at t={t}')
    plt.gcf().tight_layout()

The accepted simulations closely match the observed data well, as we can see from data plots:

[ ]:
def plot_data(sumstat, weight, ax, **kwargs):  # noqa: ARG001
    """Plot a single trajectory"""
    for i in range(3):
        ax.plot(sumstat['t'], sumstat['u'][:, i], color=f'C{i}', alpha=0.1)


for t in [0, h.max_t]:
    _, ax = plt.subplots()
    pyabc.visualization.plot_data_callback(
        h,
        plot_data,
        t=t,
        ax=ax,
    )
    ax.plot(obs['t'], obs['u'])
    ax.set_title(f'Simulations at t={t}')