# Fast Risk Battle Simulator

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')

• What application do you have in mind? – Stefan Pochmann Jan 17 at 17:15
• Reinforcement Learning environment. Don't want the game environment speed to be a limiting factor. – mcgowan_wmb Jan 17 at 17:19
• – superb rain Jan 17 at 18:15
• Will not do this again. I rolled back thinking you had not rolled back. – mcgowan_wmb Jan 17 at 18:25

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])

• As far as I can tell, "pick random outcome" for 3 vs 2 would be x = random.choices([0, 1, 2], cum_weights=[2275, 4886, 7776]). The defense loses x soldiers, the attackers loses (2-x). It should be many orders of magnitude faster than what OP wrote. – Eric Duminil Jan 18 at 16:04
• @EricDuminil Yeah, that looks about right. I computed the numbers as well, but already deleted the code. Had a version with just r = randrange(7776) and two boundary comparisons instead of choices. Did you get the numbers with math somehow, or did you evaluate all 7776 possibilities and counted (that's what I did)? I doubt it'll be "many orders of magnitude", but then again, I don't know how Numba is affected by all this. – superb rain Jan 18 at 17:32
• @superbrain Simple brute-force, because there aren't that many possibilities. I took the liberty to write an answer as a complement to yours. I might have been too optimistic, but I sure hope to see a good improvement since now, only a dict lookup and one single random number are needed, compared to random arrays generation and sorting. – Eric Duminil Jan 18 at 17:43
• @mcgowan_wmb: Good to know. Do you need every intermediate step e.g. ['(12, 7)', '(12, 5)', '(12, 3)', '(10, 3)', '(8, 3)', '(8, 1)', '(8, 0)', 'attack_wins'], or just to know if the attack wins? – Eric Duminil Jan 19 at 9:40
• @Eric Duminil the only important information is the end state eg (8,0) – mcgowan_wmb Jan 19 at 12:30

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.

• You'll also want to investigate whether views or copies are faster when sorting - l[(2,1, 0)] vs l[[2, 1, 0]] – ruds Jan 17 at 18:01
• Thank you! I will try all of these changes and report back. – mcgowan_wmb Jan 17 at 18:11

# Preprocessing

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.

• @superbrain: Thanks again for the bug finding. Are the last sentences correct, now? – Eric Duminil Jan 18 at 20:14
• Yes, looks good now. – superb rain Jan 18 at 20:22
• @Eric Duminil Could your Markov chain approach be used to compute the exact probabilities of each battle outcome, including remaining troops? – mcgowan_wmb Jan 19 at 22:14
• @mcgowan_wmb: Yes, it works. The absorbing states need to be changed (so no "attack wins" anymore, but simply "4 vs 0") and the code needs to be slightly modified. For N=50, it takes less than a minute and outputs a (99 * 2401) matrix, with the probability distribution of each end state, for each initial state. – Eric Duminil Jan 19 at 23:27
• @mcgowan_wmb For example, here are the possible outcomes of '12 vs 14' : ['0 vs 2 : 7.1 %', '0 vs 3 : 6.7 %', '0 vs 4 : 7.2 %', '0 vs 5 : 6.4 %', '0 vs 6 : 6.3 %', '0 vs 7 : 5.1 %', '3 vs 0 : 7.6 %', '4 vs 0 : 7.3 %', '5 vs 0 : 6.4 %', '6 vs 0 : 5.4 %'] I filtered the ones with > 5%. It's getting late here, I might update the code tomorrow. – Eric Duminil Jan 19 at 23:29

# 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:


## 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.

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.)

• I'm no sure I understand. n_atk_dice is 3 at most, and n_def_dice is 2 at most, right? Your sort example, which you apparently got from another answer, doesn't make sense with 200 dice, does it? – Eric Duminil Jan 18 at 15:28
• Also, you cannot sort the arrays more than 3 elements at a time. Let's say you have an army of 200, against one single defender. It makes a huge difference if you roll 199 times 1 followed by one 6 at the end, or if you roll a 6 directly at the beginning. – Eric Duminil Jan 18 at 15:57
• @EricDuminil: n_atk_dice would be an 3x200 array. And n_def_dicewould be an 2x200 array. If you have an army of 200 vs one: instead of executing a loop (with at most 200 iterations) you could pregenerate the numbers and have them on an array. So, inside the loop there will not be any call to randint. – user3808394 Jan 19 at 15:03
• The main idea is that it's a lot slower calling randint 200 times to generate a number than calling it one time to generate 200 numbers. – user3808394 Jan 19 at 15:06
• Did you test your idea, and benchmarked it? It also means that you'll sometimes generate 5*200 random numbers, for a single dice battle. – Eric Duminil Jan 19 at 15:18