Fast Risk Battle Simulator

Question

I have written some code in python to determine the outcome(including the remaining troops of both sides) of a Risk battle fought until one side is defeated. When using the Numba's njit decorator the code can simulate one battle's outcome in ~ 0.2ms. This is not fast enough for the application I have in mind. Ideally it would be one to two orders of magnitude faster. Any advice on how to improve performance would be greatly appreciated.

Function Code:

@njit
def _attack_(n_atk,n_def):
    n_atk_dies = 0
    n_def_dies = 0
    while n_atk >0 and n_def>0:
        if n_atk>=3:
            n_atk_dies = 3
        if n_atk==2:
            n_atk_dies = 2
        if n_atk ==1:
            n_atk_dies = 1
        if n_def >=2:
            n_def_dies = 2
        if n_def ==1:
            n_def_dies = 1   
        atk_rolls = -np.sort(-np.random.randint(1,7,n_atk_dies))
        def_rolls = -np.sort(-np.random.randint(1,7,n_def_dies))
        comp_ind = min(n_def_dies,n_atk_dies)
        comp = atk_rolls[:comp_ind]-def_rolls[:comp_ind]
        atk_losses =len([1 for i in comp if i<=0])
        def_losses = len([1 for i in comp if i>0])
        n_atk -= atk_losses
        n_def -= def_losses
    return [n_atk,n_def]

Benchmarking Code:

atk_list = list(np.random.randint(1,200,10000))
def_list = list(np.random.randint(1,200,10000))
cProfile.run('list(map(_attack_,atk_list,def_list))',sort='tottime')

Answer 1

You play with up to 199 dice per player. So most of the time, the attacker will attack with 3 dice and the defender will attack with 2 dice. As long as that's the case, you could simply pick one of the three possible outcomes (with precalculated probabilities) and apply it. Outline:

def _attack_(n_atk,n_def):
    while n_atk >= 3 and n_def >= 2:
        # pick random outcome and update the dice numbers
    while n_atk >0 and n_def>0:
        # same as before
    return [n_atk,n_def]

You could even precompute the result distributions for all \$199^2\$ possible inputs and randomly pick the result according to the given input's distribution. Outline:

def _attack_(n_atk,n_def):
    return random.choices(results[n_atk, n_def], cum_weights=cum_weights[n_atk, n_def])

Answer 2

I haven't benchmarked any of the following but I suspect that they can help you out. I'd check each alteration separately.

First, selecting the number of dice (die is singular, dice is plural) as you have is confusing and probably slow.

n_atk_dice = min(n_atk, 3)
n_def_dice = min(n_def, 2)

Second, using np.sort is (a) overkill for the problem and (b) probably slow as you need to cross into another library to do it. What do I mean by "overkill"? General sorting algorithms have to work for any length of list. But there are optimal sorting algorithms for any fixed length of list. Pulling up my Art of Computer Programming Volume 3, we can write the following:

def sort2(l):
    if l[0] < l[1]:
        return l[(1, 0)]
    return l

def sort3(l):
    if l[0] < l[1]:
        if l[1] < l[2]:
            return l[(2, 1, 0)]
        elif l[0] < l[2]:
            return l[(1, 2, 0)]
        else:
            return l[(1, 0, 2)]
    elif l[1] < l[2]:
        if l[0] < l[2]:
            return l[(2, 0, 1)]
        else:
            return l[(0, 2, 1)]
    else:
        return l

sort = [None, lambda l: l, sort2, sort3]
...
    atk_rolls = sort[n_atk_dice](np.random.randint(1, 7, n_atk_dice))
    def_rolls = sort[n_def_dice](np.random.randint(1, 7, n_def_dice))

Third, use numpy to count instead of creating a new list to count:

def_losses = (comp > 0).count_nonzero()

and only iterate over the list once:

atk_losses = comp_ind - def_losses

Fourth, consider doing a similar trick for comparisons as you did for sorting. Instead of using the fully general numpy comparisons, array subtraction, etc, realize that there are two cases: comp_ind is either 1 or 2. If you write specific code for each case you may get faster code and possibly avoid using numpy altogether.

