fixed a minor bug to match update to stdpopsim repo

n_t/plots.py

@@ -1,40 +1,35 @@
 """
 Code for generating plots.
 """
-
-import numpy as np
+import pandas
 import matplotlib
+from matplotlib import pyplot as plt
 # Force matplotlib to not use any Xwindows backend.
 matplotlib.use('Agg')
-from matplotlib import pyplot as plt
-import pandas
 
-def plot_stairway_Ne_estimate(infile, outfile):
 
+def plot_stairway_Ne_estimate(infile, outfile):
+    """
+    figure of N(t) for single run of stairwayplot
+    """
     nt = pandas.read_csv(infile, sep="\t", skiprows=5)
     nt = nt[nt['year'] > 10]
     f, ax = plt.subplots(figsize=(7, 7))
     ax.set(xscale="log", yscale="log")
-    # JK: turned off these limits as they weren't working for me.
-           # ylim=(np.min(nt['Ne_2.5%']-np.std(nt['Ne_2.5%'])),
-           #       np.max(nt['Ne_97.5%'])+np.std(nt['Ne_97.5%'])))
     ax.plot(nt['year'], nt['Ne_median'], c="red")
     ax.plot(nt['year'], nt['Ne_2.5%'], c='grey')
     ax.plot(nt['year'], nt['Ne_97.5%'], c='grey')
     f.savefig(outfile, bbox_inches='tight')
 
-def plot_compound_Ne_estimate(infiles, outfile):
 
+def plot_compound_Ne_estimate(infiles, outfile):
+    """
+    figure of N(t) for multiple runs of stairwayplot
+    """
     f, ax = plt.subplots(figsize=(7, 7))
     ax.set(xscale="log", yscale="log")
     for infile in infiles:
         nt = pandas.read_csv(infile, sep="\t", skiprows=5)
         nt = nt[nt['year'] > 10]
-        # JK: turned off these limits as they weren't working for me.
-               # ylim=(np.min(nt['Ne_2.5%']-np.std(nt['Ne_2.5%'])),
-               #       np.max(nt['Ne_97.5%'])+np.std(nt['Ne_97.5%'])))
         ax.plot(nt['year'], nt['Ne_median'], c="red")
-#        ax.plot(nt['year'], nt['Ne_2.5%'], c='grey')
-#        ax.plot(nt['year'], nt['Ne_97.5%'], c='grey')
-
     f.savefig(outfile, bbox_inches='tight')

n_t/simulations.py

@@ -4,12 +4,10 @@
 
 import msprime
 from stdpopsim import homo_sapiens
-import sys
+
 
 def homo_sapiens_Gutenkunst(path, seed, chrmStr, sample_size=20):
     chrom = homo_sapiens.genome.chromosomes[chrmStr]
-    recomb_map = chrom.recombination_map()
-
     model = homo_sapiens.GutenkunstThreePopOutOfAfrica()
     # model.debug()
 
@@ -18,7 +16,7 @@ def homo_sapiens_Gutenkunst(path, seed, chrmStr, sample_size=20):
     ts = msprime.simulate(
         samples=samples,
         recombination_map=chrom.recombination_map(),
-        mutation_rate=chrom.mean_mutation_rate,
+        mutation_rate=chrom.default_mutation_rate,
         random_seed=seed,
         **model.asdict())
     ts.dump(path)