SODECL: An Open Source Library for Calculating Multiple Orbits of a System of Stochastic Differential Equations in Parallel

TL;DR

SODECL is an open-source library that leverages multi-core CPUs and GPUs to efficiently compute multiple orbits of stochastic differential equations, significantly reducing computation time for complex models.

Contribution

The paper introduces SODECL, a novel software library that utilizes parallel computing devices to accelerate the simulation of stochastic differential equations.

Findings

Up to 6.7x acceleration using 32 CPU cores

Up to 4.5x acceleration on high-end GPU

Scalability depends on model size and hardware

Abstract

Stochastic differential equations (SDEs) are widely used to model systems affected by random processes. In general, the analysis of an SDE model requires numerical solutions to be generated many times over multiple parameter combinations. However, this process often requires considerable computational resources to be practicable. Due to the embarrassingly parallel nature of the task, devices such as multi-core processors and graphics processing units (GPUs) can be employed for acceleration. Here, we present {\bf SODECL} (\url{https://github.com/avramidis/sodecl}), a software library that utilises such devices to calculate multiple orbits of an SDE model. To evaluate the acceleration provided by SODECL, we compared the time required to calculate multiple orbits of an exemplar stochastic model when one CPU core is used, to the time required when using all CPU cores or a GPU. In addition, to assess scalability, we investigated how the model size affected execution time on different parallel compute devices. Our results show that when using all 32 CPU cores of a high-end high-performance computing node, the task is accelerated by a factor of up to $\simeq$6.7, compared to when using a single CPU core. Executing the task on a high-end GPU yielded accelerations of up to $\simeq$4.5, compared to a single CPU core.

Tables

Table 1. Table 1: List of named constants used by the SODECL OpenCL functions.

Named constant	Description
_numeq_	Number of equations in the SDE system
_numnoi_	Number of noise variables in the SDE system
_numpar_	Number of parameters in the SDE system
_m_dt_	SDE integrator time step in seconds
_numsteps_	Number of SDE integrator steps per OpenCL kernel call

Table 2. Table 2: Arguments of the SODECL Python interface.

Argument	Variable	Description
1	platform	OpenCL platform number
2	device	OpenCL device number of the selected platform
3	kernel	Path of the file with the OpenCL function defining the SDE system
4	initx	Initial conditions for each orbit of the ODE system
5	params	Parameter sets for the SDE system for all orbits
6	integrator	SODE integrator
7	orbits	Number of orbits to be calculated
8	nequat	Number of equations of the SDE system
9	nparams	Number of parameters of the SDE system
10	nnoi	Number of noise processes
11	dt	SDE solver time step
12	tspan	Integration time span
13	ksteps	Number of SDE integrator steps executed in each OpenCL kernel call
14	localgroupsize	OpenCL local group size

Table 3. Table 3: Hardware used for the speed evaluation tests.

ID	Type	Model
I4790K	CPU	Intel Core i7-4790K
2X6142	CPU	2 x Intel Xeon Gold 6142
W8100	GPU	AMD Firepro W8100
P100	GPU	NVIDIA Tesla P100-PCIE

Table 4. Table 4: Local group sizes used for each compute device to generate the runtime results shown in Tables  S1 - S6 and plotted in Figure 7 .

Model size, $N$	I4790K	2X6142	W8100	P100
5	8	32	256	8
10	8	32	16	8
15	8	32	32	8

Table 5. Table S1: Runtimes in seconds (s) for the stochastic Kuramoto model as a function of the number of orbits for system size N = 5 𝑁 5 N=5 , when using different parallel computer hardware implementations. The hardware used in each case is shown in Table  3 .

Orbits	2X6142 (1C)	2X6142 (32C)	I4790K (4C)	W8100	P100
512	0.5044	0.449	0.728	0.552	0.201
5120	0.7463	0.401	0.898	0.617	0.363
25600	2.3734	0.814	1.612	1.424	0.935
40960	3.6019	1.163	2.254	1.973	1.400
81920	6.8686	1.961	3.800	3.573	2.512
163840	13.390	3.623	6.779	6.772	4.719

Table 6. Table S2: Runtimes in second (s) for the stochastic Kuramoto model as a function of the number of orbits for system size N = 10 𝑁 10 N=10 , when using different parallel computer hardware implementations. The hardware used in each case is shown in Table  3 .

Orbits	2X6142 (1C)	2X6142 (32C)	I4790K (4C)	W8100	P100
512	0.441	0.301	0.714	0.708	0.280
5120	1.444	0.499	1.138	0.864	0.513
25600	5.956	1.385	3.160	2.784	1.714
40960	9.676	2.042	4.746	4.211	2.566
81920	18.308	3.821	8.784	7.981	4.828
163840	36.641	7.471	16.939	15.589	9.315

Table 7. Table S3: Runtimes in second (s) for the stochastic Kuramoto model as a function of the number of orbits for system size N = 15 𝑁 15 N=15 , when using different parallel computer hardware implementations. The hardware used in each case is shown in Table  3 .

Orbits	2X6142 (1C)	2X6142 (32C)	I4790K (4C)	W8100	P100
512	0.8247	0.495	0.772	1.271	0.323
5120	2.8846	0.759	1.605	1.679	0.672
25600	12.6576	2.223	5.256	5.968	2.940
40960	19.9707	3.288	7.987	9.128	4.566
81920	39.6412	6.185	15.279	17.736	9.056
163840	79.5124	11.953	29.882	34.700	17.656

Table 8. Table S4: Speedups obtained for the stochastic Kuramoto model as a function of the number of orbits for system size N = 5 𝑁 5 N=5 , when using each of the compute devices listed in Table  3 . Speedups are calculated as the runtime for one core of the Intel Xeon 6142 CPU divided by the runtime of the compute device. In the case of the CPUs (2X6142 and 4790K), all cores were used.

Orbits	2X6142	I4790K	W8100	P100
512	1.122	0.692	0.914	2.498
5120	1.861	0.831	1.209	2.051
25600	2.915	1.471	1.666	2.537
40960	3.095	1.597	1.826	2.572
81920	3.501	1.807	1.922	2.733
163840	3.696	1.975	1.977	2.837

Table 9. Table S5: Speedups obtained for the stochastic Kuramoto model as a function of the number of orbits for system size N = 10 𝑁 10 N=10 , when using each of the compute devices listed in Table  3 . Speedups are calculated as the runtime for one core of the Intel Xeon 6142 CPU divided by the runtime of the compute device. In the case of the CPUs (2X6142 and 4790K), all cores were used.

Orbits	2X6142	I4I90K	W8100	P100
512	1.462	0.617	0.622	1.574
5120	2.895	1.269	1.671	2.814
25600	4.301	1.884	2.139	3.474
40960	4.738	2.039	2.298	3.771
81920	4.791	2.084	2.294	3.792
163840	4.905	2.163	2.350	3.934

Table 10. Table S6: Speedups obtained for the stochastic Kuramoto model as a function of the number of orbits for system size N = 15 𝑁 15 N=15 , when using each of the compute devices listed in Table  3 . Speedups are calculated as the runtime for one core of the Intel Xeon 6142 CPU divided by the runtime of the compute device. In the case of the CPUs (2X6142 and 4790K), all cores were used.

Orbits	2X6142	I4790K	W8100	P100
512	1.664	1.068	0.649	2.552
5120	3.796	1.796	1.718	4.289
25600	5.692	2.408	2.121	4.305
40960	6.073	2.500	2.188	4.373
81920	6.409	2.594	2.235	4.377
163840	6.652	2.661	2.291	4.503

Table 11. Table S7: Runtimes for the stochastic Kuramoto model as a function of the number of orbits for system size N = 5 𝑁 5 N=5 , using SODECL and parallel MATLAB implementations on all cores of the i7-4970K CPU. Runtimes are given in seconds (s). The last column (MATLAB/SODECL) shows the ratio of the runtimes of the MATLAB and SODECL implementations (the speedup). Speedup values greater than one indicate that SODECL was faster.

Orbits	SODECL	MATLAB	MATLAB/SODECL
512	0.728	0.131	0.180
5120	0.898	1.229	1.369
25600	1.612	6.149	3.813
40960	2.254	9.841	4.364
81920	3.800	19.588	5.154
163840	6.779	39.213	5.784

Table 12. Table S8: Runtimes for the stochastic Kuramoto model as a function of the number of orbits for system size N = 10 𝑁 10 N=10 , using SODECL and parallel MATLAB implementations on all cores of the i7-4970K CPU. Runtimes are given in seconds (s). The last column (MATLAB/SODECL) shows the ratio of the runtimes of the MATLAB and SODECL implementations (the speedup). Speedup values greater than one indicate that SODECL was faster.

Orbits	SODECL	MATLAB	MATLAB/SODECL
512	0.714	0.2966	0.415
5120	1.1383	2.8749	2.525
25600	3.1609	14.3615	4.543
40960	4.7469	23.0054	4.846
81920	8.784	45.8372	5.218
163840	16.939	91.1109	5.378

Table 13. Table S9: Runtimes for the stochastic Kuramoto model as a function of the number of orbits for system size N = 15 𝑁 15 N=15 , using SODECL and parallel MATLAB implementations on all cores of the i7-4970K CPU. Runtimes are given in seconds (s). The last column (MATLAB/SODECL) shows the ratio of the runtimes of the MATLAB and SODECL implementations (the speedup). Speedup values greater than one indicate that SODECL was faster.

Orbits	SODECL	MATLAB	MATLAB/SODECL
512	0.772	0.529	0.6861
5120	1.606	5.050	3.144
25600	5.256	25.322	4.817
40960	7.987	40.642	5.087
81920	15.279	80.467	5.266
163840	29.882	161.250	5.396

Equations

$$d{X}_i(t)=f_i(t,\mathbf{X}(t),\mathbf{p})dt+\sum _{j=1}^m{g}_{ij}(t,\mathbf{X}(t),\mathbf{p})d{W}_j(s),$$

$$X_i^{k+1}=X_i^k+f_i(t,{\mathbf{X}}^k,\mathbf{p})\Delta t+\sum _{j=1}^m{g}_{ij}(t,{\mathbf{X}}^k,\mathbf{p})\sqrt{\Delta t}{N}_j(0,1),$$

$$\displaystyle\begin{array}[]{lcc}r_{1}=\sqrt{-2\ln U_{1}}\cos(2\pi U_{2}),\\ r_{2}=\sqrt{-2\ln U_{1}}\sin(2\pi U_{2}).\end{array}$$

$$k=t_{\mathrm{total}}/(\Delta t\cdot s),$$

$$d{\theta }_i(t)={\omega }_i+\frac{K}{N}\sum _{j=1}^N\sin ({\theta }_j(t)-{\theta }_i(t))+p_{i}d{W}_i(t),$$

$${\theta }_i^{k+1}$$

$$r(t)e^{i\Phi (t)}=\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j(t)}.$$

Full text

Here, we present SODECL (https://github.com/avramidis/sodecl), a software library that utilises such devices to calculate multiple orbits of an SDE model. To evaluate the acceleration provided by SODECL, we compared the time required to calculate multiple orbits of an exemplar stochastic model when one CPU core is used, to the time required when using all CPU cores or a GPU. In addition, to assess scalability, we investigated how the model size affected execution time on different parallel compute devices.

Our results show that when using all 32 CPU cores of a high-end high-performance computing node, the task is accelerated by a factor of up to $\simeq$6.7, compared to when using a single CPU core. Executing the task on a high-end GPU yielded accelerations of up to $\simeq$4.5, compared to a single CPU core.

K****eywords stochastic differential equations  $\cdot$ CPU  $\cdot$ GPU  $\cdot$ HPC  $\cdot$ OpenCL

1 Introduction

1.1 Modelling dynamic systems using SDEs

Noise affects the behaviour of a vast number of physical and biological phenomena, such as gene expression [1, 2], the transmission of nerve impulses [3, 4] and the dynamics of electronic circuits [5, 6, 1]. Systems that are affected by noise are commonly modelled using sets of stochastic differential equations (SDEs) [7]. In such models, each equation represents the rate of change of a system variable (e.g. voltage) with time and can be dependent on one or more parameters (e.g. capacitance). Moreover, one or more equations are also affected by a noise process. The It$\hat{\rm o}$ form of a first-order coupled SDE system is

[TABLE]

where $t\geq t_{0}$ is time, $\mathbf{X}=\left(X_{1},\ldots,X_{n}\right)$ is the state vector, $\mathbf{X}_{0}=\mathbf{X}(t_{0})$ specifies the initial state of the system, $\mathbf{p}=\left(p_{1},\ldots,p_{d}\right)$ is a vector of parameter values, $\{W_{1}(t),\ldots,W_{m}(t)\}$ are independent, scalar Wiener processes, and the functions $\{f_{i}(t,\mathbf{X},\mathbf{p}):1\leq i\leq n\}$ and $\{g_{ij}(t,\mathbf{X},\mathbf{p}):1\leq i\leq n,1\leq j\leq m\}$ are the drift and diffusion coefficients, respectively [8, 9, 10]. Setting the diffusion coefficients $g_{ij}$ to 0 reduces the system to a set of deterministic, ordinary differential equations (ODEs).

1.2 Optimising the parameters of SDE models

By varying the parameters $\mathbf{p}$ of the SDE system (1), different classes of dynamical behaviour can be observed (e.g. transitions between bursting and spiking in neural models). However, in many cases, the parameters are difficult to measure experimentally. Consequently, optimisation methods are often employed to find the particular parameter combinations that most closely reproduce the experimentally-measured behaviour of the system of interest [11]. Performing this parameter optimisation step in a robust, automated fashion is a critical step in the construction and analysis of biological models [12, 13, 14, 15, 16, 17], because determining the optimal parameter values enables alternative models to be systematically ranked and experimentally testable predictions to be formulated [18, 19, 20, 21, 14, 22, 23, 16]. The assumptions made in constructing a given model can then be rigorously assessed, and insights obtained into how the model could be modified so as to improve the accuracy of its predictions [20, 21, 14, 22, 16].

However, optimising the parameters $\mathbf{p}$ of system (1) can incur a high computational cost, depending on the forms of $f_{i}$ and $g_{ij}$, the system size $n$, and the number of noise terms $m$. Large systems (high $n$ values), multiple noise processes (high $m$ values) and computationally demanding drift and diffusions terms (e.g. terms containing transcendental functions) can all increase the computational cost considerably. Furthermore, an additional cost is associated with the generation of the random numbers required to simulate the noise processes. The computational cost can also be significantly increased by the particular parameter optimisation method employed; for example, when using a population-based method – such as an evolutionary algorithm (EA) or a particle swarm optimiser (PSO) – the equations have to be integrated multiple times over different parameter values to explore the underlying fitness landscape [24, 25, 26]. Additionally, a resampling approach might be required to mitigate the effects of noise and uncertainty in model evaluation, thereby increasing the computational load further [27, 17].

In general, the evaluation of an SDE model for the optimisation methods mentioned above involves individual tasks with no interdependence or communication. These tasks can therefore be parallelised in a straightforward manner, and hence the optimisation problem is referred to as embarrassingly or pleasingly parallel [28].

1.3 Accelerating the optimisation process using multi-core and many-core devices

The significant computational demands of such embarrassingly parallel tasks means that high-performance computing (HPC) clusters (i.e. multiple connected computers) are required to obtain results within a reasonable time frame. Recently, competitively-priced CPU models that provide between 6 and 16 physical cores have become commercially available. In the case of workstation CPUs, there are now models with 32 cores.

Moreover, current GPU models, which are commonly used for graphics generation in personal computers and workstations, can also be leveraged for scientific computing [29, 30]. They contain a large number of processing units – each of which is relatively slow compared to the cores of a CPU – and are optimised for single instruction, multiple data (SIMD) processing. This makes them well-suited to large-scale parallelisation. Indeed, the collective processing power of a GPU can exceed that of a single CPU for certain tasks. Moreover, depending on the motherboard model, multiple GPUs can be installed on a single desktop/workstation computer, boosting the potential processing power even further.

In a previous study, we demonstrated that a desktop computer with a GPU was able to robustly optimise the parameters of a spiking neuron model to experimental eye movement data using a multi-objective EA [16]. As part of this work, we showed that the most computationally intensive part of the optimisation process was the numerical integration of the model over the multiple parameter combinations comprising the EA population at each generation [16]. This task can be executed in an embarrassingly parallel fashion, enabling the current trend for heterogenous HPC architectures to be exploited. Accordingly, in our spiking neuron model, numerical integration was executed on the GPU, while all the EA operations (mutation, crossover, fitness evaluation etc.) were executed on the CPU [16]. Utilising the GPU in this manner yielded a speedup of up to $\approx$20, compared with a high-end CPU.

1.4 CUDA and OpenCL

Programs that are executed on a GPU have to be written using a specific programming framework. The two most commonly used frameworks are CUDA (Compute Unified Device Architecture) and OpenCL (Open Computing Language) [31]. The CUDA framework is a proprietary architecture specifically designed to be run on NVIDIA compute devices. By contrast, OpenCL is a royalty-free standard for general purpose parallel programming across CPUs and GPUs, giving software developers portable and efficient access to the power of these heterogeneous processing platforms [32]. OpenCL includes a cross-platform intermediate language for writing functions (kernels), which are executed on OpenCL-supported devices, together with an application programming interface (API) that is used to coordinate the parallel computations across these devices. A simple OpenCL program involves the following steps: identifying the OpenCL device; compiling the OpenCL code that will be run on the device; copying the data to the device; performing the computation; copying the results back from the device.

A software library equipped with OpenCL functionality can use CPUs and GPUs, removing the need to develop different libraries for each processor type.

1.5 A new library for numerically integrating SDE models on OpenCL-supported devices

Here, we present SODECL, a C++ library that uses OpenCL to calculate multiple orbits of an SDE (or ODE) system in an embarrassingly parallel way. We focus on presenting the library for solving SDE models, for which orbits are calculated using the Euler-Maruyama method [33, 9]. ODE orbits can be computed using any of the following integration methods: Euler, Runge-Kutta, Implicit Euler or Implicit Midpoint. SODECL has previously been used to fit both ODE and SDE versions of an oculomotor control model to experimentally recorded timeseries [34, 16].

In the following sections, we describe the design principles of the SODECL library, its organisation and the numerical algorithms used. Moreover, we describe the experimental protocols that were implemented to measure the execution speed of SODECL for different compute devices (i.e. CPUs and GPUs), and to assess the numerical stability and accuracy of our SDE solver. As part of the speed test experiments, we compare the performance of SODECL on a high-end multi-core CPU with a MATLAB program that also calculates multiple SDE orbits in parallel using the same integration method (i.e. Euler–Maruyama). Lastly, we outline directions for future development of the library.

2 The SODECL library

2.1 Design principles

We chose to implement SODECL using OpenCL, because it does not pose significant constraints on the computer hardware used. As mentioned above, whilst a user must run CUDA executables on an NVIDIA GPU, OpenCL executables can run on both NVIDIA GPUs and AMD/Intel GPUs. Moreover, OpenCL executables can run on Intel and AMD CPUs. This gives OpenCL a key advantage over CUDA in terms of disseminating research methods and results, and promoting scientific collaborations.

We designed SODECL to be relatively portable, modular, and easy to extend. We tested our library on the following operating systems: Windows 10 64-bit; Ubuntu 18.04 64-bit; macOS Sierra (10.12). The SODECL library can easily be integrated into any C++ source code by including the main header of the library and linking it to an OpenCL library. Moreover, SODECL can be extended by adding OpenCL functions for other SDE integration methods (e.g. the Milstein method [8]). The source code is released under the MIT License and is under version control with git at https://github.com/avramidis/sodecl.

2.1.1 The integration method

SODECL uses the Euler-Maruyama method [33, 9] to integrate the equations in (1). The Euler-Maruyama approximation $\left\{\mathbf{X}^{k}=\left(X_{1}^{k},\ldots,X_{n}^{k}\right):k\geq 1\right\}$ to the true solution $\{\mathbf{X}(t):\mathbf{X}(t_{0})=\mathbf{X}_{0},t\geq t_{0}\}$ is defined by the recursion

[TABLE]

where $\mathbf{X}^{1}=\mathbf{X}_{0}$, $\Delta t$ is the time step (so that $\mathbf{X}^{k}$ is the approximation to $\mathbf{X}((k-1)\Delta t)$ for $k>1$) and $N_{j}(0,1)$ are independent, normally distributed random variables with zero mean and unit variance.

2.1.2 Noise generator

To generate the noise variables $N_{j}(0,1)$ used in each step of (2), SODECL utilises the Random123 library [35]. Random123 generates uniformly distributed random numbers, which are then converted to normally distributed random numbers using the Box-Muller algorithm [36]:

[TABLE]

Here, $\{U_{1},U_{2}\}$ are two independent, uniformly distributed random variables over [0,1] and $\{r_{1},r_{2}\}$ are two independent, normally distributed random variables with zero mean and unit variance. In the case where $j>2$, Eqs. (5) are called multiple times.

2.1.3 OpenCL code generation

SODECL creates and builds the OpenCL source code string at runtime. This allows the user to change the $f_{i}$ and $g_{ij}$ functions in (1) without the need to compile the source code. Moreover, SODECL allows the implementation of wrappers – e.g., for MATLAB or Python, the use of a script/function that runs the SODECL executable, or the use of a MATLAB .mex file without it being recompiled. The procedure for formulating the OpenCL source code string is performed in the four steps outlined below.

Step 1

Initially, SODECL appends the definitions of five named constants to the beginning of the empty OpenCL source code string (see Table 1). The values of these constants depend on the general properties of the SDE system (numbers of equations, noise terms and parameters), the integrator’s time step and the number of steps to be executed by the integrator in each OpenCL kernel call. These values are passed to the SODECL library by the user. The remainder of the OpenCL source code string is generated using the code located in three separate files, described in steps 2-4 below.

Step 2

The first file, integrator_caller.cl, contains the kernel function that is called by the host to run on an OpenCL device. This function is defined in Figure 1. The parameters of the kernel function are pointers that show the location of the arrays in the global memory of the device. These arrays contain the time points of the SDE system at which to approximate the solution, the values of the SDE system’s dependent variables and the parameter combinations for which the different orbits are to be calculated. One additional array of counters is passed to the kernel for use by the Random123 library. The kernel generates the independent, normally distributed numbers required for the integration of the SDE system. Subsequently, the same kernel calls the function that integrates the SDE system for one time step (see step 3). The total number of calls is defined by the constant numsteps. Having integrated the SDE system numsteps times, the kernel copies the last integrated state values to the global memory of the OpenCL device, to be used in the next kernel step and accessed by the host.

Step 3

The second file used in the creation of the OpenCL source is named stochastic_euler.cl. It contains the function system_integrator (defined in Figure 2), that implements the Euler-Muryama integration method. system_integrator itself calls functions that calculate the deterministic and stochastic components of the SDE system (see step 4). Once these have been evaluated, stochastic_euler.cl calculates the values of all dependent variables at the next time step.

Step 4

The third file contains two functions called sode_system and sode_system_stoch, which calculate the drift and diffusion terms ${f_{i}}$ and $g_{ij}$ in the numerical scheme (2), respectively. Their definitions are shown in Figure 3 and they have to be specified by the user (an example is shown in Figure 5).

2.1.4 SDE integrator execution algorithm

The logical flow of the SODECL integration execution algorithm is predominately a loop. During each loop iteration, the corresponding values of the dependent variables for all orbits are copied to the host device, and the copied values are saved in parallel to an array in preparation for the next kernel call. The OpenCL kernel is called to integrate the model for numsteps steps of length m_dt, which means that the output of the integration is saved to the array every (numsteps $\times$ m_dt) of integration time. This allows the data to be stored with a frequency of $1/$(numsteps $\times$ m_dt) Hz. Once the integration over the required time span is completed for all parameter combinations, the algorithm ends. The number of loop iterations, $k$, is calculated using the equation

[TABLE]

where $t_{\rm total}$ is the time span of the integration, $\Delta t$ is the integrator time step and $s$ is the number of steps the integrator performs in each OpenCL kernel call (i.e. the number of parallel operations, which could include multiple parameter combinations and/or initial conditions). The values of $t_{\rm total}$, $\Delta t$ and $s$ are supplied by the user.

The storing frequency of the results plays an important role in the speedup provided by the GPU compared to the CPU. This is due to the bottleneck caused by the transfer to and from the GPU memory, and the time required for storing the data to the host memory for further analysis. A low frequency (250-1000 Hz) avoids substantial data storage in the host memory. Also, the model size is an important factor, since it sets the number of model integrations that are to be saved on the computer’s RAM. Moreover, a larger number of returned data points can increase the analysis time; a low number, however, may be insufficient for an accurate analysis [16].

2.1.5 Python interface

The Python package interface function for SODECL is defined in Figure 4, whilst Table 2 describes the arguments passed to the SODECL executable by the function. The Python library pybind11 (https://github.com/pybind/pybind11) was used to implement the interface. A MATLAB function interface is also provided with SODECL, although it should be noted that this is currently experimental.

2.2 SODECL speed evaluation

2.2.1 Protocol

To evaluate the execution speed of SODECL, we used a stochastic variant of the well-established Kuramoto model for describing the dynamics of a population of weakly coupled phase oscillators [37, 38, 39, 40, 41]. The Kuramoto model is commonly used to model synchronisation of biological oscillators [39, 40], and is relatively straightforward to implement and customise (e.g. the oscillator population size can be easily varied). The SDE Kuramoto system we used here has the form

[TABLE]

where $\{\theta_{1},\ldots,\theta_{N}:-\pi\leq\theta_{i}<\pi\}$ specify the phase of each oscillator, $\omega_{i}$ is the free-running frequency of the $i$th oscillator, $\{W_{1}(t),\ldots,W_{N}(t)\}$ are independent, scalar Wiener processes, $p_{i}$ controls the strength of the noise effect on the $i$th oscillator and $K$ is the strength of the mean-field coupling. The system parameters are thus $\mathbf{p}=(K,\omega_{1},\ldots,\omega_{N},p_{1},\ldots,p_{N})$.111In terms of Eq. (1), $X_{i}=\theta_{i}$, $n=m=N$, $f_{i}=\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i})$ and $g_{ij}=p_{i}\delta_{ij}$. The initial value of each phase variable $\theta_{i}$ was taken to be uniformly distributed over $[-\pi,\pi)$, the frequencies $\omega_{i}$ to be uniformly distributed over [0.01, 0.03] and the noise strengths $p_{i}$ to be uniformly distributed over [0.001, 0.003]. The following equations define the integration scheme for (7) using the Euler-Maruyama method (cf. (2) above):

[TABLE]

We set the coupling constant $K$ in (7) to 1 and integrated the model for 400 s with a time step $\Delta t$ of $0.05$ s. Moreover, we used double-precision floating-point data type for better accuracy. The SODECL implementation of the stochastic Kuramoto model is shown in Figure 5.

We compared the execution speed of SODECL when using different multi-core and many-core compute devices with the execution speed when using only one core of an HPC CPU. The compute devices used to evaluate SODECL performance are listed in Table 3. The CPU with ID I4790K is a high-end desktop CPU, whereas the 2X6142 is an HPC node with two many-core CPUs. The W8100 and P100 are workstation/server GPUs, with high double-precision compute capabilities required for scientific computing. To gain a better understanding of how the SODECL runtime is affected by the size of the system, we varied the number of equations in the Kuramoto model (i.e. we varied the number of oscillators $N$ in the network). Moreover, we examined how the number of orbits $M$ being integrated affects the performance of SODECL on the different compute devices. Eight independent runs of the solver were executed for each choice of $N$ and $M$.

We also compared the execution speed of SODECL on a CPU against an equivalent parallel MATLAB implementation running on the same CPU. In order to further accelerate this implementation, described in Figure 6, we used the MATLAB Coder to convert MATLAB functions into C++ code. The parallelisation of MATLAB was performed automatically with OpenMP (Open Multi-Processing), an application programming interface for parallel applications that uses the cores of CPUs and coprocessors [42].

2.2.2 Results

The comparative execution speeds and speedups obtained using SODECL with different compute devices, numbers of orbits and system sizes are shown in Figure 7 (see Tables S1-S6 for the corresponding numerical values). For each compute device, the speedup was calculated relative to one core of the 2X6142 node (i.e. as the ratio of the 2X6142 single core runtime to the compute device runtime). The local group sizes used for each device and system size ($N$) are shown in Table 4. The local group size is the number of work items that will run in parallel and can communicate. In the case of embarrassingly parallel applications, although there is no communication, the group size can have an effect on the runtime. The values shown were found by doing exploratory tests, since – to our knowledge – there is no analytical formula for calculating the optimal local group size.

For the CPUs tested, the results indicate that when using all 32 cores of the 2X6142 node, SODECL is between $\simeq$1.12 and $\simeq$6.65 times faster than one core of the same node across different system sizes and orbit numbers. By contrast, when using all cores of the I4790K CPU and comparing its speed to one core of 2X6142, we observe speedups ranging from $\simeq$0.69 to $\simeq$2.66 as the system size is increased from 5 to 15 and the number of orbits is increased from 512 to 163840. The greater speedup obtained with 2X6142 is due to its larger number of cores and newer architecture.

For the GPUs tested, SODECL yields speedups in the range 0.6-4.5 compared to using one core of 2X6142, depending on the size of the Kuramoto model and the number of orbits. Larger numbers of model equations and orbits allow for greater use of GPU resources, which translates to improved speedup. When comparing the GPU runtimes to those obtained when all cores of 2X6142 are utilised, it can be seen that the GPUs are generally slower. This could, for example, be due to higher initialisation time overheads.

Finally, Figure 8 compares the runtimes obtained with the SODECL and MATLAB implementations of the Kuramoto model solver on the I4790K CPU (the corresponding numerical values are given in Tables  S7-S9). The MATLAB implementation is only faster than SODECL for the smallest system size and lowest number of orbits (see Figures  8B, E & H). In all other cases, SODECL is up to 5.8 times faster.

2.3 Accuracy and numerical stability tests

2.3.1 Protocol

To assess the accuracy and numerical stability of SODECL, we generated a set of further simulations with $N=100$, $\omega_{i}\sim\mathcal{U}(0.2,0.4)$, $\phi_{i}(0)\sim\mathcal{U}(-\pi,\pi)$ and $p_{i}\sim\mathcal{U}(0.01,0.03)$, for two values $K_{1}=0.02$ and $K_{2}=0.2$ of the global coupling strength $K$. These values were chosen to lie either side of the Kuramoto transition that occurs in the deterministic system in the continuum limit ($N\to\infty$) at $K=K_{c}=\frac{2}{5\pi}\approx$ 0.1273. As $K$ is increased through this critical value, macroscopic mutual entrainment (MME) – a collective synchronised rhythm with a common frequency – is observed in the oscillator population [37, 38]. For both values of $K$, 64 independent integrations were carried out over 400 s (around 19 cycles of the mean period 20.94 s), using the following integration time steps: $\Delta t=\frac{2^{l-5}}{5},1\leq l\leq 5$. For each realisation, the degree of synchronisation of the simulated oscillator population was quantified using the complex order parameter [38, 41], defined below:

[TABLE]

In the above, the radius $r(t)$ measures the phase coherence and $\Phi(t)$ measures the collective phase, with $r$ values of 0 and 1 corresponding to complete desynchronisation and complete synchronisation of the population, respectively [38, 41]. In the continuum limit of the deterministic model, increasing the coupling strength $K$ through the critical value $K_{c}$ causes the steady state coherence $\lim_{t\rightarrow\infty}r(t)$ to increase rapidly from 0 to 1 as the population becomes globally synchronised [38, 39, 41]. For finite $N$, values of $K$ less than $K_{c}$ yield $\mathcal{O}(N^{-1/2})$ fluctuations in $r(t)$, while values of $K$ greater than $K_{c}$ result in $r(t)$ saturating at a value $\lim_{t\rightarrow\infty}r(t)<1$, with $\mathcal{O}(N^{-1/2})$ fluctuations [38].

2.3.2 Results

Figure 9 plots how the mean phase coherence $\langle r(t)\rangle$ and the standard deviation of the phase coherence $\sigma\left(r(t)\right)$ vary with time for the two $K$ values when $\Delta t=0.05$. As expected, for $K=K_{1}<K_{c}$, $\langle r(t)\rangle$ exhibits small-amplitude oscillations close to 0, while for $K=K_{2}>K_{c}$, $\langle r(t)\rangle$ asymptotes to a value close to 1. Kymographs showing the temporal evolution of the population phases for a typical realisation generated with each $K$ value are plotted in Figure 10. These are consistent with the phase coherence plots – the oscillators are desynchronised throughout the integration interval for $K=K_{1}$, while for $K=K_{2}$, the oscillators synchronise after around 50 s and remain synchronised thereafter. Finally, Figure 11 plots how the values of $\langle r(t)\rangle$ and $\sigma\left(r(t)\right)$ at the end of the integration interval $t_{MAX}=400$ vary with the integration timestep $\Delta t$. It can be seen that for both $K$ values, $\langle r(t_{MAX})\rangle$ and $\sigma\left(r(t_{MAX})\right)$ exhibit a weak dependence on $\Delta t$, suggesting that the integration algorithm is numerically stable [10].

3 Conclusions

We have described the design and performance of SODECL, an open source C++ library that uses OpenCL to calculate multiple orbits of a SDE system in an embarrassingly parallel way. The advantage of SODECL is that it can use the parallel capabilities of multi-core CPUs and GPUs to reduce computation time, thereby allowing large numbers of simulations to be performed simultaneously. This capability is critical when using population-based metaheuristics for data-fitting, such as genetic algorithms, as it facilitates the optimisation and analysis of SDE models in a practical time frame [43]. SODECL was designed to compile and run on Windows, Linux and macOS operating systems, to be user friendly, fast and fairly easy to extend.

Although execution speed was not the main priority in our design of SODECL, we showed that it was faster in almost all cases when integrating the stochastic Kuramoto system compared to a parallel MATLAB implementation. We also showed that when SODECL utilises a GPU or many-core CPUs, there is considerable speedup compared to a high-end desktop CPU. However, the speedup can be affected by the type and size of the SDE system being integrated. It would therefore be instructive to obtain performance benchmarks for other canonical stochastic models.

In our previous work [16], SODECL achieved a speedup of $\simeq$20 on the W8100 GPU compared to all cores of the I4790K CPU. By contrast, the greatest speedup obtained with the same GPU against I4790K in this study was $\simeq$1.16. We believe that the difference in speedups between our previous and current studies is that in [16], the SODECL OpenCL code could not be vectorised efficiently to facilitate better utilisation of the CPU, for the following reasons. Firstly, the system was a deterministic model with six heterogeneous equations (i.e. equations containing different functional forms). Similar equations, if vectorised, could run in parallel more efficiently. Secondly, two of these equations included a case mathematical function. Case functions are implemented with one or more if statements that cause branch divergence. Most of the time, branch divergence cannot be predicted and a compiler cannot vectorise the code for branches in which each case implements a different mathematical operation. Thirdly, for the numerical integration of the model, the implicit, midpoint Euler method was used. This method includes a loop with a variable iteration number. All of the above factors hinder the compiler from efficiently vectorising the execution of the code on a parallel device.

There are a number of aspects of SODECL that we did not examine, which could provide interesting avenues for future work. Firstly, we only measured SODECL execution speed using double-precision, rather than single-precision, floating-point data type. Single-precision can be used if high accuracy is not necessary, and has the advantage that compute devices yield significantly better execution speeds when using lower precision. Secondly, we did not measure the OpenCL initialisation and finalisation time in SODECL. This would indicate whether different compute devices require more time for these tasks, perhaps explaining why the MATLAB implementation is comparable to SODECL for small numbers of orbits. Thirdly, we did not evaluate the effects of branch divergence on different compute devices using the Euler-Maruyama method. This potentially would have shown which devices are less affected. However, to assess the extent to which performance is affected by equations with many branches, a more extensive study would need to be carried out involving a broader range of SDE models. Lastly, we did not examine the memory requirements for different orbit numbers and model sizes. A follow-up study could examine the impact that the high memory requirements of larger models (i.e. more than 100 equations) could have on the performance of different compute devices.

We remark that the future development of SODECL would involve incorporating additional methods for numerically integrating SDE systems, such as the Milstein and Runge-Kutta schemes [8]. Here, we have implemented the Euler-Maruyama numerical integration scheme, which despite its simplicity, is still widely used, in part owing to the relative ease of implementation. We note, however, that the inefficiency of the Euler-Maruyama method compared to more sophisticated numerical schemes is offset by the significant acceleration conferred by parallelising the computation of trajectories across multiple CPU cores or GPU compute units. The incorporation of additional solvers would involve modifying the SODECL system_integrator function so as to implement the chosen schemes (see Fig. 2).

Also, additional optimisations could be explored for both GPUs and CPUs, such as using newer versions of OpenCL. This could potentially further accelerate SODECL on certain devices. Furthermore, we believe that splitting the integration algorithm over multiple OpenCL kernels could potentially increase the execution speed.

Acknowledgments

This work was financially supported by the Engineering and Physical Sciences Research Council (grant numbers EP/K040987/1, EP/N017846/1 and EP/N014391/1). Computational experiments were performed using resources provided by the Cambridge Service for Data Driven Discovery (CSD3, https://www.hpc.cam.ac.uk/high-performance-computing) and the University of Exeter High-Performance Computing (HPC) facility. CSD3 is operated by the University of Cambridge Research Computing Service (http://www.hpc.cam.ac.uk), funded by EPSRC Tier-2 capital grant EP/P020259/1, the STFC DiRAC HPC Facility (http://www.dirac.ac.uk) and the University of Cambridge. CSD3 and DiRAC are part of the UK National e-Infrastructure.

Supplementary Tables