SiSyPHE: A Python package for the Simulation of Systems of interacting mean-field Particles with High Efficiency

Over the past decades, the study of systems of particles has become an important part of many research areas, from theoretical physics to applied biology and computational mathematics. One of the main motivations in mathematical biology is the modelling of large animal societies and the emergence of complex patterns from simple behavioral rules, e.g., flocks of birds, fish schools, ant colonies, etc. In the microscopic world, particle systems are used to model a wide range of phenomena, from the collective motion of spermatozoa to the anarchical development of cancer cells. Within this perspective, there are at least three important reasons to conduct large scale computer simulations of particle systems. First, numerical experiments are essential to calibrate the models and test the influence of each parameter in a controlled environment. For instance, the renowned Vicsek model (Vicsek et al., 1995) is a minimal model of flocking, which exhibits a complex behavior, studied numerically in particular in (Chaté et al., 2008). Secondly, particle simulations are used to check the validity of macroscopic models that describe the statistical behavior of particle systems. These models are usually based on partial differential equations (PDE) derived using phenomenological considerations that are often difficult to justify mathematically (Degond et al., 2021; Degond & Motsch, 2008; Dimarco & Motsch, 2016). Finally, inspired by models in biology, there is an ever growing literature on the design of algorithms based on the simulation of artificial particle systems to solve tough optimization problems (Grassi & Pareschi, 2020; Kennedy & Eberhart, 1995; Pinnau et al., 2017; Totzeck, 2021) and to construct new more efficient Markov Chain Monte Carlo methods (Cappé et al., 2004; Clarté et al., 2021; Del Moral, 1998, 2013; Doucet et al., 2001). The simulation of systems of particles is also at the core of molecular dynamics (Leimkuhler & Matthews, 2015), although the present library is not specifically written for this purpose. The SiSyPHE library builds on recent advances in hardware and software for the efficient simulation of large scale interacting mean-field particle systems, both on the GPU and on the CPU. The versatile object-oriented Python interface of the library is designed for the simulation and comparison of new and classical many-particle models of collective dynamics in mathematics and active matter physics, enabling ambitious numerical experiments and leading to novel conjectures and results.


Summary
Over the past decades, the study of systems of particles has become an important part of many research areas, from theoretical physics to applied biology and computational mathematics. One of the main motivations in mathematical biology is the modelling of large animal societies and the emergence of complex patterns from simple behavioral rules, e.g., flocks of birds, fish schools, ant colonies, etc. In the microscopic world, particle systems are used to model a wide range of phenomena, from the collective motion of spermatozoa to the anarchical development of cancer cells. Within this perspective, there are at least three important reasons to conduct large scale computer simulations of particle systems. First, numerical experiments are essential to calibrate the models and test the influence of each parameter in a controlled environment. For instance, the renowned Vicsek model (Vicsek et al., 1995) is a minimal model of flocking, which exhibits a complex behavior, studied numerically in particular in (Chaté et al., 2008). Secondly, particle simulations are used to check the validity of macroscopic models that describe the statistical behavior of particle systems. These models are usually based on partial differential equations (PDE) derived using phenomenological considerations that are often difficult to justify mathematically (Degond et al., 2021;Degond & Motsch, 2008;Dimarco & Motsch, 2016). Finally, inspired by models in biology, there is an ever growing literature on the design of algorithms based on the simulation of artificial particle systems to solve tough optimization problems (Grassi & Pareschi, 2020;Kennedy & Eberhart, 1995;Pinnau et al., 2017;Totzeck, 2021) and to construct new more efficient Markov Chain Monte Carlo methods (Cappé et al., 2004;Clarté et al., 2021;Del Moral, 1998, 2013Doucet et al., 2001). The simulation of systems of particles is also at the core of molecular dynamics (Leimkuhler & Matthews, 2015), although the present library is not specifically written for this purpose. The SiSyPHE library builds on recent advances in hardware and software for the efficient simulation of large scale interacting mean-field particle systems, both on the GPU and on the CPU. The versatile object-oriented Python interface of the library is designed for the simulation and comparison of new and classical many-particle models of collective dynamics in mathematics and active matter physics, enabling ambitious numerical experiments and leading to novel conjectures and results.

Statement of need
A major difficulty in the simulation of systems of particles is the high computational cost, typically quadratic in the number of particles, which prevents large scale experiments. The implementation of SiSyPHE is based on recent libraries originally developed for machine learning purposes to significantly accelerate tensor (array) computations, namely the PyTorch package (Paszke et al., 2019) and the KeOps library (Charlier et al., 2021). On a GPU, the SiSyPHE library speeds up both traditional Python and low-level implementations by one to three orders of magnitude for systems with up to several millions of particles.
In addition, to the best of our knowledge, only model-specific packages such as Motsch (2016) are available. The SiSyPHE library includes, within a common framework, the implementation of many classical models and their variants as well as recent models for which no implementation was previously available. All the models detailed in the Example gallery of the documentation are directly taken from the literature on collective dynamics in mathematics and active matter physics. Moreover, the SiSyPHE library is designed in such a way that new custom models can easily be added in order to facilite the study and comparison of models from a research perspective.
The development of the SiSyPHE library was initially motivated by the study of body-oriented particles (Degond et al., 2019). The (formal) derivation of a macroscopic PDE model from the particle system has lead to a novel conjecture which postulates the existence of a class of so-called bulk topological states in (Degond et al., 2021). The quantitative comparison between this theoretical prediction and the numerical simulation of the particle system in a suitable regime (with more than 10 6 particles) has confirmed the existence of these new states of matter. The study of their physical properties which are observed in the numerical experiments but not readily explained by the PDE model is an ongoing work.

A typical example
A typical model that is implemented in the SiSyPHE library is the variant of the Vicsek model introduced by Degond & Motsch (2008) and defined by the system of 2N Stratonovich Stochastic Differential Equations where the kernel K : [0, +∞) → [0, +∞) is a smooth nonnegative function vanishing at infinity which models the visual perception of the particles; in the Vicsek model, the vision of the particles depends on the distance between them. With the target given by Equation 2, each particle tries to adopt the average orientation of its neighbors, which is a typical flocking behavior.
On a computer, the time-continuous system given by Equation 1 needs to be discretized first. For the Vicsek model, a natural discretization method is the (geometric) Euler-Maruyama scheme (Kloeden & Platen, 1992;Piggott & Solo, 2016). In general, the discretization method depends on the model considered as illustrated in the Example gallery. Then, at each time step, the most expensive operation is the computation of the target given by Equation 2, which requires O(N ) operations for each of the N particles. The total simulation cost is thus O (N 2 T ) where T is the total number of iterations. Within the framework of the KeOps library on which SiSyPHE is based, the computation of the target Equation 2 is called a kernel operation, which is efficiently carried out using a symbolic definition of the N × N interaction matrix whose (i, j)-entry is K(|X j t − X i t |). The computation of the target is then understood as a symbolic matrix-vector product between the interaction matrix and the vector of orientations.