Finally, investigate random number generation. This may actually be your chokepoint. See e.g. https://eli.thegreenplace.net/2018/slow-and-fast-methods-for-generating-random-integers-in-python/. If you can get an all-base-Python solution (no numpy) with fast random numbers, or generate large batches of random numbers at once, you may get a performance boost.

Answer 3

Preprocessing

Be sure to read superb_rain's excellent answer.

The following code is meant as a preprocessing. It simply iterates over every single attack & defense configuration, and calculates the cumulative probability that the attacker manages to kill 0, 1 or 2 defense soldiers.

It is slow, but it's not a problem since it's only supposed to run once and output a probability distribution that you can use in your code.

from itertools import product
from collections import Counter

D6 = [1, 2, 3, 4, 5, 6]

probabilities = {}

for n_atk_dice in range(1, 4):
    for n_def_dice in range(1, 3):
        attacker_kills = Counter()
        for attack in product(D6, repeat=n_atk_dice):
            attack = sorted(attack, reverse=True)
            for defense in product(D6, repeat=n_def_dice):
                defense = sorted(defense, reverse=True)
                result = sum(1 for a,d in zip(attack, defense) if a > d)
                attacker_kills[result] +=1
        probabilities[(n_atk_dice, n_def_dice)] = [
            attacker_kills[0],
            attacker_kills[0] + attacker_kills[1],
            attacker_kills[0] + attacker_kills[1] + attacker_kills[2]
            ]

import pprint
pprint.pprint(probabilities)

probabilities is now:

{(1, 1): [21, 36, 36],
 (1, 2): [161, 216, 216],
 (2, 1): [91, 216, 216],
 (2, 2): [581, 1001, 1296],
 (3, 1): [441, 1296, 1296],
 (3, 2): [2275, 4886, 7776]}

That's all the information you need, and you won't need to sort any array anymore. You could save the code in a separate script, and simply use this literal definition in your code. You can check the values with other sources (e.g. https://web.stanford.edu/~guertin/risk.notes.html).

Use

With random.choices

As an example, for 3 vs 2:

import random
random.choices([0, 1, 2], cum_weights=probabilities[(3, 2)])
# => Either [0], [1] or [2]

There will be a 2275/7776 chance that the attack loses 2 soldiers, a (4886 - 2275)/7776 chance that both sides lose 1 soldier, and a (7776-4886)/7776 chance that the defense loses 2 soldiers.

For 1 vs 2:

random.choices([0, 1, 2], cum_weights=probabilities[(1, 2)])
# => [0] or [1]

Since the cumulative weights are [161, 216, 216], there's 0% chance that the defense loses 2 soldiers. Either the attack loses a soldier (with a 161/216 probability) or the defense loses one (with a 55/216 probability).

With random.random()

You could also define probabilities this way:

total = sum(attacker_kills.values())
probabilities[(n_atk_dice, n_def_dice)] = [
    attacker_kills[0] / total,
    (attacker_kills[0] + attacker_kills[1]) / total
]

The output becomes:

{(1, 1): [0.5833333333333334, 1.0],
 (1, 2): [0.7453703703703703, 1.0],
 (2, 1): [0.4212962962962963, 1.0],
 (2, 2): [0.44830246913580246, 0.7723765432098766],
 (3, 1): [0.3402777777777778, 1.0],
 (3, 2): [0.2925668724279835, 0.628343621399177]}

You can then simply get a random number between 0.0 and 0.99999999999 (with random.random()), and compare it to the two values.

• If the random number is smaller than the first number : defense loses no soldier.
• If the random number is between both numbers, defense loses 1 soldier.
• If the random number is larger than the second number : defense loses 2 soldiers.

If defense loses x soldiers, the attacker loses min(n_atk_dice, n_def_dice) - x soldiers.

Markov chains, part 1

Just for fun, I tried to use the non-cumulated probabilities to define a Markov chain:

# pip install PyDTMC
# see https://pypi.org/project/PyDTMC/
from pydtmc import MarkovChain, plot_graph

N = 30

probabilities = {
 (1, 1): [0.5833333333333334, 0.4166666666666667, 0.0],
 (1, 2): [0.7453703703703703, 0.25462962962962965, 0.0],
 (2, 1): [0.4212962962962963, 0.5787037037037037, 0.0],
 (2, 2): [0.44830246913580246, 0.32407407407407407, 0.22762345679012347],
 (3, 1): [0.3402777777777778, 0.6597222222222222, 0.0],
 (3, 2): [0.2925668724279835, 0.3357767489711934, 0.37165637860082307]}

states = [(a, b) for a in range(N) for b in range(N)]

#NOTE: A sparse matrix would be a good idea for large N
p = [[0 for _ in range(N * N + 2)] for _ in range(N * N + 2)]

p[-1][-1] = 1.0 # Absorbing state 'defense wins'
p[-2][-2] = 1.0 # Absorbing state: 'attack_wins'

def state_id(a, b):
    return a * N + b

for a, b in states:
    i = state_id(a, b)
    if a == 0:
        p[i][-1] = 1.0 # defense wins
    elif b == 0:
        p[i][-2] = 1.0 # attack wins
    else:
        a_dice = min(a, 3)
        b_dice = min(b, 2)
        d = min(a_dice, b_dice)
        p0, p1, p2 = probabilities[(a_dice, b_dice)]
        p[i][state_id(a-d,b)] = p0
        p[i][state_id(a-(d-1),b-1)] = p1
        if d == 2:
            p[i][state_id(a-(d-2),b-2)] = p2

state_names = [f'{a} vs {b}' for a, b in states] + ['attack_wins', 'defense_wins']
    
mc = MarkovChain(p, state_names)

The resulting graphs are interesting:

import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [15, 10]
plot_graph(mc, nodes_type=False, dpi=200)

mc.absorption_probabilities is a matrix containing the probabilities to win or lose depending on the initial state.

You can get the probability that the attacker wins for a given start (e.g. 10 vs 7), without any loop:

mc.absorption_probabilities[0][state_id(10, 7)]
# 0.7998329909591375

And by iterating over every initial state, it's possible to recreate the matrix presented in this answer:

import numpy as np
prob_matrix = np.array([[mc.absorption_probabilities[0][state_id(b,a)]
                         for a in range(N)] for b in range(N)])

It's also possible to see the intermediate steps during a random battle:

mc.walk(20, initial_state='12 vs 9')

It outputs:

['12 vs 7',
 '11 vs 6',
 '10 vs 5',
 '10 vs 3',
 '9 vs 2',
 '8 vs 1',
 '8 vs 0',
 'attack_wins', ....

See this answer for more examples.

Answer 4

Markov Chains, Part 2

Goal

My above answer was getting too long. The goal here will be to use a Markov chain to calculate the probability of reaching an absorbing state (end state, e.g. 3 vs 0 -> attacker wins or 0 vs 2 -> defense wins) from a transient state (initial state, e.g. 12 vs 7).

Code

Here's a slightly modified code:

# pip install PyDTMC
# see https://pypi.org/project/PyDTMC/
from pydtmc import MarkovChain

N = 30

probabilities = {
 (1, 1): [0.5833333333333334, 0.4166666666666667, 0.0],
 (1, 2): [0.7453703703703703, 0.25462962962962965, 0.0],
 (2, 1): [0.4212962962962963, 0.5787037037037037, 0.0],
 (2, 2): [0.44830246913580246, 0.32407407407407407, 0.22762345679012347],
 (3, 1): [0.3402777777777778, 0.6597222222222222, 0.0],
 (3, 2): [0.2925668724279835, 0.3357767489711934, 0.37165637860082307]}

states = [(a, d) for a in range(N) for d in range(N)]

#NOTE: A sparse matrix would be a good idea for large N
p = [[0 for _ in range(N * N)] for _ in range(N * N)]

def state_id(a, d):
    return a * N + d

def transient_state_id(a, d):
    return (a - 1) * (N - 1) + (d - 1)

for a, d in states:
    i = state_id(a, d)
    if a == 0 or d == 0:
        p[i][i] = 1.0 # someone wins -> absorbing state
    else:
        a_dice = min(a, 3)
        d_dice = min(d, 2)
        m = min(a_dice, d_dice) # how many dice are compared
        p0, p1, p2 = probabilities[(a_dice, d_dice)]
        p[i][state_id(a - m, d)] = p0
        p[i][state_id(a - m + 1, d - 1)] = p1
        if m == 2:
            p[i][state_id(a - m + 2, d - 2)] = p2

state_names = [f'{a}vs{d}' for a, d in states]
    
mc = MarkovChain(p, state_names)

Probabilities

mc.absorption_probabilities is the matrix with the desired information. It has 2 * N - 1 rows (1 for each end state) and (N - 1)**2 columns (1 for each initial state).

For example, with N=30, the matrix is (59, 841). For N=199, the shape would be (397, 39204). I'm not sure how long it would take to generate this matrix, though.

The matrix could be saved to a file once it is generated:

import numpy as np
with open(f'risk_markov_{N}.npy', 'wb') as f:
    np.save(f,  mc.transient_states)
    np.save(f, mc.absorbing_states)
    np.save(f, mc.absorption_probabilities)

# with open(f'risk_markov_{N}.npy', 'rb') as f:
    # print(np.load(f))
    # print(np.load(f))
    # print(np.load(f))

Probability distribution

Here's the code to get the probability distribution of 14 vs 5:

ps = mc.absorption_probabilities[:,transient_state_id(14, 5)]
# array([0.        , 0.04169685, 0.08077598, 0.07195455, 0.07046567,
#  0.05591131, 0.04746686, 0.03114066, 0.02060625, 0.00910449,
#  0.00328452, 0.        , 0. ....

ps is an array with 2 * N - 1 elements.

It's possible to select states with more than 5% probability:

['%s : %.1f %%' % (s, p*100) for s,p in zip(mc.absorbing_states, ps) if p > 0.05]

It returns:

['7vs0 : 6.1 %',
 '8vs0 : 7.6 %',
 '9vs0 : 10.5 %',
 '10vs0 : 12.2 %',
 '11vs0 : 14.5 %',
 '12vs0 : 14.7 %',
 '13vs0 : 12.4 %',
 '14vs0 : 9'.1 %']

Diagrams

# 3 vs 5 -> -2
def get_diff(vs):
    a, d = vs.split('vs')
    return int(a.strip()) - int(d.strip()) 

import pandas as pd
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = [15, 10]

attack, defense = 12, 7
ps = mc.absorption_probabilities[:,transient_state_id(attack, defense)]
df = pd.DataFrame(ps * 100, mc.absorbing_states, columns=['p'])
df['difference'] = df.index.map(get_diff)
df = df.sort_values(by='difference', ascending=False)
df = df[df.p > 0]
df['color'] = df.difference.map(lambda d: 'green' if d > 0 else 'red')
win = df[df.difference > 0].p.sum()
df.p.plot(kind='bar',
          color=df.color,
          title = 'Initial state : %d vs %d (%.1f %% win)' % (attack, defense, win))
plt.ylabel('Probability [%]')
plt.xlabel('End state')
plt.show()

It displays, for 12 vs 7:

or for 9 vs 14:

Just for fun, it's also possible to plot them all in a large plot:

import numpy as np
import matplotlib.pyplot as plt

n = 10
fig, subplots = plt.subplots(n, n, sharex='col', sharey='row')
for (i, j), subplot in np.ndenumerate(subplots):
    attack, defense = i + 1, j + 1
    ps = mc.absorption_probabilities[:,transient_state_id(attack, defense)]
    df = pd.DataFrame(ps * 100, mc.absorbing_states, columns=['p'])
    df['difference'] = df.index.map(get_diff)
    df = df.sort_values(by='difference', ascending=False)
    df = df[df.p > 0]
    df['color'] = df.difference.map(lambda d: 'green' if d > 0 else 'red')
    win = df[df.difference > 0].p.sum()
    df.p.plot(kind='bar', color=df.color,
              ax=subplots[i,j]
             )

fig.suptitle("Risk: Attack vs Defense")
plt.show()

Generating the matrix for N=100 took half an hour, required ~10GB of RAM and wrote a 16 MB npy binary file.

The script crashed for N=199 on my laptop because no memory was left. I guess that ~64GB should be enough.

Answer 5

If you profile the code without using njit, you will find that most of the time is consumed in calls to randint and sort:

atk_rolls = sort[n_atk_dice](np.random.randint(1, 7, n_atk_dice))
def_rolls = sort[n_def_dice](np.random.randint(1, 7, n_def_dice))

So, instead of generating and sorting those random numbers for each loop iteration, you could pregenerate them before: saving them on an array with length 200. (For the worst case, you will need 200 rows.)