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

# # Cellular Automaton
# 
# John Conway's [Game of Life](https://en.wikipedia.org/wiki/Conway's_Game_of_Life) is something lots of people have probably seen. Its behavior is intriguing because it exhibits complex behavior with simple rules. I want to use his game as a starting point for an exploration of cellular automaton.
# 
# I will assume some basic knowledge of programming (arrays, loops, etc.) but will try and make the code as clear as possible because I think reading code can sometimes explain the details quite well. Oh, and lots of pictures, because I like pictures. 
# 
# Let's get started.

# In[427]:


# this makes plots show up as inline images
get_ipython().run_line_magic('matplotlib', 'inline')
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd # we only use pandas for its pretty table displays
from itertools import islice, product, takewhile
from functools import partial
from scipy.signal import convolve2d
from collections import Counter
from random import random


# ## The grid

# A [cellular automaton](https://en.wikipedia.org/wiki/Cellular_automaton) is a grid of cells where each cell has a state. For the Game of Life, we will use a square grid and each cell has 2 states: dead (0, white) or alive (1, black).

# In[428]:


grid = lambda n: np.zeros((n, n), dtype='int')


# In[429]:


# this grid is all dead
print(grid(8))


# In[430]:


# this grid as all alive
print(grid(8) + 1) # <- the `+ 1` gets added to all the cells


# In[431]:


# make a grid with (approximately) `p` percent of cells alive
# `p` takes values [0.00, 1.00]
random_grid = lambda n, p: (np.random.ranf((n, n)) < p).astype('int')


# In[432]:


# this grid will have ~10% starting cells alive
g = random_grid(8, 0.10)
print('There are {} alive cells --> {} percent alive'.format(g.sum(), g.sum() / g.size))
print(g)


# In[433]:


# let's get our pictures going
def show_grid(g):
    im = plt.imshow(g, cmap=plt.cm.Greys, interpolation='None')
    plt.tick_params(labelbottom='off', labelleft='off', left='off', right='off', top='off', bottom='off')
    return im


# In[434]:


g = random_grid(8, 0.1)
print(g)
show_grid(g);


# ## Iteration
# 
# A nice way to think about the Game of Life is as a function. It takes a grid as an input and returns a new grid as an output. In math, we would write this as something like:
# 
# $$step : G \rightarrow G$$
# $$g \mapsto g'$$
# 
# where `step` is the function we will define and `G` just means, "the set of all possible grids"
# 
# The Game of Life rules are defined by
# 
#   1. Any live cell with fewer than two live neighbours dies, as if caused by under-population.
#   2. Any live cell with two or three live neighbours lives on to the next generation.
#   3. Any live cell with more than three live neighbours dies, as if by over-population.
#   4. Any dead cell with exactly three live neighbours becomes a live cell, as if by reproduction.
# 
# Neighbors are the set of 8 cells immediately adjacent to some cell.

# In[435]:


# the black cells are neighbors of the one middle white cell
grid1 = np.array([[0, 0, 0, 0, 0],
                  [0, 1, 1, 1, 0],
                  [0, 1, 0, 1, 0],
                  [0, 1, 1, 1, 0],
                  [0, 0, 0, 0, 0]])
show_grid(grid1);


# One thing we need to consider is what the neighbors of the upper left cell would be? One possibility is to assume anything past the grid is dead. Another is that the cells wrap around like those old asteroid games. Let's call the first mode `clip` and the second `wrap`

# In[436]:


def n_live_neighbors(g, i, j, mode):
    """Count the live neighbors of (i, j) under the given `mode`,
    `mode` should be one of `clip` or `wrap`"""
    N = g.shape[0]
    live = 0
    for di in [-1, 0, 1]:
        for dj in [-1, 0, 1]:
            if di == 0 and dj == 0: # skip the cell at (i, j)
                continue
            if mode == 'wrap':
                live += g[(i + di) % N, (j + dj) % N]
            else: # mode == 'clip'
                if not (i + di) in range(N) or (j + dj) not in range(N):
                    continue # adding 0 is equivalent to ignoring the out of bound cells
                live += g[i + di, j + dj]
    return live


# In[437]:


# note that array indexing starts at 0
print('The middle cell (2, 2) has {} live neighbors'.format(n_live_neighbors(grid1, 2, 2, 'clip')))
print('The bottom left cell (4, 0) has {} live neighbors'.format(n_live_neighbors(grid1, 4, 0, 'clip')))


# Now we can write a function which produces the next state of each cell.

# In[438]:


def step(g, mode):  
    """Take one step in the Game of Life, 
    `mode` should be one of `clip` or `wrap`"""
    N = g.shape[0]
    g_prime = grid(N) # this is an all dead grid, so we only have to worry about setting live cells
    for i in range(N):
        for j in range(N):
            live = n_live_neighbors(g, i, j, mode)
            if g[i, j] == 1: # this cell is alive
                if live == 2 or live == 3:
                    g_prime[i, j] = 1
                # otherwise the cell dies (is already set to 0 by default)
            else: # this cell is dead
                if live == 3:
                    g_prime[i, j] = 1
                # otherwise it stays dead (already set to 0 by default)
    return g_prime

clip_step = lambda g: step(g, 'clip')
# the above is equivalent to partial(step, mode='clip') in case you're curious
wrap_step = lambda g: step(g, 'wrap')


# In[439]:


g = random_grid(16, 0.1)
g_prime_clip = clip_step(g)
g_prime_wrap = wrap_step(g)

plt.figure(figsize=(16,9))

plt.subplot(131)
plt.title('Original state')
show_grid(g);

plt.subplot(132)
plt.title('Next state mode=clip')
show_grid(g_prime_clip);

plt.subplot(133)
plt.title('Next state mode=wrap')
show_grid(g_prime_wrap);

# there may or may not be differences between the two modes in this random example


# So let's just stick with `clip` mode for now and look what happens when we advance many steps in a row. Since `step` gives us back a grid of the same size, we can simply call `step(step(g))` to take two steps.

# In[440]:


g = random_grid(16, 0.1)
g_prime = clip_step(g)
g_prime_prime = clip_step(g_prime)

plt.figure(figsize=(16,9))

plt.subplot(131)
plt.title('Original state $t = 0$')
show_grid(g);

plt.subplot(132)
plt.title('Next state $t = 1$')
show_grid(g_prime);

plt.subplot(133)
plt.title('Next Next state $t = 2$')
show_grid(g_prime_prime);


# To see what happens as `t` gets larger, lets get a better visual going, like a movie! But first a function which takes a grid `g` and will yield `n` intermediate grids at each time step.

# In[441]:


def step_n_try1(g0, N, mode):
    yield g0
    g = g0
    for t in range(N):
        g = step(g, mode)
        yield g


# As an aside, this was the first version of the function I wrote; it does exactly what I want but is very specific to the Game of Life and the `step` function and its `mode` option. Let's try for something more general.

# In[442]:


def step_n_try2(f, x0, N):
    yield x0
    x = x0
    for t in range(N):
        x = f(x)
        yield x


# Now `step_n_try2` only cares about an arbitrary (but single argument) function `f` and the initial value, `x0`. This bit of abstraction would let us reuse this function in another context. But we can still do better. We'll first define a a function `step_many` to take steps forever.

# In[443]:


def step_many(f, x0):
    yield x0
    x = x0
    while True:
        x = f(x)
        yield x


# Finally we can see that the standard Python function `itertools.islice` will let us take the first `N` of those. This gives a general set of functions which can easily be applied to other problems.

# In[444]:


step_n = lambda f, x0, N: islice(step_many(f, x0), N)


# In[445]:


# Don't worry about anything in here, just connecting pipes together to make movies show up in this notebook
# adapted from source: http://nbviewer.ipython.org/gist/edrex/9044756/mpl_animation_html.ipynb
from tempfile import NamedTemporaryFile
from IPython.display import HTML
from base64 import b64encode
import matplotlib.animation as animation

VIDEO_TAG = """<video controls autoplay loop>
 <source src="data:{0}">
 Your browser does not support the video tag.
</video>"""

def anim_to_html(anim, fps):
    if not hasattr(anim, '_encoded_video'):
        with NamedTemporaryFile(suffix='.m4v') as f:
            anim.save(f.name, fps=fps, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
            video = open(f.name, "rb").read()
        anim._encoded_video = 'video/mp4;base64,' + b64encode(video).decode('utf-8')
    # prevent figure displayed as a PNG below the animation
    plt.close()
    
    return HTML(VIDEO_TAG.format(anim._encoded_video))

def grid_anim(gs):
    fig = plt.figure()
    def frames():
        for t, g in enumerate(gs):
            # TODO: get this to work with the title
            # yield (show_grid(g), plt.title('$t = {}$'.format(t)))
            yield (show_grid(g), )
    return animation.ArtistAnimation(fig, list(frames()), blit=True)

def show_grid_anim(gs, fps=2):
    return anim_to_html(grid_anim(gs), fps)


# In[446]:


# here we have 100 steps of a random grid at 3 frames per second
show_grid_anim(step_n(clip_step, random_grid(32, 0.1), 50), fps=3)


# Its interesting to see that some grids seem to get stuck, like the one above. While this makes for a boring movie, it is very interesting because we just discovered a fix point!
# 
# A [fixed point](https://en.wikipedia.org/wiki/Fixed_point_%28mathematics%29) of a function `f` is an `x` where
# 
# $$ f(x) = x $$
# 
# A fixed point everyone knows from algebra would be
# 
# $$f(x) = x^2$$
# 
# where `1` is a fixed point of `f` because
# 
# $$1^2 = 1$$
# 
# And so a fixed point in the Game of Life is a grid `g` which steps to itself.
# 
# (I'll note that fixed points appear in programming language theory, for example the [Y combinator](https://en.wikipedia.org/wiki/Fixed-point_combinator#Fixed_point_combinators_in_lambda_calculus), completely tagent to the topic of this notebook but interesting to note the connection nonetheless)

# Lets try and find some fix points of the Game of Life! First with some brute force on *all* 5x5 grids.

# In[447]:


# is x a fixpoint of f
# due to numpy internals, this is not as general as it could be, but it will suffice
is_fixpoint = lambda f, x: np.all(f(x) == x)
# don't worry about the insides of this, just know it yields every possible board state of size n
nxn_grids = lambda n: map(lambda l: np.array(l).reshape((n, n)), product(*[range(2)] * n ** 2))


# In[448]:


# a quick test, there should be 16 possible 2x2 boards
for i, g in enumerate(nxn_grids(2)):
    plt.subplot(4, 4, i + 1)
    show_grid(g)


# An `n x n` grid will have
# $$ 2 ^{n ^ 2}$$
# possible cell configurations since there are `n x n` grid cells and each can take one of 2 values. This grows very very very quickly.

# In[449]:


ns = np.arange(2, 8)
pd.DataFrame(2 ** (ns ** 2), index=ns)


# Lets print out all configurations in `4x4` grids which are fixpoints.

# In[450]:


# only keep the 4x4 grids which are fixpoints
fixpoints = filter(lambda g: is_fixpoint(clip_step, g), nxn_grids(4))
plt.figure(figsize=(16, 9))
n = 0
for i, g in enumerate(fixpoints):
    n += 1
    plt.subplot(10, 16, i + 1)
    show_grid(g)
print('{}% of possible 4x4 grids are fixpoints'.format(100 * n / (2 ** (4 ** 2))))


# These fixpoints are not the only interesting patterns in the Game of Life. Another is repetitive sequences of grids which end up back at the same point after `n` steps. A 2-cycle is something where (read `|` as "such that"):
# 
# $$ x \mid f(f(x)) = x $$
# 
# But it may take a lot of steps of `f` to get back to the original `x`. We can say that `x` is an `n` length cycle if (read f ^ n as n applications of f):
# 
# $$ x \mid f^{n}(x) = x $$
# 
# This notion of cycles is related to something like the [Lorenz System](https://en.wikipedia.org/wiki/Lorenz_system).

# Lets try and find some of these cycles!

# In[451]:


def possible_cycle(it):
    """return None if `it` only has unique values
    otherwise, return the first instance of a value which is repeated"""
    # the call to `str` is a cheat since numpy.ndarray is unhashable
    cache = set(str(next(it))) # store every grid we've seen in our iterations
    for el in it:
        if str(el) in cache: # have we seen this already?
            return el
        else:
            cache.add(str(el))
    return None

def cycle(it):
    """We assume the first element in `it` forms a finite cycle,
    yield the elements which form a cycle"""
    x0 = next(it)
    yield x0
    for el in it:
        if np.all(el == x0): # this breaks generality, but it will suffice
            return el
        else:
            yield el


# Lets print the first non-simple cycle we find in 4x4 grids.

# In[452]:


MAX_STEPS = 100 # how long we'll look for a cycle
for g in nxn_grids(4):
    g_prime = possible_cycle(step_n(clip_step, g, MAX_STEPS))
    if g_prime is not None: # we have a cycle beginning/ending at g_prime
        gs = list(cycle(step_many(clip_step, g_prime)))
        if len(gs) == 1:
            continue # skip simple cycles aka fixed-points
        gs.append(gs[0]) # add the first grid to the end to make a seamless video
        break
show_grid_anim(gs, fps=1)


# In[453]:


def n_cycle(it, n, max_steps=100):
    for g in it:
        g_prime = possible_cycle(step_n(clip_step, g, max_steps))
        if g_prime is not None: # we have a cycle beginning/ending at g_prime
            gs = list(cycle(step_many(clip_step, g_prime)))
            if len(gs) != n:
                continue # skip simple cycles aka fixed-points
            gs.append(gs[0]) # add the first grid to the end to make a seamless video
            return g, gs


# In[454]:


g, gs = n_cycle(nxn_grids(4), 3)
show_grid_anim(gs, fps=1)


# Turns out, its rather difficult to find n-cycles with such simple methods as the number of possible starting grids to search through is enormous. See the Wikipedia page on Game of Life for examples of oscillators with periods of up to 15!

# We'll be transitioning away from the Game of Life (not before two tangents though), so I'll include a nice long video of 100 steps of a random 32x32 grid to bask in its glory. If you're interested in larger simulations, there are programs specialized in efficient simulation of large grids and for many time steps.

# In[455]:


show_grid_anim(step_n(clip_step, random_grid(32, 0.5), 100))


# The following is something which struck me as awesome the first time I saw it as it connected two otherwise disparate areas of study with such simplicity. I first saw this [here](http://kamicut.cc/2014/05/29/conways-game-of-life-julia.html).
# 
# Remember that the Game of Life is computed for each cell based on its neighboring cells. If you imagine a 3x3 grid centered over and hovering above the current cell, the cell's neighbors are exactly those which line up with the hovering grid. If we want to count the number of live cells in the neighbors, we can fill the hovering grid with all 1s (and the middle cell is a 0) and compute a dot-product-like sum with the original grid. I think this image sums it up best (the hovering grid is called a kernel):
# 
# ![Convolution image](https://developer.apple.com/library/ios/documentation/Performance/Conceptual/vImage/Art/kernel_convolution.jpg)
# 
# This process is a fundamental process in image processing, you can read more [here](https://en.wikipedia.org/wiki/Kernel_%28image_processing%29)
# 
# I bring this up because it greatly simplifies our function. We no longer have to worry about choosing the proper indices, but only about producing a new "image".

# In[456]:


def step_convolution(g, mode):
    """Use the same modes as before"""
    # The source I reference uses a slightly different kernel and is a bit trickier, see it if you're interested
    kernel = np.array([[1, 1, 1], [1, 0, 1], [1, 1, 1]])
    if mode == 'clip':
        live = convolve2d(g, kernel, mode='same', boundary='fill', fillvalue=0)
    else: # mode == 'wrap'
        live = convolve2d(g, kernel, mode='same', boundary='wrap')
    # live[i, j] contains the number of live neighbors of the cell g[i, j]
    # the following are the combinations of what will yield a live cell
    return (((g == 1) & (live == 2)) | ((g == 1) & (live == 3)) | ((g == 0) & (live == 3))).astype('int')
 


# In[457]:


show_grid_anim(step_n(lambda g: step_convolution(g, 'clip'), random_grid(32, 0.5), 100))


# So you can see we get the same result, but in a different way; neat!

# The last tangent before moving past Game of Life is the mode we've been ignoring this whole time: `wrap`.
# 
# Remember that `wrap` treats a grid as if cells on the left edge are next to the cells on the right edge and likewise for the top and bottom edges. First let's see an example of this mode in action

# In[458]:


show_grid_anim(step_n(lambda g: step_convolution(g, 'wrap'), random_grid(32, 0.2), 100))


# A grid which exhibits this wrapping property behaves just like a [torus](https://en.wikipedia.org/wiki/Torus) (aka donut)!
# 
# Let's view the Game of Life running on a torus! We can map a cartesian grid onto the surface of a torus as follows:
# $$ u : [0, 2\pi]$$
# $$ v : [0, 2\pi]$$
# 
# $$x(u, v) = (R + r\cos(u))\cos(v)$$
# $$y(u, v) = (R + r\cos(u))\sin(v)$$
# $$z(u, v) = r\sin(u)$$
# 
# Where `R` is the distance from center of the tube to the center of the torus
# and `r` is the radius of the tube

# In[459]:


R = 4
r = 3
u = np.linspace(0, 2 * np.pi, 200)
v = np.linspace(0, 2 * np.pi, 200)
uu, vv = np.meshgrid(u, v)
x = (R + r * np.cos(uu)) * np.cos(vv)
y = (R + r * np.cos(uu)) * np.sin(vv)
z = r * np.sin(uu)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, z, alpha=0.2);


# In[460]:


def grid_to_torus(g, R, r):
    N = g.shape[0]
    u = np.linspace(0, 2 * np.pi, N)
    v = np.linspace(0, 2 * np.pi, N)
    x = []
    y = []
    z = []
    for i in range(N):
        for j in range(N):
            if g[i, j] == 1:
                uu = u[i]
                vv = v[j]
                x.append((R + r * np.cos(uu)) * np.cos(vv))
                y.append((R + r * np.cos(uu)) * np.sin(vv))
                z.append(r * np.sin(uu))
    return x, y, z

def torus_surface(R, r):
    u = np.linspace(0, 2 * np.pi, 200)
    v = np.linspace(0, 2 * np.pi, 200)
    uu, vv = np.meshgrid(u, v)
    x = (R + r * np.cos(uu)) * np.cos(vv)
    y = (R + r * np.cos(uu)) * np.sin(vv)
    z = r * np.sin(uu)
    return x, y, z
 


# In[461]:


# add some 3d versions of our show function, they are slightly optimized to
# only plot the surface once, which makes them repetitive

def clear_ticks_labels(ax):
    ax.set_xticklabels([])
    ax.set_xticks([])
    ax.set_yticklabels([])
    ax.set_yticks([])
    ax.set_zticklabels([])
    ax.set_zticks([])

def show_grid3(g, R=4, r=3):
    ax = plt.gca(projection='3d')
    clear_ticks_labels(ax)
    x, y, z = torus_surface(R, r)
    ax.plot_surface(x, y, z, alpha=0.1)
    x, y, z = grid_to_torus(g, R, r)    
    return ax.scatter(x, y, z)

def grid3_anim(gs, R=4, r=3):
    fig = plt.figure()
    ax = plt.gca(projection='3d')
    clear_ticks_labels(ax)
    x, y, z = torus_surface(R, r)
    ax.plot_surface(x, y, z, alpha=0.1)
    def frames():
        for t, g in enumerate(gs):
            yield (ax.scatter(*grid_to_torus(g, R, r)), )
    return animation.ArtistAnimation(fig, list(frames()), blit=True)

def show_grid3_anim(gs, fps=2):
    return anim_to_html(grid3_anim(gs), fps)


# Let's test it!

# In[462]:


show_grid3(random_grid(32, 0.1));


# In[463]:


show_grid3_anim(step_n(wrap_step, random_grid(32, 0.2), 50))


# From here, we will depart from the Game of Life and venture into more general territory, starting with a stochastic, or random element.
# 
# Instead of considering a fixed rule of state transitions, we can think of assigning a probability of transitioning from one state to the next based on some information from the neighbors of each cell. We'll only consider `clip` mode from here on out (assume `0` is the default value for the edges) and will generalize the stepping method. In anticipation of something in the future, we will write this function a bit more generally. We will allow the grid to have any number of states, and the transition function should take a cell's current state, the tally of the states of its neighbors, and produce a new state for that cell.

# In[464]:


def iter_neighbors(g, i, j):
    """Iterate over the 8 neighbors of g[i, j]"""
    N = g.shape[0]
    counter = Counter() # counter simply tallies how many of each kind it sees
    for di in [-1, 0, 1]:
        for dj in [-1, 0, 1]:
            if di == 0 and dj == 0: # skip the cell at (i, j)
                continue
                live += g[(i + di) % N, (j + dj) % N]
            if not (i + di) in range(N) or (j + dj) not in range(N):
                yield 0
            else:
                yield g[i + di, j + dj]

# return a tally of the number of each state for the 8 neighbors            
count_neighbors = lambda g, i, j: Counter(iter_neighbors(g, i, j))

# stoch == stochastic
def stoch_step(g, transition):
    """`g` is a grid and `transition` is a function which takes the current state and a
    Counter object and returns a new state"""
    N = g.shape[0]
    g_prime = grid(N)
    for i in range(N):
        for j in range(N):
            g_prime[i, j] = transition(g[i, j], count_neighbors(g, i, j))
    return g_prime


# To see this is a generalization, we can easily write the rules for Game of Life using this more general, `stoch_step` just without any randomness.

# In[465]:


def gol_transition(state, neighbors):
    if state == 1 and (neighbors[1] == 2 or neighbors[1] == 3):
        return 1
    if state == 0 and (neighbors[1] == 3):
        return 1
    return 0


# In[466]:


show_grid_anim(step_n(lambda g: stoch_step(g, gol_transition), random_grid(16, 0.2), 20))


# And we can now easily add some randomness to the Game of Life, this transition function is the same as before, but with a only a 99% chance (as opposed to 100%) of following the original rules.

# In[467]:


def stoch_gol_transition(state, neighbors):
    if state == 1 and (neighbors[1] == 2 or neighbors[1] == 3):
        return int(random() < 0.99) 
    if state == 0 and (neighbors[1] == 3):
        return int(random() < 0.99)
    return int(random() < 0.01)


# In[468]:


show_grid_anim(step_n(lambda g: stoch_step(g, stoch_gol_transition), random_grid(16, 0.1), 20))


# As you can see, we get some strange behavior. Now lets add more states to model disease transmissions! This is a very rudimentary example but the idea can easily be extended to more complex models.

# In[469]:


HEALTHY = 0
VACCINATED = 1
INFECTED = 2

def disease_transition(state, neighbors, pTransmission, pRecovery):
    if state == VACCINATED:
        return VACCINATED
    if state == HEALTHY:
        if random() < (neighbors[INFECTED] * pTransmission / sum(neighbors.values())):
            return INFECTED
        else:
            return HEALTHY
    if state == INFECTED:
        if random() < pRecovery:
            return VACCINATED
        else:
            return INFECTED
        
def disease_grid_init(N, pVaccinated, pInfected):
    pHealthy = 1 - pVaccinated - pInfected
    return np.random.choice(np.arange(3), size=N ** 2, p=[pHealthy, pVaccinated, pInfected]).reshape((N, N))
    


# In[470]:


def disease_anim(pTransmission, pRecovery, pVaccinated, pInfected, N, T):
    transition = partial(disease_transition, pTransmission=pTransmission, pRecovery=pRecovery)
    g = disease_grid_init(N, pVaccinated, pInfected)
    return grid_anim(step_n(lambda g: stoch_step(g, transition), g, T))


# In[471]:


anim_to_html(disease_anim(0.2, 0.05, 0.05, 0.01, 32, 100), fps=2)


# The above shows healthy cells in white, vaccinated cells in gray, and infected cells in black. Infected cells have a chance to recover and end up the same as the vaccinated ones (since they've developed an immunity after having the disease).
# 
# Once we have this framework for modeling, we can explore what would happen if we, say, have no vaccines and no good cure for a disease.

# In[472]:


anim_to_html(disease_anim(0.3, 0.01, 0.0, 0.01, 32, 100), fps=2)


# Bad news! this disease spread really fast and infected almost everyone!
# 
# What if we had 40% of people vaccinated, people have good hygiene (low transmission) and a decent treatment exists?

# In[473]:


anim_to_html(disease_anim(0.1, 0.05, 0.4, 0.01, 32, 20), fps=2)


# The disease quickly dies out! Yay!
# 
# As a parting video, we will go larger scale! Thanks for reading.

# In[474]:


anim_to_html(disease_anim(0.2, 0.05, 0.2, 0.01, 64, 100), fps=4)