<jupyter_start><jupyter_text><NAME>2012, 2020 The Unfinished Game ... of Risk![](components_.jpg)[<NAME>](https://web.stanford.edu/~kdevlin/)'s [book](https://www.amazon.com/Unfinished-Game-Pascal-Fermat-Seventeenth-Century/dp/0465018963) [*The Unfinished Game*](https://wordplay.blogs.nytimes.com/2015/12/14/devlin-unfinished-game/) describes how Fermat and Pascal discovered the rules of probability that guide gambling in games. The question they confront is: what if a gambling game is interrupted, but one player is in the lead by a certain score. How much of the pot should the leader get?My friends and I faced a similar question when a game of [Risk](https://www.ultraboardgames.com/risk/game-rules.php) ran on too long (as they often do) and we were unable to finish. Player **A** had just cashed in cards and added troops to a single large force in Brazil that was poised to make a sweeping attack on player **B**, whose territories were situated in Africa and Asia in such a way that **A** could attack North Africa first and then go from one territory to the next without ever having to branch off.  We wrote down the number of **A**'s armies in Brazil, **72**, and the number of armies in **B**'s successive territories: **22, 8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1.** What is the probability that **A** can capture all these territories?______ Monte Carlo SimulationWe can answer the question with a [**Monte Carlo simulation**](https://en.wikipedia.org/wiki/Monte_Carlo_method) that follows the rules of Risk and uses random numbers to roll the dice. First some preliminaries:<jupyter_code>from typing import List, Tuple
import matplotlib.pyplot as plt
import numpy as np
import random

attackers   = 72
territories = (22, 8, 2, 2, 2, 7, 1, 1, 3, 1, 2, 3, 5, 1)

die = (1, 2, 3, 4, 5, 6)
Die = int # a type alias

random.seed(42) # For reproducability<jupyter_output><empty_output><jupyter_text>How many territories and how many total defenders?<jupyter_code>len(territories), sum(territories)<jupyter_output><empty_output><jupyter_text>It looks fairly even: the attackers have 12 more armies, 72 to 60, but they will have to leave an army in each of the 14 territories along the way.  Battless, Invasions, and CampaignsIn Risk a **battle** consists of a one-time roll of some dice and the resulting loss of armies. The attacker will roll three dice if possible (but can roll a maximum of the number of armies in the attacking territory minus one) and the defender will roll two dice if possible (or only one if they have only one army remaining). We compare the highest die on each side, with the defender losing an army if the attacker's die is higher, and the attacker losing an army if tied or lower. Then if both sides rolled at least two dice, we do the same comparison with the second highest die on each side. The function `deaths` returns a tuple of (number of attacking armies that die, number of defending armies that die).<jupyter_code>def deaths(attack_dice: List[Die], defend_dice: List[Die]) -> Tuple[int, int]:
    """The number of (attackers, defenders) who die due to this roll."""
    dead = [0, 0]
    for a, d in zip(sorted(attack_dice, reverse=True), 
                    sorted(defend_dice, reverse=True)):
        dead[a > d] += 1
    return tuple(dead)
def test1() -> bool:
    # Four examples from www.ultraboardgames.com/risk/game-rules.php
    assert deaths([6, 1, 1], [3])    == (0, 1) # 6 beats 3, so defender loses 1
    assert deaths([6, 2, 1], [3, 2]) == (1, 1) # 6 beats 3, but 2 ties 2, so 1 loss each
    assert deaths([3, 3],    [4, 3]) == (2, 0) # 3 loses to 4, and attacker 3 loses to 3
    assert deaths([6],       [5, 4]) == (0, 1) # 6 beats 5, so defender loses 1
    return True
    
test1()<jupyter_output><empty_output><jupyter_text>An **invasion** consists of a series of battles until either the defenders are all defeated (in which case the attacker can move armies into the captured territory, but must leave one army behind), or the attackers have fewer than two armies left, and can no longer attack. A **campaign** consists of a sequence of invasions, designed to capture all the defender's territories.In the function `campaign`, the variable `attackers` tracks the number of attackers invading the current territory. (The total number of armies for the attacking side will in general be more than `attackers`, because there will be one army left behind in each captured territory.) If `verbose` is true, we'll `say` what is happening along the way. In the end, return the number of `attackers` in the final territory minus the total number of defenders still alive. This will be positive if all the defenders have been defeated, and nonpositive if not.<jupyter_code>def campaign(attackers: int, territories: List[int], verbose=True) -> int:
    """Given a number of armies for attacker and a list of territory army sizes for defender,
    randomly play to the end, and return remaining attackers minus remaining defenders."""
    for t, defenders in enumerate(territories, 1):
        if verbose: say(attackers, 'attack', defenders, t)
        while attackers > 1 and defenders > 0:
            A, D = deaths(roll(min(3, attackers - 1)), 
                          roll(min(2, defenders)))
            attackers -= A
            defenders -= D
        if attackers == 1:
            if verbose: say(attackers, "can't beat", defenders, t)
            break
        attackers -= 1 # Capture a new territory, leave one behind
    return attackers - (defenders + sum(territories[t:]))

def say(attackers, action, defenders, t):
    """Say what is happening."""
    print(f'{attackers:3} armies {action} {defenders:2} defenders in territory {t}')
        
def roll(n) -> List[Die]: 
    """A random roll of `n` dice, sorted largest first."""
    return sorted((random.choice(die) for _ in range(n)), reverse=True)<jupyter_output><empty_output><jupyter_text>The AnswerLet's answer the question of who wins the unfinished game:<jupyter_code>campaign(attackers, territories)<jupyter_output>72 armies attack 22 defenders in territory 1
 48 armies attack  8 defenders in territory 2
 44 armies attack  2 defenders in territory 3
 40 armies attack  2 defenders in territory 4
 35 armies attack  2 defenders in territory 5
 33 armies attack  7 defenders in territory 6
 23 armies attack  1 defenders in territory 7
 22 armies attack  1 defenders in territory 8
 20 armies attack  3 defenders in territory 9
 19 armies attack  1 defenders in territory 10
 17 armies attack  2 defenders in territory 11
 16 armies attack  3 defenders in territory 12
 15 armies attack  5 defenders in territory 13
 13 armies attack  1 defenders in territory 14<jupyter_text>The attackers won, capturing all the territories and moving 12 remaining armies into the final territory.But that was just one simulation; other simulations could have different results. Let's summarize, say, 100,000 simulations, see how long it takes, and plot a histogram of the resulting scores:<jupyter_code>%time scores = [campaign(attackers, territories, False) for _ in range(100000)]
def summary(scores):
    """Summarize the scores from simulations."""
    P = np.average([s > 0 for s in scores])    
    N, avg = len(scores), np.average(scores)
    title = f'Attacker wins {P:.0%}. Average armies left: {avg:.1f}'
    plt.ylabel(f'Frequency (out of {N:,d} runs)'); 
    plt.xlabel('Armies left (positive for attacker)')
    plt.hist(scores, bins=max(scores) - min(scores) + 1)
    plt.plot([0, 0], [0, N/25], 'r:'); plt.title(title)
    
summary(scores)<jupyter_output><empty_output><jupyter_text>We see that the attackers win 81% of the time, and the scores look roughly like a bell-shaped curve, but with a non-normal pattern on the left side. What's causing the non-normal pattern? Note that the number of defenders in the final four territories are `(2, 3, 5, 1)`, and the width of the four spikes on the left are `(1, 2, 4, 1)`, one less than the number of defenders (except for the last spike). I think there are spikes rather than a smooth curve because it is doubly-difficult to capture a territory: you need to have two armies left, not just one, so that you can leave one behind.   Exact ProbabilitiesBy repeatedly running a simulation, we can get approximate probabilities. But what if we wanted **exact probabilities** (at least exact under the assumption that the dice are exactly fair and that floating point roundoff is not an issue)? I'll start by defining the function `outcomes`, so that, for example, `outcomes(3, 2)` returns a Counter of all possible outcomes of rolling 3 attacker dice versus 2 defender dice; each outcome is a pair of `(attacker_deaths, defender_deaths)` and has a count of how often it occurred.<jupyter_code>from functools   import lru_cache
from collections import Counter
import itertools

@lru_cache()
def outcomes(num_A: int, num_D: int) -> List[Tuple[int, int]]:
    """All (equiprobable) outcomes of (attacker_deaths, defender_deaths)
    from rolling num_A dice versus num_D dice."""
    return Counter(deaths(dice[:num_A], dice[num_A:])
                   for dice in itertools.product(die, repeat=num_A + num_D))
outcomes(3, 2)<jupyter_output><empty_output><jupyter_text>The result for `outcomes(3, 2)` says that the attackers lose two armies in 2,275 out of the $6^5 = 7,776$ different outcomes; the defenders lose two in 2,890 outcomes; and they each lose one in the remaining 2,611 outcomes. I found a [web page](http://datagenetics.com/blog/november22011/) that lists the results for all six possible battles (attackers roll 1, 2, or 3 dice; defenders roll 1 or 2); let's verify that our `outcomes` agree with theirs:<jupyter_code>def test2() -> bool:
    # See http://datagenetics.com/blog/november22011/
    assert outcomes(1, 1) == {(1, 0):   21, (0, 1):   15}
    assert outcomes(2, 1) == {(1, 0):   91, (0, 1):  125}
    assert outcomes(3, 1) == {(1, 0):  441, (0, 1):  855}
    assert outcomes(1, 2) == {(1, 0):  161, (0, 1):   55}
    assert outcomes(2, 2) == {(2, 0):  581, (1, 1):  420, (0, 2):  295}
    assert outcomes(3, 2) == {(2, 0): 2275, (1, 1): 2611, (0, 2): 2890}
    return True
    
test2()<jupyter_output><empty_output><jupyter_text>Now we know exact outcomes for a single roll of the dice; what about an invasion where $A$ attackers keep attacking until they defeat $D$ defenders (or all but one attacker dies trying)? The function `winP(A, D)` gives the probability of the attackers winning. It works recursively. The two base cases are that the probability is zero if the attackers don't have at least two armies and the probability is one if the defenders have no armies. In the recursive case we observe the possible outcomes for the first battle, and then compute the win probability for the invasion as the average, over every possible outcome of the dice rolls, of the win probability for the number of remaining attackers and defenders after the battle, weighted by the number of times that the outcome occurs. (We use `np.average` because it accepts an optional `weights` argument.)<jupyter_code>@lru_cache(None)
def winP(attackers: int, defenders: int) -> float:
    """The probability that `attackers` can invade and defeat all the `defenders`."""
    if attackers <= 1:
        return 0
    elif defenders == 0:
        return 1
    else:
        battle = outcomes(min(3, attackers - 1), min(2, defenders))
        return np.average([winP(attackers - A, defenders - D) for (A, D) in battle],
                          weights=list(battle.values()))<jupyter_output><empty_output><jupyter_text>Let's try a simple example: `winP(2, 1)` should be `15/36`, the same as the probability of the attacker prevailing in `outcomes(2, 1)`:<jupyter_code>winP(2, 1) == 15/36<jupyter_output><empty_output><jupyter_text>What's the probability that 12 attackers successfully invade 10 defenders?<jupyter_code>winP(12, 10)<jupyter_output><empty_output><jupyter_text>Let's make a chart, with the number of defenders varying from 1 to 60, and the number of attackers separated into eight cases (depicted as eight lines), where in each case there are Δ more attackers than defenders, for Δ = -5, -2, 0, 1, 2, 5, and 10:<jupyter_code>def chart(Ds=range(1, 61), deltas=(-5, -2, -1, 0, 1, 2, 5, 10)):
    plt.figure(figsize=(9, 6)); plt.grid()
    plt.title('Each line: attackers with Δ more armies than defenders')
    plt.xlabel('Number of Defenders'); plt.ylabel('Win Probability for Attackers')
    for delta in reversed(deltas):
        Ps = [winP(max(0, D + delta), D) for D in Ds]
        plt.plot(Ds, Ps, '.-', label=f'Δ={delta}')
    plt.legend()
chart()<jupyter_output><empty_output><jupyter_text>Note that the purple line (fourth from bottom), where the number of attackers exactly equals the number of defenders, gives a low win probability for a small attacking force, but reaches 50% for 12-on-12, and 73% for 60-on-60. The red line, where the attackers have one more army than the defenders, dips from one to two defenders but is over 50% for a 6-on-5 attack. Similarly, the green line, where the attackers have a surplus of two armies, dips sharply from 75% to 66% as the number of defenders goes from 1 to 2, dips slightly more for 3 and 4 defenders, and then starts to rise. So overall, an attacker does not need a big advantage in armies as long as there are many armies on both sides. Even when the attacker is at a disadvantage in numbers (as in the bottom grey line where the attacker has five fewer armies), the attacker can still have an advantage in win probability; `winP(55, 60)` is about 57%.  Simulation versus Exact ComputationLet's see how the exact computation compares with a simulation:<jupyter_code>A, D = 32, 30

exact = winP(A, D)
simul = np.average([campaign(A, [D], False) > 0 for _ in range(10000)])

exact, simul<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09366887799250072, 'fraction_numerical': 0.032016152293048744, 'mean_word_length': 4.157679434734102, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 66, 'https://': 5, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '138906', 'license': 'MIT', 'n_tokens_mistral': 4521, 'n_tokens_neox': 3917, 'n_words': 1922, 'path': 'ipynb/risk.ipynb', 'repo_name': 'filmy33/pytudes', 'size': '13880'}