#!/usr/bin/env python
# coding: utf-8

# <div align="right" style="text-align:right"><a href="http://norvig.com">Peter Norvig</a><br>March 2019</div>
# 
# #  Baseball Simulation
# 
# The [538 Riddler for March 22, 2019](https://fivethirtyeight.com/features/can-you-turn-americas-pastime-into-a-game-of-yahtzee/) asks us to simulate baseball using probabilities from a 19th century dice game called *Our National Ball Game*. The Riddler description of the rules said *you can assume some standard baseball things* but left many things unspecified, so I [looked up](http://baseballgames.dreamhosters.com/BbDiceHome.htm) the original rules of *Our National Ball Game*, which are shown below and which, it turns out, contradict some of the *standard baseball things* assumed by 538. I'll go with the rules as stated below.
# 
# 
# |RULES FOR PLAYING "OUR NATIONAL BALL GAME"|DICE ROLL OUTCOMES|
# |-----|-----|
# | ![](http://baseballgames.dreamhosters.com/bbOurNationalBallGameMcGillDelanyHoag86diE2.jpg) |  ![](http://baseballgames.dreamhosters.com/bbOurNationalBallGameMcGillDelanyHoag86ddE.jpg) |
# 
# *Note: There is a great history of dice baseball [here](https://baseballgames.dreamhosters.com/BbDiceHome.htm).*
# # Design Choices
# 
# 
# - Exactly one thing happens to each batter. I'll call that an **event**.
# - To clarify: the dice roll `1,1` has probability 1/36, whereas `1,2` has probability 2/36, because it also represents `2,1`.
# - The "One Strike" dice roll is not an event; it is only *part* of an event. From the probability of a "One Strike" dice roll, 7/36, I compute the probability of three strikes in a row, `(7/36)**3 == 0.00735`, and call that a strikeout event. 
# - I'll represent events with the following 11 one letter **event codes**:
#   - `1`, `2`, `3`, `4`: one-, two-, three-, and four-base (home run) hits. Runners advance same number of bases.
#   - `B`: base on balls. Runners advance only if forced.
#   - `D`: double play. Batter and runner nearest home are out; others advance one base.
#   - `E`: error. Batter reaches first and all runners advance one base.
#   - `F`, `K`, `O`: fly out, strikeout, foul out. Batter is out, runners do not advance.
#   - `S`: called "out at first" in rules, but actually a sacrifice. Batter is out, runners advance one base.
# 
# 
# #  Implementation

# In[1]:


import matplotlib.pyplot as plt
import random
from statistics import mean, stdev
from collections import Counter
from itertools import islice


# In[2]:


event_codes = {
    '1': 'single',        '2': 'double',      '3': 'triple',   '4': 'home run',
    'B': 'base on balls', 'D': 'double play', 'E': 'error',
    'F': 'fly out',       'K': 'strikeout',   'O': 'foul out', 'S': 'out at first'}


# I'll define the function `inning` to simulate a half inning and return the number of runs scored.  Design choices for `inning`:
# 
# - I'll keep track of runners with a set of occupied bases; `runners = {1, 3}` means runners on first and third.
# - I'll keep track of the number of `runs` and `outs` in an inning, and return the number of `runs` when there are three `outs`.
# - Each event follows four steps. If `runners = {1, 3}` and the event is `'2'` (a double), then the steps are:
#   - The batter steps up to the plate. The plate is represented as base `0`, so now `runners = {0, 1, 3}`.
#   - Check if the event causes runner(s) to be out, and if the inning is over. In this case, no.
#   - Advance each runner according to `advance(r, e)`. In this case, `runners = {2, 3, 5}`.
#   - Remove the runners who have `scored` and increment `runs` accordingly. In this case, runner `5` has scored, so we increment `runs` by 1 and end up with `runners = {2, 3}`.
# - I want `inning` to be easily **testable**: I want to say `assert 2 = inning('1KO4F')`.
# - I also want `inning` to be capable of simulating many independent random innings.  So the interface is to accept an *iterable* of event codes. That could be string, or a generator, as provided by `event_stream()`.
# - I want `inning` to be **loggable**: calling `inning(events, verbose=True)` should produce printed output for each event.
# - `advance(r, e)` says that a runner advances `e` bases on an `e` base hit; one base on an error, sacrifice, or double play; and one base on a base on balls only if forced.
# - A runner on base `r` is `forced` if all the lower-numbered bases have runners.
# - `ONBG` is defined as a generator of random events with the probabilities from "Our National Ball Game".

# In[3]:


def inning(events, verbose=True) -> int:
    """Simulate a half inning based on events, and return number of runs scored."""
    outs = runs = 0     # Inning starts with no outs and no runs,
    runners = set()     # and with nobody on base
    def out(r) -> int: runners.remove(r); return 1
    def forced(r) -> bool: return all(b in runners for b in range(r))
    def advance(r, e) -> int: 
        return int(e if e in '1234' else (e in 'ESD' or (e == 'B' and forced(r))))
    for e in events:
        if verbose: show(outs, runs, runners, e)
        runners.add(batter)      # Batter steps up to the plate
        if e == 'D' and len(runners) > 1: # Double play: batter and lead runner out
            outs += out(batter) + out(max(runners))
        elif e in 'DSKOF':       # Batter is out
            outs += out(batter)  
        if outs >= 3:            # If inning is over: return runs scored
            return runs 
        runners = {r + advance(r, e) for r in runners} # Runners advance
        runs += len(runners & scored)                  # Tally runs
        runners = runners - scored                     # Remove runners who scored
        
