Reverse Engineering xkcd's 'Frequency'

import requests
import lxml.html
import gifparse

import pandas as pd
import numpy as np

from IPython.display import HTML
from markdown import markdown

def md(x, *args, **kwargs):
    return HTML(markdown(x.format(*args, **kwargs)))

def commafy(n):
    if type(n) == int:
        return "{:,d}".format(n)
    return str(n)

import matplotlib as mpl
import matplotlib.pyplot as plt
import mplstyle
%matplotlib inline

def gif_to_html(raw_bytes):
    """Convert the raw bytes of a GIF into an HTML <img>."""
    b64 = raw_bytes.encode("base64")
    tmpl = '<img src="data:image/gif;base64,%s">'
    return tmpl % b64

def get_image_urls():
    """Fetch the URLs for all GIFs from the xkcd's 'Frequency' piece."""
    URL = "http://xkcd.com/1331/"
    html = requests.get(URL).content
    dom = lxml.html.fromstring(html)
    images = dom.cssselect("#comic td img")
    srcs = [ i.attrib["src"] for i in images ]
    return srcs

### NOTE
# This analysis uses a custom GIF parser, which you can see here: 
#
#   https://github.com/jsvine/gifparse
# 
# I wrote the parser in part because the PIL library — the standard 
# tool for this sort of thing — seemed to be having trouble with 
# the particular format of the GIFs in "Frequency." (I also just wanted
# to learn the spec.) As a result, `gifparse` is an incomplete library.
# I'm hoping to beef it up soon, but feel free to send pull requests
# in the meantime.
###
class Scenario(object):
    """A scenario from xkcd's 'Frequency' piece.
    
    The corresponding GIF is downloaded upon instantiation.
    """
    def __init__(self, url):
        self.url = url
        self.name = url.split("/")[-1].split(".")[0]
        self.raw = requests.get(url).content
        self.gif = gifparse.parse(self.raw)
    def __repr__(self):
        return "<Scenario:{}>".format(self.name)

# Instantiate 50 scenarios, one for each GIF in the comic
image_urls = get_image_urls()
scenarios = map(Scenario, image_urls)

# Note: GIF delays are, by spec, multiples of 1/100th of a second.
freq = pd.DataFrame({
    "scenario": s,
    "name": s.name,
    "loop_seconds": round(float(s.gif.total_delay) / 100, 2)
} for s in scenarios).set_index("name")
freq["per_second"] = (1 / freq["loop_seconds"]).apply(lambda x: round(x, 4))
freq["per_hour"] = (60*60 / freq["loop_seconds"]).apply(lambda x: round(x, 2))
freq["heartbeats"] = (freq["loop_seconds"] / freq.ix["heartbeat"]["loop_seconds"])\
    .apply(lambda x: round(x, 4))

grid = pd.DataFrame(np.split(freq["per_hour"].values, 10))
grid_log = grid.apply(np.log10)

def create_heatmap():
    pixel_coeff = 1.5 / mplstyle.get("figure.dpi")
    grid_dim = (120 * 5 + 120, 50 * 10)
    mplstyle.set({
    "figure.figsize": tuple(pixel_coeff * x for x in grid_dim)
    })
    fig, ax = plt.subplots()
    heatmap = ax.pcolor(grid_log, 
        cmap=mpl.cm.Blues,
        edgecolors="white")
    ax.invert_yaxis()
    ax.set_xticks([])
    ax.set_yticks([])
    for row_i in range(len(grid)):
        for col_i in range(len(grid.ix[row_i])):
            index = row_i * 5 + col_i
            name = freq.index[index]
            per_hour = freq["per_hour"][index]
            ax.text(col_i + 0.5, row_i + 0.5, name, 
                ha="center", 
                va="center", 
                fontsize=12,
                color="white" if per_hour > 50 else "black")
    max_freq = int(np.log10(freq["per_hour"].max())) + 1
    ticks = range(0, max_freq)
    cbar = plt.colorbar(heatmap, ticks=ticks)
    cbar.set_ticklabels([ commafy(pow(10, int(x))) for x in ticks ])
    cbar.set_label("Occurances per Hour",
        rotation=270,
        fontsize=12)
create_heatmap()

def plot_distribution():
    mplstyle.set({ "figure.figsize": (8, 6) })
    ticks = range(-1, 6)
    ax = grid_log.unstack().hist(bins=ticks, color=mpl.cm.Blues(mpl.cm.Blues.N / 2))
    ax.set_xticklabels([ commafy(pow(10, x)) for x in ticks ])
    ax.set_title("Distribution of Frequencies\n", fontweight="bold")
    ax.set_xlabel("\nOccurances per Hour")
    ax.set_ylabel("Number of Scenarios\n")
    pass
plot_distribution()

HTML("".join(gif_to_html(freq.ix[name]["scenario"].raw) 
             for name in freq.sort("loop_seconds").index[:5]))

HTML("".join(gif_to_html(freq.ix[name]["scenario"].raw) 
             for name in freq.sort("loop_seconds").index[-5:]))

old_faithful_mins = freq.ix["oldfaithful"]["loop_seconds"] / 60
md("""(You'll have to wait {} minutes for Old Faithful.)""".format(int(old_faithful_mins)))

def display_twinned_gifs():
    time_grouped = freq.groupby("loop_seconds")
    twinned_times = time_grouped.size()[time_grouped.size() > 1].index
    html = ""
    for time in twinned_times:
        matching_scenarios = freq[freq["loop_seconds"] == time]
        html += "<h4>Every {} seconds:</h4>".format(time)
        html += "".join(gif_to_html(x.raw) for x in matching_scenarios["scenario"])
    return HTML(html)
display_twinned_gifs()

def get_delays(gif):
    return [ x["gce"].delay for x in gif.frames ]
freq["zeros"] = freq["scenario"].apply(lambda x: sum(pd.Series(get_delays(x.gif)) == 0))
freq["loop_seconds_adj"] = freq["loop_seconds"] + freq["zeros"] * 0.1
freq["per_hour_adj"] = (60*60 / freq["loop_seconds_adj"]).apply(lambda x: round(x, 2))

signals = freq.ix[[ "turnsignal1", "turnsignal2" ]]
t1, t2, = signals["loop_seconds"]
t1_adj, t2_adj = signals["loop_seconds_adj"]

HTML("".join(gif_to_html(x.raw) for x in signals["scenario"]))

signal_beat = 1 / np.abs((1.0 / t1) - (1.0 / t2))
md(u"""
According to the data embedded *Frequency*'s GIFs, 
Randall Munroe's turn signal blinks once every {:.2f} seconds, 
while the turn signal in front of him blinks once every {:.2f} seconds — just out of sync.
How often will the two signals blink together? Using the equation for 
[beat frequency](http://hyperphysics.phy-astr.gsu.edu/hbase/sound/beat.html), we can calculate this.
The two signals should blink together approximately every {:.2f} seconds.
""", t1, t2, signal_beat)

signal_beat_adj = 1 / np.abs((1.0 / t1_adj) - (1.0 / t2_adj))
md(u"""
But you'll see something different. If you're using a modern browser in 2014, 
you'll observe a beat frequency closer to {:.2f} seconds — about {:.2f} seconds
longer than expected. Why is that?
""", signal_beat_adj, signal_beat_adj - signal_beat)

min_loop_seconds = freq["loop_seconds"].min()
max_loop_seconds = freq[freq["zeros"] > 0]["loop_seconds"].max()
md("""
The __zero-delay problem affects {} of the 50 GIFs__ in *Frequency*. 
When the browser renders them, the loops will take 1/10th-second longer than
it seems Munroe intended. The affected GIFs happen to be all those shorter than 
two seconds, which might be an artifact of how Munroe created them. 
(I emailed Munroe early in my research, before discovering this issue, but never heard back.
I'll update this post if I do hear from him.) For the fastest affected GIF, that's a 
{:.1f}% lengthening. For the slowest, the zero-delay problem adds {:.1f}% to its encoded loop time.
""",
    sum(freq["zeros"] > 0),
    100 * (((0.1 + min_loop_seconds) / min_loop_seconds) - 1),
    100 * (((0.1 + max_loop_seconds) / max_loop_seconds) - 1)
)

freq["needs_adj"] = freq["zeros"].apply(lambda x: "*" if x else "")
freq.sort("loop_seconds")[["loop_seconds", "per_hour", "needs_adj", "loop_seconds_adj", "per_hour_adj"]]