Background
I'm attempting to write a physics simulation code, one portion of which involves simulating the triggering system of some equipment. The equipment works as follows: environmental noise (which we may model as a Gaussian-distributed random variable) is taken in as data in numChannels
independent streams at some sampling rate samplingRate
. Any time at least channelThreshold
channels have at least one sample exceeding a threshold threshold
within a window of size windowSize
samples, a "trigger" is recorded, and the equipment stops recording data for a duration corresponding to writeDelay
samples (so that it may devote computing resources to recording and analyzing the data). In hardware this is done using a system of integrating diodes. The goal of the simulation is to imitate this on several seconds worth of mock data, so as to determine the triggering rate as a function of the parameters of environmental noise (i.e. of the Gaussian distribution).
In practice, the equipment has a sampling rate north of 40 GS/s
, hence it is impractical to store in memory several seconds worth of data. To solve this, I've concocted an equivalent algorithm which is both far more memory efficient and somewhat faster. Unfortunately, however, it is still much too slow to iterate over a useful amount of moch data (running on a modern computing cluster, it takes on the order of 2.5
days to run on a single second, and I'd like to run on multiple seconds hundreds or thousands of times using different environmental noise parameters).
Code Walkthrough
int Antenna::triggerRate(double threshold, double temperature){
std::mt19937 gen(time(0));
std::normal_distribution<double> dist(0, vrms = _vrms * sqrt(temperature));
std::vector<boost::dynamic_bitset<>>
window(numChannels, boost::dynamic_bitset<>(windowSize));
auto clear = [&](){
for(auto &channel: window)
for(int i = 0; i < windowSize; ++i)
channel.set(i, fabs(dist(gen)) > threshold);
};
clear();
int numContributingChannels = 0, numTriggers = 0;
for(int i = 0, imod = 0; i < samplingRate; ++i, imod = (imod + 1) % windowSize){
numContributingChannels = 0;
for(auto &channel: window)
if(channel.any())
++numContributingChannels;
if(numContributingChannels > channelThreshold){
++numTriggers;
i += writeDelay;
clear();
continue;
}
for(auto &channel: window)
channel.set(imod, fabs(dist(gen)) > threshold);
}
return numTriggers;
}
Note that I've excluded the Antenna
class, which is quite large and somewhat messy. I don't believe that it is important here, but I would be happy to provide it.
The code works as follows: at any given time, only a window of size windowSize
worth of data is stored in memory. Moreover, I do not store the actual data, but instead store bitsets whose bits correspond to samples, and whose values correspond to whether or not a given sample exceeds the threshold threshold
. The vector window
holds numChannels
of these, one corresponding to each channel.
I begin by populating the window using the clear()
lambda function. The clear()
lambda iterates over each bit (i.e. sample) in each channel within the window, and sets the value of that bit to fabs(dist(gen)) > threshold
(that is, I generate a Gaussian-distributed sample on the fly, and set the bit equal to 1
if that sample exceeds the threshold, and 0
otherwise).
I then enter a for
loop which iterates samplingRate
times, simulating iterating over 1 second of data. With each iteration, I first count the number of channels in the window containing at least one sample exceeding the threshold by counting the number of channels containing at least one nonzero bit (using .any()
). If this number (numContributingChannels
) exceeds the channel threshold, I add one to the trigger count, step forward in i
by the writeDelay
(simulating the pause in data collection), and repopulate the window using the clear()
lambda. If this number does not exceed the channel threshold, I replace the "oldest" samples with new data using:
for(auto &channel: window)
channel.set(imod, fabs(dist(gen)) > threshold);
The index of the "oldest" sample is given by imod
(or i % windowSize
), which is calculated with each iteration. Replacing the oldest sample is equivalent to appending the new samples (next iterating in the time series) to the end of the window, and moving forward in time by one sample.
After iterating through 1 second worth of data, the number of triggers counted is returned. Here's an animation of a "toy model", where numChannels = 4
, windowSize = 6
, and channelThreshold = 3
.
Goals
As mentioned previously, the code runs much to slowly to be able to iterate over a useful amount of mock data. I'll need to significantly improve upon this in order for it to be useful for my purposes. I performed profiling with gprof
, and found that on the order of 30%
of execution time is spent on random number generation, and nearly 70%
on the outermost for
loop.
How can this be improved?