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

# In[1]:


get_ipython().run_line_magic('matplotlib', 'inline')
import matplotlib.pyplot as plt
import pandas as pd
import scipy


# In[14]:


mu = 1.25e-8
gen = 30
afrDat = pd.read_csv("/home/training/share/MSMC-tutorial-files/results/AFR.msmc2.final.txt", delim_whitespace=True)
eurDat = pd.read_csv("/home/training/share/MSMC-tutorial-files/results/EUR.msmc2.final.txt", delim_whitespace=True)
plt.step(afrDat["left_time_boundary"]/mu*gen, (1/afrDat["lambda"])/(2*mu), label="AFR")
plt.step(eurDat["left_time_boundary"]/mu*gen, (1/eurDat["lambda"])/(2*mu), label="EUR")
plt.ylim(0,40000)
plt.xlabel("years ago");
plt.ylabel("effective population size");
plt.gca().set_xscale('log')
plt.legend()


# In[20]:


mu = 1.25e-8
gen = 30
crossPopDat = pd.read_csv("/home/training/share/MSMC-tutorial-files/results/EUR_AFR.combined.msmc2.final.txt", delim_whitespace=True)
plt.step(crossPopDat["left_time_boundary"]/mu*gen, 2 * crossPopDat["lambda_01"] / (crossPopDat["lambda_00"] + crossPopDat["lambda_11"]))
plt.xlim(1000,500000);
plt.ylim(0, 1.2)
plt.xlabel("years ago");
plt.ylabel("relative cross coalescence rate");
plt.gca().set_xscale('log')


# In[22]:


[1,2,3]


# In[ ]: