Improving execution time of physics, data acquisition triggering system simulation

Question

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?

Answer 1

Use math to simplify the problem

You should be able to come up with a closed-form solution for the problem by using probability theory. This way you won't need to run a simulation at all, and you will be able to get the answer in microseconds.

But first let's assume you still need a simulation. You can still use math to reduce the number of steps your code needs to get an answer. For example, instead of simulating each individual timestep, most of which nothing interesting happens (either there is no trigger, or you are in the window), I would try to calculate the time to the next event instead.

First, given a fixed threshold, you don't need to sample the normal distribution, you just need to know the probability that a random draw would result in a value higher than threshold. You can then, using the binomial distribution, calculate the probability that at least channelThreshold channels would get a value higher than threshold. That will be the chance to trigger at each time step. So you could then use a std::bernouilli_distribution to ask each time step "did I trigger?".

But it gets better than that. If you have a fixed probability to trigger an event each time step, then it is approximately a Poisson process, and the time between triggers follow an exponential distribution. So instead of simulating each time step, just draw a random value from the inverse CDF of the exponential distribution, and that's how long you have to wait for the next trigger.

Once you trigger, you can just continue calculating times to the next trigger, and you just ignore them while you are still inside the window. If your window is very large compared to the time between triggers, then you could even optimize that by calculating the distribution function of "the time of the next trigger after a window".

So your code should look like:

int Antenna::triggerRate(double threshold, double temperature) {
    double lambda = /* calculate rate of raw triggers */;

    std::random_device rd;
    std::mt19937 gen(rd());
    std::uniform_real_distribution dist(0.0, 1.0);
    auto nextEvent() = [&](){ return -std::log(dist(gen)) / lambda; };

    int numTriggers = 0;
    int t = nextEvent(); // first trigger

    while (t < samplingRate) {
        ++numTriggers;

        // Find next trigger after the window
        auto windowEnd = t + windowSize;
        do {
            t += nextEvent();
        } while (t < windowEnd);
    }

    return numTriggers;
}

I'm leaving the calculation of the lambda parameter as an excercise for the reader.

How precise should your answer be?

You only simulate one second, so the answer you get will just be an approximation. How precise it is depends on how many triggers happen within one second. If only a few triggers happen, your answer will be very imprecise, if you have millions of triggers happening, your answer might be more precise than you need. Thus, instead of fixing your simulation to 1 second, determine first what an acceptable error margin is, then calculate how many triggers you approximately need to see to get an answer within that margin. This in itself might also avoid wasting CPU time.

Use std::uint64_t instead of int

You use int in your code, but if you really want to calculate one second at 40 GS/s, then you need at least 36 bits to count up to 40e9. Even on 64-bit systems, int is typically only 32 bits, and 1 of them is used for the sign of the value. Use std::uint64_t to make sure your counters can count far enough.

Closed-form solution

Once a trigger happens and your equipment starts recording, you are ignoring any other triggers for windowSize samples. Then you wait for the next trigger after the window ends. You might think it will be difficult to calculate that, because while you have a distribution for the time between events, how do you account for the time that is already spent between the last event inside the window and the end of the window? That's the neat part, you don't!

In your problem every sample there is a probability an event will trigger, and that probability is independent of what happened in the previous sample. So it doesn't matter which time point I choose: regardless of whether there was an event at that time or not, the exponential distribution will give me the time to the next event. So if you have calculated lambda, then the average time between triggers being recorded is just windowSize + 1 / lambda.