def event_stream(events, strikes=0):
    """A generator of random baseball events."""
    while True:
        yield 'K' if (random.random() < strikes ** 3) else random.choice(events)

def show(outs, runs, runners, event):
    """Print a representation of the current state of play."""
    bases = ''.join(b if int(b) in runners else '-' for b in '321')
    print(f'{outs} outs {runs} runs  {bases}  {event} ({event_codes[event]})')
    
ONBG   = event_stream('2111111EEBBOOSSSSSSSFFFFFD334', 7/36) # Our National Ball Game
batter = 0             # The batter is not yet at first base
scored = {4, 5, 6, 7}  # Runners in these positions have scored


# # Examples and Tests
# 
# Let's peek at some random innings:

# In[4]:


inning(ONBG)


# In[5]:


inning(ONBG)


# Let's also test some historic innings. I'll take some of the Red Sox innings from their 2004 playoff series against the Yankees.

# In[6]:


# 7th inning in game 1: 5 runs (Homer by Varitek)
# (But not a perfect reproduction, because our simulation doesn't have passed balls.)
assert 5 == inning('K2O1214K')


# In[7]:


# 4th inning in game 6: 4 runs (Homer by Bellhorn)
assert 4 == inning('SS2114F')


# In[8]:


# 2nd inning in game 7: 4 runs (Grand Slam by Damon)
assert 4 == inning('S1BB4BFS')


# That looks good to me.
# 
# # Simulation
# 
# Now, simulate a hundred thousand innings, and then sample from them to simulate a hundred thousand nine-inning games (for one team), and show histograms of the results, labelled with statistics:

# In[9]:


def simulate(N=100000, inning=inning, events=ONBG) -> None:
    innings = [inning(events=events, verbose=False) for _ in range(N)]
    games   = [sum(random.sample(innings, 9)) for _ in range(N)]
    hist(innings, 'Runs/inning (for one team)')
    hist(games,   'Runs/game (for one team)')
    
def hist(nums, title): 
    """Plot a histogram and show some statistics."""
    plt.hist(nums, ec='black', bins=max(nums)-min(nums), align='left')
    plt.xlabel(title)
    plt.title(f'μ: {mean(nums):.2f}, σ: {stdev(nums):.2f}, max: {max(nums)}')
    plt.show()


# In[10]:


get_ipython().run_line_magic('time', 'simulate()')


# So, about 13 runs per game (per team). This shows that the dice game is not very realistic with respect to current-day baseball. It is true that games were higher-scoring 130 years ago, and perhaps a dice game is more fun when there is a lot of action.
# 
# # Real Major League Baseball Stats
# 
# Could I make the game reflect baseball as it is played today? To do so I would need:
# 1. A source of major league baseball (MLB) statistics.
# 2. A way to convert those statistics into the format expected by the function `inning`.
# 3. Possibly some modifications to  `inning`, depending on how the conversion goes.
# 
# [Baseball-reference.com](https://www.baseball-reference.com) has lots of stats, in particular 
# [MLB annual batting stats](https://www.baseball-reference.com/leagues/MLB/bat.shtml) and
# [fielding stats](https://www.baseball-reference.com/leagues/MLB/field.shtml); I'll use the stats for the complete 2019 season. The batting stats have most of what we need, and the fielding stats give us double plays and errors.
# 
# I start by defining two utility functions that can be useful for any tabular data: `cell_value`, which converts a table cell entry into an `int`, `float`, or `str` as appropriate; and `header_row_dict`, which creates a dict of `{column_name: value}` entries. The function `mlb_convert` then converts this format (a dict keyed by `H/2B/3B/HR` etc.) into the event code format (a string of `'1234...'`). As part of the conversion  I'll add hit-by-pitch (`HBP`) into the "base on balls" category, and I'll record all otherwise  unaccounted-for outs under the "fly out" (`F`) category (runners do not advance). With this understood, we won't need to change the function `inning` at all. (It is true that `mlb_convert` returns a very long string, equal in length to the number of plate appearances over the whole MLB season. But that takes up less space than storing one photo, so I'm not going to worry about it.)

# In[11]:


def cell_value(entry, types=(int, float, str)):
    """Convert a cell entry into the first type that doesn't raise an error."""
    for typ in types:
        try:
            return typ(entry)
        except ValueError:
            pass
    
def header_row_dict(header, row, sep=None, value=cell_value) -> dict:
    """Parse a header and table row into a dict of `{column_name: value(cell)}`."""
    return dict(zip(header.split(sep), map(value, row.split(sep))))

