

In [1]:

    
import sys
sys.path.append('../..')
from IPython.display import Image



In [3]:

    
import functools
import math
import os
import matplotlib.pyplot as plt
%matplotlib inline
import plot
import trout
pref = 'trt'
import seaborn as sns
sns.set_style('white')









    



/core/software/conda/lib/python3.5/site-packages/matplotlib/__init__.py:1350: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)



In [4]:

    
nsnps = [100, 200, 400]

Figure 1

Violin plot of point estimators. With and without logquad correction. Red dot: 10 and 90% of CONFIDENCE intervals

Nb = 50



In [5]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_lt_comp(case, 50, "Newb", "None")









    



---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-5-d6cea99fdfa0> in <module>()
----> 1 case = trout.load_file(pref, 40, mydir='../..')
      2 fig = plot.do_lt_comp(case, 50, "Newb", "None")

/scratch/1/tiago_antao/AgeStructureNe/trout.py in load_file(pref, cut, mydir)
     71         f = open("%s/output/%s.txt" % (mydir, pref))
     72     else:
---> 73         f = open("%s/output/%s-%d.txt" % (mydir, pref, cut))
     74     case = {}
     75     f.readline()

FileNotFoundError: [Errno 2] No such file or directory: '../../output/trt-40.txt'



In [5]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_lt_comp(case, 50, "Newb", "NbNe")

Nb = 100



In [6]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_lt_comp(case, 100, "Newb", "None")



In [7]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_lt_comp(case, 100, "Newb", "NbNe")

Figure 2

Different sample sizes

NB = 25



In [4]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 86, nsnps, 'allsnps-', 'None')



In [5]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 86, nsnps, 'allsnps-', 'NbNe')

Nb = 50



In [6]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 174, nsnps, 'allsnps-', 'None')



In [7]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 174, nsnps, 'allsnps-', 'NbNe')

Nb = 100



In [8]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 353, nsnps, 'allsnps-', 'None')



In [9]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 353, nsnps, 'allsnps-', 'NbNe')

Nb = 200



In [10]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 713, nsnps, 'allsnps-', 'None')



In [11]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_nb(case, 'Newb', 713, nsnps, 'allsnps-', 'NbNe')

Figure 3

Nb = 50



In [4]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_cohort(case, "bulltrout", 174, 50, 'None')



In [5]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_cohort(case, "bulltrout", 174, 50, 'NbNe')

Nb = 100



In [6]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_cohort(case, "bulltrout", 353, 50, 'None')



In [7]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_cohort(case, "bulltrout", 353, 50, 'NbNe')

Figure 4



In [ ]:

    
#Cannot do with 200, too little sibs



In [25]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_rel(case, "bulltrout", 174, 50, 'NbNe')



In [26]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.do_rel(case, "bulltrout", 86, 50, 'NbNe')

Figure 5

100 SNPs



In [4]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 86, [100], [25, 50, 100])



In [7]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 174, [100], [25, 50, 100])



In [8]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 353, [100], [25, 50, 100])



In [9]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 713, [100], [25, 50, 100])

200 SNPs



In [10]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 86, [200], [25, 50, 100])



In [11]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 174, [200], [25, 50, 100])



In [12]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 353, [200], [25, 50, 100])



In [13]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 713, [200], [25, 50, 100])

400 SNPs



In [14]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 86, [400], [25, 50, 100])



In [15]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 174, [400], [25, 50, 100])



In [16]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 353, [400], [25, 50, 100])



In [17]:

    
case = trout.load_file(pref, 40, mydir='../..')
fig = plot.compare_correction_ci(case, 'bulltrout', 713, [400], [25, 50, 100])

Figure 6



In [5]:

    
case = trout.load_file(pref, 40, mydir='../..')
linear = [("bullt2", 606), ("bullt2", 1214), ("bullt2", 1822),
          ("bullt2", 2433), ("bullt2", 3040),
          ("bullt2", 4565), ("bullt2", 6090),]
fig = plot.do_nb_linear(case, linear, "NbNe", functools.partial(trout.correct_ci))



In [6]:

    
fig = plot.do_nb_linear(case, linear, "NbNe", functools.partial(trout.correct_ci, jcorr=0.2))

Figure S1



In [6]:

    
plot.do_hz_comp(pref, '../..', "bulltrout", 174)



In [7]:

    
plot.do_hz_comp(pref, '../..', "bulltrout", 353)

Pcrit



In [20]:

    
case = trout.load_file(pref, 25, mydir='../..')
plot.do_pcrit(case, "bulltrout", 353, True)



In [21]:

    
def get_limits(nb):
    top = nb + 2 * math.sqrt(nb / 2)
    bottom = nb - 2 * math.sqrt(nb / 2)
    return top, bottom

ratios = []
nbs = [10, 20, 30, 40, 50, 75, 100, 200]
for nb in nbs:
    top, bottom = get_limits(nb)
    print(nb, top)
    ratios.append(top / nb)
fig, ax = plt.subplots()
ax.plot(nbs, ratios, '.-')









    



10 14.47213595499958
20 26.32455532033676
30 37.74596669241483
40 48.94427190999916
50 60.0
75 87.24744871391589
100 114.14213562373095
200 220.0






    Out[21]:





[<matplotlib.lines.Line2D at 0x7f1d4b477e10>]

Tables



In [ ]:

    
case = trout.load_file(pref, 40, mydir='../..')
cis = [("WFrog", "wfrog", 148), ("WFrog", "wfrog", 297),
       ("WFrog", "wfrog", 597),
       ("Mosq", "mosquito", 102), ("Mosq", "mosquito", 167),
       ("Mosq", "mosquito", 210), ("Mosq", "mosquito", 428),
       ("SynSweed", "synseaweed", 81), ("SynSweed", "synseaweed", 164),
       ("SynSweed", "synseaweed", 207),
       ("BT-Std", "bulltrout", 86), ("BT-Std", "bulltrout", 174),
       ("BT-Std", "bulltrout", 353), ("BT-Std", "bulltrout", 713),
       ("WCT-S", "shepard", 513), ("WCT-S", "shepard", 1030),
       ("WCT-F", "fraley", 637), ("WCT-F", "fraley", 1278),
       ("BT-Long", "bullt2", 606), ("BT-Long", "bullt2", 1214),
       ("BT-Long", "bullt2", 1822),
       ("BT-Long", "bullt2", 2433), ("BT-Long", "bullt2", 3040),
       ("BT-Long", "bullt2", 4565), ("BT-Long", "bullt2", 6090),
       ("BT-Pred", "bullpred", 189), ("BT-Pred", "bullpred", 381),
       ("BT-Pred", "bullpred", 765)]
medians = plot.do_table_ci(trout.Nbs, case, '../../output/', cis, 100, 50)
##plot.do_table_ci(case, '../../output/', cis, 100, 50, 20)
#plot.do_table_ci(case, '../../output/', cis, 100, 25)
#plot.do_table_ci(case, '../../output/', cis, 200, 50)



In [ ]:

    
corr_nbs = {}
for (model, N0), val in trout.Nbs.items():
    nb, real = medians[(model, N0)]
    corr_nbs[(model, N0)] = real



In [ ]:

    
case = trout.load_file(pref, 40, mydir='../..')
corr_medians = plot.do_table_ci(corr_nbs, case, '/tmp/', cis, 100, 50)



In [ ]:

    
case = trout.load_file(pref, 40, mydir='../..')
medians = plot.do_table_ci(trout.Nbs, case, '/tmp/', cis, 100, 50, flex_nb=True)



In [ ]: