Risk Simulation Using Python

RiskSimulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Risk Simulations\n",
    "\n",
    "The following are Python simulations based on the Java simulation written by Thomas Glezen.  The first simulation uses base Python lists.  The second uses matrices and vectors provided by the [NumPy](https://docs.scipy.org/doc/numpy/index.html).\n",
    "\n",
    "## Base Python\n",
    "\n",
    "This version uses a loop in which *listcomps* are used to generate the die rolls for each attack."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "numOfAttacks = 1000000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Attacker wins 1,106,639 times.\n",
      "Defender wins 1,893,361 times.\n",
      "Defender/Attacker = 1.71\n"
     ]
    }
   ],
   "source": [
    "import random as r\n",
    "\n",
    "aWins = dWins = 0\n",
    "for i in range(numOfAttacks):\n",
    "    a = sorted([r.randint(1,6) for i in range(3)])\n",
    "    d = sorted([r.randint(1,6) for i in range(3)])\n",
    "    aWins += sum([a[i] >  d[i] for i in range(3)])\n",
    "    dWins += sum([a[i] <= d[i] for i in range(3)])\n",
    "\n",
    "print(\"Attacker wins {:,d} times.\".format(aWins))\n",
    "print(\"Defender wins {:,d} times.\".format(dWins))\n",
    "\n",
    "fraction = 1.0 * dWins / aWins;\n",
    "print(\"Defender/Attacker = {:.2f}\".format(fraction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NumPy\n",
    "\n",
    "The **NumPy** library sacrifices flexibility for speed.  It also supports vector and matrix operations.  Check it out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The NumPy integer random number generator takes a lower range, upper range, and a pair of dimensions for the matrix that will hold the random numbers.  In our case, we create a row for each attack and a column for each die roll (3 for attacker, 3 for defender)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "diceRolls = np.random.randint(1, 7, (numOfAttacks, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `np.sort` function sorts an array along an axis (either row axis (0) or column axis (1)).  We use the first three columns for the attacker and the last three for the defender."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "aSortedRolls = np.sort(diceRolls[:,0:3], axis=1)\n",
    "dSortedRolls = np.sort(diceRolls[:,3:6], axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With all the rolls matched up, we can compare them element-wise.  The expression `aSortedRolls > dSortedRolls` is a matrix of booleans (True/False).  When you sum a matrix of booleans, you add one for each `True` and zero for each `False`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Attacker wins 1,106,943 times.\n",
      "Defender wins 1,893,057 times.\n",
      "Defender/Attacker = 1.71\n"
     ]
    }
   ],
   "source": [
    "aWins = np.sum(aSortedRolls >  dSortedRolls)\n",
    "dWins = np.sum(aSortedRolls <= dSortedRolls)\n",
    "\n",
    "print(\"Attacker wins {:,d} times.\".format(aWins))\n",
    "print(\"Defender wins {:,d} times.\".format(dWins))\n",
    "\n",
    "fraction = 1.0 * dWins / aWins;\n",
    "print(\"Defender/Attacker = {:.2f}\".format(fraction))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Performance\n",
    "\n",
    "The NumPy library requires one to think in terms of matrices and vectors rather than loops.  The NumPy libraries are backed by **C** libraries that are not flexible with respect to types, but are highly performant.\n",
    "\n",
    "The **listcomp** version uses little memory.  It reuses the same memory for each iteration of the loop.  It also doesn't require NumPy which is not always easy to install.  But it takes over 3 1/2 minutes to run 10,000,000 iterations.\n",
    "\n",
    "The **NumPy** version executes 10,000,000 simulations in under 3 seconds!  But it doesn't use memory efficiently.  The rolls for all simulations are stored at one time.  This is fine if you have enough memory."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python (p3)",
   "language": "python",
   "name": "p3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}