STROOPWAFEL: an algorithm to efficiently sample rare astrophysical events.

Floor Broekgaarden 1

  • 1 Center For Astrophysics | Harvard & Smithsonian, United States


Gravitational-wave observations of binary black hole mergers are rapidly providing new insights into the physics of massive stars and the evolution of binary systems. Making the most of expected near-future observations for understanding stellar physics will rely on comparisons with binary population synthesis models. However, the vast majority of simulated binaries never produce binary black hole mergers, which makes calculating such populations computationally inefficient.
In this meeting I will present our adaptive importance sampling algorithm, STROOPWAFEL, that we wrote to improve the computational sampling efficiency of population studies of rare events. I will present its performance compared to traditional Monte Carlo sampling from birth distributions and will discuss the similarities of the code with playing the board game battleships.
At the end of the presentation I will discuss some statistical challenges that we are currently facing in our effort to further optimize the STROOPWAFEL code, for which I would love to get input and ideas from the audience.


Illustration of a violent merger between two binary black holes billion light years away. The gravitational-waves produced by such an event were for the first time obeserved by the detectors of LIGO on September 14, 2015 as they passed by Earth. Video credit: LIGO/SXS/R.Hurt and T. Pyle

Gravitational-wave observations of double compact object mergers are rapidly providing new insights into the physics of massive stars and the evolution of binary systems [1,2,3].  The distributions of the inferred masses and spins of the black holes and neutron stars contain valuable information about their origin. Distinguishing different theories for their formation and learning about the complex physical processes that govern the lives of their possible massive-star progenitors requires comparing observed populations with theoretical predictions.

However, theoretical simulations of the population of merging double compact objects are challenging because gravitational-wave events represent an extremely rare outcome of binary evolution. From a thousand massive binary systems typically only of order one, or less, yields a double compact object [4]. A meaningful comparison with population observations requires simulating a statistically significant sample of events.  Modeling binaries over cosmological volumes whilst exploring the uncertain physics requires running computationally expensive binary population synthesis simulations with billions of binaries that cost millions of CPU hours.

To make it feasible to simulate such large numbers of systems, all present-day simulations pay a high price. In many studies computational speed is ensured by using highly approximate algorithms that treat the physical processes in a simplified way. Another way to keep the computational costs reasonable is to limit the total number of simulations, restricting the exploration of the impact of the uncertain physical input assumptions beyond a few variations ​


To overcome this issue we developed the algorithm STROOPWAFEL: Sampling The Rare Outcome Of Populations With AIS* For Efficient Learning. We have designed STROOPWAFEL to improve the efficiency of simulations of rare astrophysical events, and so to enable accurate simulations of populations of extremely rare events at reasonable computational cost.


*AIS: Adaptive Importance Sampling

STROOPWAFEL works similar to playing the game Battleships

The STROOPWAFEL method explained by comparison with playing the game Battleships.


The number of simulated binaries NT that forms a black hole - neutron star system as a function of the total number of binaries Nbinaries sampled for the traditional sampling method (gray dashed line) and the sampling method presented in this study (solid coloured line).  The duration of the exploratory phase is shown with a hashed gray area. In the background the standard Poisson fractional uncertainties of 0.3, 1 and 3% are shown with dashed lines.

We find that STROOPWAFEL can increase the number of binaries of a target population that are found in a simulation by factors of about 25 – 200 when using our STROOPWAFEL sampling algorithm compared to simulations with traditional sampling. The  figure on the left showcases this for the black hole-neutron star simulation. The figure presents the number of systems formed from the target population as a function of total number of sampled binaries, for both our sampling method and traditional birth distribution Monte Carlo sampling. For this simulation the gain is a factor 53.  

Contour plots of the locations in initial mass  m1,i and initial separation ai space of the black hole-neutron star mergers found in each simulation when using traditional birth distribution Monte Carlo sampling (left panel) and our sampling method STROOPWAFEL  (right panel). Contours represent a constant density of binaries of the target population found per unit area in initial mass-separation space. The colour gradient indicates the number of samples per area ∆S, the size of which is shown with a black rectangle. If the density is below the level of our lowest contour we plot the individual points. The total number of hits NT found in each simulation is quoted in parentheses. The metallicity assumed in all simulations is Z = 0.001.

The increase in computational efficiency from STROOPWAFEL leads to finding substantially more events of the target population which naturally enables a much higher resolution mapping of both the input and output parameter spaces. An example is shown above, which demonstrates how the parameter space is explored in far greater detail using our sampling method compared to traditional birth distribution Monte Carlo sampling. The figure shows the location of the target population in the initial parameter space of primary mass m1,i and separation  ai at birth. With our  STROOPWAFEL sampling method we obtain more detailed contours and more contour levels that map the initial parameter space with higher resolution. This leads to better knowledge of the initial conditions of a binary system that yield a binary of the target population. Physically the structures seen in the input parameter space correspond to the assumed physics of the different formation channels leading to compact-object mergers. More details are discussed in  [5,6,7].

Future outlook

The STROOPWAFEL code is publicly available. The algorithm is described and published in Broekgaarden et al., 2019. All data used in this study is publicly available, as well as jupyter notebooks to reproduce all figures and the code to reproduce the paper and tables.   

STROOPWAFEL is currently implemented in COMPAS and has already been used in studies including in work on population predictions for common-envelope events in the formation of double neutron stars by Vigna-Gómez et al., 2019 and for a study in forming binary black hole mergers in the upper black hole mass gap with super-Eddington accretion by van Son et al, 2020. Several other studies using STROOPWAFEL are in prep. 

At the moment Lokesh Khandelwal and Floris Kummer are further improving the STROOPWAFEL code.

Lokesh Khandelwal is finishing the publication of an object-orientated public version of STROOPWAFEL that is flexible and easy to use for any population synthesis code (and other sampling code). In addition, Lokesh Khandelwal implemented independently many essential improvements including more flexible sampling, improved boundary conditions for the adaptive sampling distribution and an extension to STROOPWAFEL sampling in 11 dimensions (instead of three). Lokesh is currently preparing his work for publication in a scientific journal.  Contact:

Floris Kummer is independently rewriting STROOPWAFEL such that it can take into account observational biases. In addition Floris Kummer is working on the integration of both STROOPWAFEL and cosmic star formation history models to population synthesis simulations to make predictions for gravitational-wave observations whilst improving to the metallicity grid-sampling problem. This is an important issue in population synthesis predictions where most current simulations predictions suffer in resolution as an effect from sparse grid-based sampling in birth metallicities. Contact: