import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

Problem 3.5

Plot $\Delta_sp$ versus $p$, $p$ versus generation, and $\overline{w}$ versus $p$ for 200 generations, where $s$ = 0.1.
Do this for $h$ = -0.5, $h$ = 0, $h$ = 0.25, $h$ = 0.5, $h$ = 1, and $h$ = 1.5.

def calcp(h, pi):
    gen = 0
    s = 0.1
    
    p = np.empty(200)
    deltap = np.empty(200)
    w = np.empty(200)
    
    p[0] = pi
    
    while gen < 200:
        q = 1 - p[gen]
        w[gen] = 1 - 2*p[gen]*q*h*s - q*q*s
        deltap[gen] = (p[gen]*q*s*(p[gen]*h + q*(1 - h)))/w[gen]
        if gen != 199:
            p[gen + 1] = p[gen] + deltap[gen]
        gen = gen + 1
    
    return p, deltap, w

p, deltap, w = calcp(-0.5, 0.01)
h = [0, 0.25, 0.5, 1, 1.5]
for x in h:
    px, deltapx, wx = calcp(x, 0.01)
    p = np.vstack([p, px])
    deltap = np.vstack([deltap, deltapx])
    w = np.vstack([w, wx])

fig, ax = plt.subplots(facecolor='w', edgecolor='black')
plt.plot(p[0], color = 'red', linewidth='1', label='Overdominant')
plt.plot(p[1], color = 'orange', linewidth='1', label='Recessive')
plt.plot(p[2], color = 'green', linewidth='1', label='Incomplete dominance')
plt.plot(p[3], color = 'blue', linewidth='1', label='Additive')
plt.plot(p[4], color = 'purple', linewidth='1', label='Dominant')
plt.plot(p[5], color = 'black', linewidth='1', label='Underdominant')
ax.set_ymargin(0.0)
ax.set_xmargin(0.0)
ax.set_ylabel('$p$');
ax.set_xlabel('Generation');
ax.set_ylim(0, 1);
ax.set_xlim(0, 200);
plt.legend();

fig, ax = plt.subplots(facecolor='w', edgecolor='black')
plt.plot(p[0], deltap[0], color = 'red', linewidth='1', label='Overdominant')
plt.plot(p[1], deltap[1], color = 'orange', linewidth='1', label='Recessive')
plt.plot(p[2], deltap[2], color = 'green', linewidth='1', label='Incomplete dominance')
plt.plot(p[3], deltap[3], color = 'blue', linewidth='1', label='Additive')
plt.plot(p[4], deltap[4], color = 'purple', linewidth='1', label='Dominant')
plt.plot(p[5], deltap[5], color = 'black', linewidth='1', label='Underdominant')
ax.set_ymargin(0.0)
ax.set_xmargin(0.0)
ax.set_ylabel('$\Delta_ps$');
ax.set_xlabel('$p$');
ax.set_ylim(-0.05, 0.05);
ax.set_xlim(0, 1);
plt.legend();

fig, ax = plt.subplots(facecolor='w', edgecolor='black')
plt.plot(p[0], w[0], color = 'red', linewidth='1', label='Overdominant')
plt.plot(p[1], w[1], color = 'orange', linewidth='1', label='Recessive')
plt.plot(p[2], w[2], color = 'green', linewidth='1', label='Incomplete dominance')
plt.plot(p[3], w[3], color = 'blue', linewidth='1', label='Additive')
plt.plot(p[4], w[4], color = 'purple', linewidth='1', label='Dominant')
plt.plot(p[5], w[5], color = 'black', linewidth='1', label='Underdominant')
ax.set_ymargin(0.0)
ax.set_xmargin(0.0)
ax.set_ylabel('$\overline{w}$');
ax.set_xlabel('$p$');
#ax.set_ylim(0.9, 1);
ax.set_xlim(0, 1);
plt.legend();

Same graphs, but now the initial frequency is 0.9

p, deltap, w = calcp(-0.5, 0.9)
h = [0, 0.25, 0.5, 1, 1.5]
for x in h:
    px, deltapx, wx = calcp(x, 0.01)
    p = np.vstack([p, px])
    deltap = np.vstack([deltap, deltapx])
    w = np.vstack([w, wx])

fig, ax = plt.subplots(facecolor='w', edgecolor='black')
plt.plot(p[0], color = 'red', linewidth='1', label='Overdominant')
plt.plot(p[1], color = 'orange', linewidth='1', label='Recessive')
plt.plot(p[2], color = 'green', linewidth='1', label='Incomplete dominance')
plt.plot(p[3], color = 'blue', linewidth='1', label='Additive')
plt.plot(p[4], color = 'purple', linewidth='1', label='Dominant')
plt.plot(p[5], color = 'black', linewidth='1', label='Underdominant')
ax.set_ymargin(0.0)
ax.set_xmargin(0.0)
ax.set_ylabel('$p$');
ax.set_xlabel('Generation');
ax.set_ylim(0, 1);
ax.set_xlim(0, 200);
plt.legend();

fig, ax = plt.subplots(facecolor='w', edgecolor='black')
plt.plot(p[0], deltap[0], color = 'red', linewidth='1', label='Overdominant')
plt.plot(p[1], deltap[1], color = 'orange', linewidth='1', label='Recessive')
plt.plot(p[2], deltap[2], color = 'green', linewidth='1', label='Incomplete dominance')
plt.plot(p[3], deltap[3], color = 'blue', linewidth='1', label='Additive')
plt.plot(p[4], deltap[4], color = 'purple', linewidth='1', label='Dominant')
plt.plot(p[5], deltap[5], color = 'black', linewidth='1', label='Underdominant')
ax.set_ymargin(0.0)
ax.set_xmargin(0.0)
ax.set_ylabel('$\Delta_ps$');
ax.set_xlabel('$p$');
ax.set_ylim(-0.05, 0.05);
ax.set_xlim(0, 1);
plt.legend();

fig, ax = plt.subplots(facecolor='w', edgecolor='black')
plt.plot(p[0], w[0], color = 'red', linewidth='1', label='Overdominant')
plt.plot(p[1], w[1], color = 'orange', linewidth='1', label='Recessive')
plt.plot(p[2], w[2], color = 'green', linewidth='1', label='Incomplete dominance')
plt.plot(p[3], w[3], color = 'blue', linewidth='1', label='Additive')
plt.plot(p[4], w[4], color = 'purple', linewidth='1', label='Dominant')
plt.plot(p[5], w[5], color = 'black', linewidth='1', label='Underdominant')
ax.set_ymargin(0.0)
ax.set_xmargin(0.0)
ax.set_ylabel('$\overline{w}$');
ax.set_xlabel('$p$');
#ax.set_ylim(0.9, 2);
ax.set_xlim(0, 1);
plt.legend();

import os
os.system('jupyter nbconvert --to html GillespieCh3.ipynb');