def mlb_convert(stats: dict) -> str:
    """Given baseball stats return a string '11...FFF'."""
    events = Counter({
        '1': stats['H'] - stats['2B'] - stats['3B'] - stats['HR'],
        '2': stats['2B'], '3': stats['3B'], '4': stats['HR'],
        'E': stats['E'],  'B': stats['BB'] + stats['HBP'],
        'K': stats['SO'], 'D': stats['DP'], 'S': stats['SH'] + stats['SF']})
    events['F'] = stats['PA'] - sum(events.values()) # All unaccounted-for outs
    return ''.join(events.elements())                # A str of events


# Below I copy-and-paste the data I need from baseball-reference.com to create the dict `mlb_stats`; convert it to the string `mlb_string`; and use that to create the event generator `mlb_stream`:

# In[12]:


mlb_stats = header_row_dict(
    "Year Tms #Bat BatAge R/G G PA AB R H 2B 3B HR RBI SB CS BB SO BA OBP SLG OPS TB GDP HBP SH SF IBB E DP",
    """2019 30 1284 27.9 4.84 4828 185377 165622 23346 41794 8485 783 6735 22358 2261 827 15806 42546 
    .252 .323 .435 .758 72050 3441 1968 774 1146 752 2882 3981""")
mlb_string = mlb_convert(mlb_stats)
mlb_stream = event_stream(mlb_string)


# We can take a look:

# In[13]:


mlb_stats


# In[14]:


mlb_string[::1000] # Just look at every 1000th character


# In[15]:


inning(mlb_stream)


# I can simulate:

# In[16]:


get_ipython().run_line_magic('time', 'simulate(events=mlb_stream)')


# That looks a *lot* more like real baseball. But MLB averaged 4.84 runs per team per game in 2019, and this is significantly lower. I think we can make some minor changes to the function `inning`&mdash;some "standard baseball things"&mdash;to make the simulation more realistic. I'm thinking of two changes:
# - The most common double play eliminates the batter and the runner on first, not the runner closest to home. 
# - On a single, a runner on second often scores. 
# 
# I'll make those two things the case for all double plays and singles.

# In[17]:


def inning2(events, verbose=False) -> int:
    """Simulate a half inning based on events, and return number of runs scored."""
    outs = runs = 0     # Inning starts with no outs and no runs,
    runners = set()     # and with nobody on base
    def out(r) -> int: runners.remove(r); return 1
    def forced(r) -> bool: return all(b in runners for b in range(r))
    def advance(r, e) -> int: 
        return ((2 if r == 2 else int(e)) if e in '1234' else 
                (e in 'ESD' or (e == 'B' and forced(r))))    
    for e in events:
        if verbose: show(outs, runs, runners, e)
        runners.add(batter)           # Batter steps up to the plate
        if e == 'D' and 1 in runners: # Double play: batter and runner on first out
            outs += out(batter) + out(1)
        elif e in 'DSKOF':            # Batter is out
            outs += out(batter)  
        if outs >= 3:                 # If inning is over: return runs scored
            return runs 
        runners = {r + advance(r, e) for r in runners} # Runners advance
        runs += len(runners & scored)                  # Tally runs
        runners = runners - scored                     # Remove runners who scored


# We show the difference with two examples. First, a triple/walk/double-play sequence scores a run under `inning2` but not `inning`:

# In[18]:


inning2('3BDK', True)


# In[19]:


inning('3BDK', True)


# Second, a double/single sequence scores a run under `inning2` but not `inning`:

# In[20]:


inning2('21FFF', True)


# In[21]:


inning('21FFF', True)


# We can simulate again and note any differences:

# In[22]:


get_ipython().run_line_magic('time', 'simulate(events=mlb_stream, inning=inning2)')


# There is a slight increase in the number of runs.
# 
# # Opportunities for Improvement
# 
# There are many problems with the code as it is. For example:
# 
# - It assumes all teams are equal. They're not.
# - It assumes all pitchers and defense are equal. They're not.
# - It assumes all batters are equal. They're not. (I think this is the main place where we get a shortfall in runs: real lineups cluster their best hitters together, and they are more apt to produce runs than a lineup of all median players.)
# - It assumes all hits are the same (runners always advance the same number of bases). They're not.
# - There's only one type of double play (batter and runner on first out) and no triple play.
# - It ignores stolen bases, pickoffs, passed balls, wild pitches, runners taking extra bases, and runners being out on attempted steals, extra bases, or sacrifices.
# - There is no strategy (offense and defense behave the same, regardless of the situation).
# - It assumes both teams bat for 9 innings. But if the home team is ahead at the bottom of the 9th, they do not bat, and if the score is tied: extra innings.
# - With two outs, or with no runners on base, there can be no sacrifice or double play; those types of events would just be regular outs. The stats say that a double play should occur in 3981 out of 185377 at bats, or about 2% of the time. In our simulation the `D` event code would come up that often, but perhaps only half the time there would be a runner and less than two outs, so we would only actually get a double play maybe 1% of the time.
# 
# 
# What can you do to make the simulation better?