How to get a running average of array values at specific indexes over multiple iterations

Accelerated Computing CUDA CUDA Programming and Performance
1

I have a device array that stores scores as int32 values and another that stores flags as int8_t. The i-th score corresponds to the i-th game in my vectorized simulator and the i-th flag is nonzero if the i-th game needs to be restarted. The arrays are updated at each step of the simulation, and the terminated games are restarted.

How could I get an estimate of the average score of terminated games over the course of a simulation consisting 100k+ steps? Runtime performance matters more than “statistical correctness”.

In my prior host-only implementation I simply used an exponential moving average with the terminated scores, updated by a single worker thread at each simulation step.

A minimal version of my current kernel would be:

__global__ void reset_score(const int game_count, const int8_t *terminated_flags, int *scores)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= game_count)
        return;

    if (terminated_flags[idx] == 0)
        return;

    scores[idx] = 0;
    return;
}
2

To understand, what you could need: At each simulation step, all games are tested, whether they are terminated. If yes, then there scores are read and an average is created from the scores. Optionally combined with some of the previous simulation steps for a(n exponential) moving average.

Probably the actual math is less of a (performance or implementation) problem, but collecting all the scores.
There are two possibilities
a) As your scores and your termination flags are in device arrays, you could just use those. Disadvantage, if it is done by some worker thread/warp or several, you have to make sure that there is no conflict with concurrent access, as the game is restarted. Of course you could make sure that the next game and its termination flag is stored at a different index.
b) As soon as a game finishes, you collect the result with atomic (add) operations. First warp wide, then SM (=block) wide, then device wide. For each step, you add up the number of terminated games, and the added scores. Then you can use any mathematical formula to create out of it suitable averages. You can weight them by terminated games per iteration, etc. Advantage of b) is that you need less global memory accesses and have no problem with concurrent accesses.
Probably both approaches a) and b) can be made to run fast enough in the background.

3

[…] First warp wide, then SM (=block) wide, then device wide […]

This likely could work and I need to find such example first.

More details:

The simulator stores all data on the device: one array for the scores, one array for the turns, one for the boards and so on. (Only certain utility variables are stored on the host.) One thread is responsible for one game (as of now), and the games are independent of one another. The i-th game will always be at position i.

One simulator step includes progressing the games based on the given actions and restarting the ones that terminate. At the beginning of the next step, all the games are ready to go.

My goal would be to get an estimate at each simulation step of the average score of terminated games.

4

You could use a custom atomic.

__global__ void reset_score(const int game_count, const int8_t *terminated_flags, int *scores, unsigned long long *my_custom_atomic)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= game_count)
        return;

    if (terminated_flags[idx] == 0)
        return;
  // I assume a non-zero flag value here corresponds to a score to be averaged
  //  before it is zeroed
    int my_score = scores[idx];
    bool done = FALSE;
    while (!done){
        unsigned long long my_atomic_value = *my_custom_atomic;
        long long my_avg_score = (long long)(my_atomic_value >> 32);
        unsigned my_avg = (unsigned)(my_atomic_val & 0x0FFFFFFFFULL);
        my_avg_score *= my_avg;
        my_avg_score += my_score;
        my_avg_score /= (my_avg+1);
        my_avg_score <<= 32;
        unsigned long long my_new_atomic_val = my_avg_score+my_avg+1;
        unsigned long long old = atomicCAS(my_custom_atomic, my_atomic_val, my_new_atomic_val);
        done = (old == my_atomic_val);}
    scores[idx] = 0;
    return;
}

(coded in browser, not tested)

You would need a global unsigned long long location pointed to by my_custom_atomic, initialized to zero. The method to extract the average score from that global location is hopefully evident from the code above.

There is an overflow hazard if the number of terminated games exceeds ~4 Billion.

Integer truncation, amongst possibly other factors, could/will impact the numerical accuracy of this computed average. A somewhat better approach might be to use floating point data for my_avg_score.

Perhaps similar to what curefab had in mind with the description of method b.

5

I had something similar, but not fully the same in mind.

As far as I understood, the games are in a kind of lockstep: Each game does a step, then we restart the terminated games, then all games do the next step.

We have (create) two arrays over the maximum number of steps:
[Those I introduced, they are different arrays than the ones mentioned in the original post]

One array contains the number of terminated games in that step.
The other array contains the sum of the scores of the terminated games in that step.

The following is a reduction operation:
We could either fill the arrays by having each thread execute an atomicAdd operation on both arrays for the current step index:
→ The number of terminated games is increased by 1 (1 is added)
→ The sum of the scores is increased by the score

To simplify, one can first reduce over the warp (__reduce_add_sync), then in shared memory for the whole block (atomicAdd), then in global memory for the whole grid.

Or use the reduce function of cooperative groups (CUDA C++ Programming Guide (Legacy) — CUDA C++ Programming Guide)

The advantage is that there is only one global memory access neccessary per block (instead of one per thread).

Compared to Robert’s solution the average is not already calculated, but can afterwards be easily calculated from that arrays.
The arithmetic average per step is (sum of scores) / (number of terminated games) and several steps can also be combined (to get a moving average as mentioned as example in the original post).

The stored value overflows more easily, when we store the full sum over all terminated games in that step but either one can remove some bits to get an approximate result,
or - as in Robert’s solution - you can get back the result of the atomic addition. This result can be used to detect overflows (the value overflowed, if the result of the addition is smaller than what you added) and do wider additions than would be possible with a single atomicAdd.

6

Based on your inputs this is what I came up with:

__global__ void reset_score(const int game_count, const int8_t *terminated_flags, int *scores, float *avg_score)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx >= game_count)
        return;

    if (terminated_flags[idx] == 0)
        return;

    // This is the interesting part
    if (idx < 1024) { // Somewhat arbitrary
        int curr_score = scores[idx];
        bool done = false;
        while (!done) {
            float old_avg = *avg_score;
            float new_avg = 5e-4 * curr_score + (1.0 - 5e-4) * old_avg; // EMA alpha 2/(N+1), with N of 3999

            float pre_update =
                __int_as_float(atomicCAS((int *)avg_score, __float_as_int(old_avg), __float_as_int(new_avg)));
            done = (pre_update == old_avg);
        }
    }

    scores[idx] = 0;
    return;
}

Here I combine Robert’s solution with my old approach. Someone on SO pointed out how a subsample of the games could provide me with the estimate I want, which should also reduce the number of per-thread memory accesses.

My project did not seem to slow down from this approach. I hope I got the float atomicCAS correctly.

Thank you Robert for the code example, it helped a lot. Thank you Curefab for bringing cooperative groups to my attention and providing a link; that was something new.

7

I first wanted to recommend (as a small numerical improvement):

float new_avg = 5e-4 * curr_score + (1.0 - 5e-4) * old_avg;
to
float new_avg = 5e-4 * (curr_score - old_avg) + old_avg;

Mathematically it is the same, but you avoid, losing precision with the (1.0 - 5e-4)

but 0.9995 is quite representatable as float. But watch out, if you use smaller constants than 5e-4.

8

Thanks for the feedback!