<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Table of Contents
<p><div class="lev1 toc-item"><a href="#Exponential-growth" data-toc-modified-id="Exponential-growth-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Exponential growth</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-12"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>PAR</a></div><div class="lev1 toc-item"><a href="#Two-epoch-model" data-toc-modified-id="Two-epoch-model-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Two epoch model</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-21"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-22"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>PAR</a></div><div class="lev1 toc-item"><a href="#Bottlegrowth" data-toc-modified-id="Bottlegrowth-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Bottlegrowth</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev3 toc-item"><a href="#Interpretation" data-toc-modified-id="Interpretation-311"><span class="toc-item-num">3.1.1&nbsp;&nbsp;</span>Interpretation</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>PAR</a></div><div class="lev3 toc-item"><a href="#Interpretation" data-toc-modified-id="Interpretation-321"><span class="toc-item-num">3.2.1&nbsp;&nbsp;</span>Interpretation</a></div><div class="lev1 toc-item"><a href="#Three-Epoch" data-toc-modified-id="Three-Epoch-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>Three Epoch</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-41"><span class="toc-item-num">4.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev3 toc-item"><a href="#Interpretation" data-toc-modified-id="Interpretation-411"><span class="toc-item-num">4.1.1&nbsp;&nbsp;</span>Interpretation</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-42"><span class="toc-item-num">4.2&nbsp;&nbsp;</span>PAR</a></div>
Step2:  Exponential growth
Step3:  ERY
Step4:  The time parameter is hitting the upper boundary that I set. The exponential growth model cannot be fit to the ery spectrum with Ludovics correction applied.
PAR
Step5:  The time parameter is hitting the upper boundary that I set. The exponential growth model can therefore not be fit to the PAR spectrum with Ludovic's correction applied.
Two epoch model
Step6:  ERY
Step7:  Still hitting upper boundary on time. The two epoch model cannot be fit to the ERY spectrum with Ludovic's correction applied.
PAR
Step8:  All hitting the upper boundary on the time parameter. The two epoch model cannot be fit to the PAR spectrum with Ludovic's correction applied.
Bottlegrowth
Step9:  ERY
Step10:  This looks like convergence.
Interpretation
Step11:  I think such a high effective population size is not realistic.
PAR
Step12:  This looks like convergence.
Interpretation
Step13:  An effective population size of 36 million is obviously to high. I therefore cannot regard this model fitting as successful.
Three Epoch
Step14:  ERY
Step15:  Reasonable convergence. Divergent parameter value combinations have the same likelihood. The optimal parameter values from the optimisation runs 16 and 10 (adjacent in the table) show that quite different demographic scenarios can have almost identical likelihood. This is not unusual.
Interpretation
Step16:  PAR
Step17:  There is no convergence. The three epoch model could not be fit to the PAR spectrum with Ludivic's correction.
<ASSISTANT_TASK:>
Python Code:
from ipyparallel import Client
cl = Client()
cl.ids
%%px --local
# run whole cell on all engines a well as in the local IPython session
import numpy as np
import sys
sys.path.insert(0, '/home/claudius/Downloads/dadi')
import dadi
%%px --local
# import 1D spectrum of ery on all engines:
fs_ery = dadi.Spectrum.from_file('ERY_modified.sfs')
# import 1D spectrum of par on all engines:
fs_par = dadi.Spectrum.from_file('PAR_modified.sfs')
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [12, 10]
pylab.rcParams['font.size'] = 14
pylab.plot(fs_ery, 'ro--', label='ery', markersize=12)
pylab.plot(fs_par, 'g>--', label='par', markersize=12)
pylab.grid()
pylab.xlabel('minor allele count')
pylab.ylabel('')
pylab.legend()
pylab.title('1D spectra - Ludivics correction applied')
def run_dadi(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    if perturb == True:
        p_init = dadi.Misc.perturb_params(p_init, fold=fold, 
                                      upper_bound=upper_bound, lower_bound=lower_bound)
        # note upper_bound and lower_bound variables are expected to be in the namespace of each engine
    # run optimisation of paramters
    popt = dadi_opt_func(p0=p_init, data=sfs, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=verbose, maxiter=maxiter, full_output=full_output, \
                                    fixed_params=fixed_params)
    # pickle to file
    import dill
    name = outname[:] # make copy of file name stub!
    for p in p_init:
        name += "_%.4f" % (p)
    with open(name + ".dill", "w") as fh:
        dill.dump((p_init, popt), fh)
    
    return p_init, popt
from glob import glob
import dill
from utility_functions import *
import pandas as pd
lbview = cl.load_balanced_view()
from itertools import repeat
Explanation: Table of Contents
<p><div class="lev1 toc-item"><a href="#Exponential-growth" data-toc-modified-id="Exponential-growth-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Exponential growth</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-11"><span class="toc-item-num">1.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-12"><span class="toc-item-num">1.2&nbsp;&nbsp;</span>PAR</a></div><div class="lev1 toc-item"><a href="#Two-epoch-model" data-toc-modified-id="Two-epoch-model-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Two epoch model</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-21"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-22"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>PAR</a></div><div class="lev1 toc-item"><a href="#Bottlegrowth" data-toc-modified-id="Bottlegrowth-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>Bottlegrowth</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev3 toc-item"><a href="#Interpretation" data-toc-modified-id="Interpretation-311"><span class="toc-item-num">3.1.1&nbsp;&nbsp;</span>Interpretation</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>PAR</a></div><div class="lev3 toc-item"><a href="#Interpretation" data-toc-modified-id="Interpretation-321"><span class="toc-item-num">3.2.1&nbsp;&nbsp;</span>Interpretation</a></div><div class="lev1 toc-item"><a href="#Three-Epoch" data-toc-modified-id="Three-Epoch-4"><span class="toc-item-num">4&nbsp;&nbsp;</span>Three Epoch</a></div><div class="lev2 toc-item"><a href="#ERY" data-toc-modified-id="ERY-41"><span class="toc-item-num">4.1&nbsp;&nbsp;</span>ERY</a></div><div class="lev3 toc-item"><a href="#Interpretation" data-toc-modified-id="Interpretation-411"><span class="toc-item-num">4.1.1&nbsp;&nbsp;</span>Interpretation</a></div><div class="lev2 toc-item"><a href="#PAR" data-toc-modified-id="PAR-42"><span class="toc-item-num">4.2&nbsp;&nbsp;</span>PAR</a></div>
End of explanation
dadi.Demographics1D.growth?
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(dadi.Demographics1D.growth)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 4]
lower_bound = [1e-4, 0]
Explanation: Exponential growth
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_ery
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/expgrowth" # set file name stub for opt. result files
fixed_params = None
# set starting values for perturbation
p0 = [1, 1]
#ar_ery = lbview.map(run_dadi, repeat(p0, 10))
ar_ery.get()
# set starting values for perturbation
p0 = [0.1, 0.1]
#ar_ery = lbview.map(run_dadi, repeat(p0, 10))
# set starting values for perturbation
p0 = [10, 0.1]
#ar_ery1 = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/expgrowth*dill"):
    ar_ery.append(dill.load(open(filename)))
get_flag_count(ar_ery, NM=True)
import pandas as pd
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 6] # increasing upper bound of time parameter
lower_bound = [1e-4, 0]
# set starting values for perturbation
p0 = [2, 1]
#ar_ery1 = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/expgrowth*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 8] # increasing upper bound of time parameter
lower_bound = [1e-4, 0]
# set starting values for perturbation
p0 = [2, 1]
ar_ery1 = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/expgrowth*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(15)
Explanation: ERY
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_par
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 300 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/PAR_expgrowth" # set file name stub for opt. result files
fixed_params = None
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 4]
lower_bound = [1e-4, 0]
# set starting values for perturbation
p0 = [1, 1]
ar_par = lbview.map(run_dadi, repeat(p0, 10))
ar_par = []
for filename in glob("OUT_1D_models/PAR_expgrowth*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(15)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 8]
lower_bound = [1e-4, 0]
# set starting values for perturbation
p0 = [1, 1]
ar_par = lbview.map(run_dadi, repeat(p0, 10))
ar_par = []
for filename in glob("OUT_1D_models/PAR_expgrowth*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(15)
Explanation: The time parameter is hitting the upper boundary that I set. The exponential growth model cannot be fit to the ery spectrum with Ludovics correction applied.
PAR
End of explanation
dadi.Demographics1D.two_epoch?
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(dadi.Demographics1D.two_epoch)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 4]
lower_bound = [1e-4, 0]
Explanation: The time parameter is hitting the upper boundary that I set. The exponential growth model can therefore not be fit to the PAR spectrum with Ludovic's correction applied.
Two epoch model
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_ery
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/ERY_twoEpoch" # set file name stub for opt. result files
fixed_params = None
# set starting values for perturbation
p0 = [1, 1]
ar_ery = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/ERY_twoEpoch*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 8]
lower_bound = [1e-4, 0]
# set starting values for perturbation
p0 = [1, 1]
ar_ery = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/ERY_twoEpoch*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: ERY
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_par
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/PAR_twoEpoch" # set file name stub for opt. result files
fixed_params = None
p0 = [1, 1]
ar_par = lbview.map(run_dadi, repeat(p0, 10))
ar_par = []
for filename in glob("OUT_1D_models/PAR_twoEpoch*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nu1_0', 'T_0', 'nu1_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: Still hitting upper boundary on time. The two epoch model cannot be fit to the ERY spectrum with Ludovic's correction applied.
PAR
End of explanation
dadi.Demographics1D.bottlegrowth?
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(dadi.Demographics1D.bottlegrowth)
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 1e4, 8]
lower_bound = [1e-4, 1e-4, 0]
Explanation: All hitting the upper boundary on the time parameter. The two epoch model cannot be fit to the PAR spectrum with Ludovic's correction applied.
Bottlegrowth
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_ery
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/ERY_bottlegrowth" # set file name stub for opt. result files
fixed_params = None
# set starting values for perturbation
p0 = [1, 1, 1]
#ar_ery = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/ERY_bottlegrowth*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'T_0', 'nuB_opt', 'nuF_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
# set starting values for perturbation
p0 = [55, 1.3, 1.5]
#ar_ery = lbview.map(run_dadi, repeat(p0, 10))
ar_ery = []
for filename in glob("OUT_1D_models/ERY_bottlegrowth*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'T_0', 'nuB_opt', 'nuF_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: ERY
End of explanation
popt = np.array( df.sort_values(by='-logL', ascending=True).iloc[0, 3:6] )
popt
# calculate best-fit model spectrum
model_spectrum = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum, fs_ery)
mu = 3e-9
L = fs_ery.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
Nref = theta/L/mu/4
Nref
print "At time {0:,} generations ago, the ERY population size instantaneously increased by almost 55-fold (to {1:,}).".format(int(popt[2]*2*Nref), int(popt[0]*Nref))
Explanation: This looks like convergence.
Interpretation
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_par
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/PAR_bottlegrowth" # set file name stub for opt. result files
fixed_params = None
%%px --local
# set lower and upper bounds to nu1 and T
upper_bound = [1e4, 1e4, 6]
lower_bound = [1e-4, 1e-4, 0]
# set starting values for perturbation
p0 = [1, 1, 1]
#ar_par = lbview.map(run_dadi, repeat(p0, 10))
ar_par = []
for filename in glob("OUT_1D_models/PAR_bottlegrowth*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'T_0', 'nuB_opt', 'nuF_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
cl[:]['maxiter'] = 100
# set starting values for perturbation
p0 = [100, 2, 1.2]
ar_par = lbview.map(run_dadi, repeat(p0, 10))
ar_par = []
for filename in glob("OUT_1D_models/PAR_bottlegrowth*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'T_0', 'nuB_opt', 'nuF_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(10)
Explanation: I think such a high effective population size is not realistic.
PAR
End of explanation
popt = np.array( df.sort_values(by='-logL', ascending=True).iloc[0, 3:6] )
popt
# calculate best-fit model spectrum
model_spectrum = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum, fs_par)
mu = 3e-9
L = fs_par.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
Nref = theta/L/mu/4
Nref
print "At time {0:,} generations ago, the PAR population size instantaneously increased by almost 124-fold (to {1:,}).".format(int(popt[2]*2*Nref), int(popt[0]*Nref))
Explanation: This looks like convergence.
Interpretation
End of explanation
dadi.Demographics1D.three_epoch?
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(dadi.Demographics1D.three_epoch)
%%px --local
# set lower and upper bounds to nuB, nuF, TB and TF
upper_bound = [1e4, 1e4, 6, 6]
lower_bound = [1e-4, 1e-4, 0, 0]
Explanation: An effective population size of 36 million is obviously to high. I therefore cannot regard this model fitting as successful.
Three Epoch
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_ery
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/ERY_threeEpoch" # set file name stub for opt. result files
fixed_params = None
# set starting values for perturbation
p0 = [10, 1, 1, 1]
#ar_ery = lbview.map(run_dadi, repeat(p0, 20))
ar_ery = []
for filename in glob("OUT_1D_models/ERY_threeEpoch*dill"):
    ar_ery.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_ery]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'TB_0', 'TF_0', 'nuB_opt', 'nuF_opt', 'TB_opt', 'TF_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: ERY
End of explanation
popt = np.array( df.sort_values(by='-logL', ascending=True).iloc[0, 4:8] )
popt
# calculate best-fit model spectrum
model_spectrum = func_ex(popt, ns, pts_l)
theta = dadi.Inference.optimal_sfs_scaling(model_spectrum, fs_ery)
mu = 3e-9
L = fs_ery.data.sum()
print "The optimal value of theta per site for the ancestral population is {0:.4f}.".format(theta/L)
Nref = theta/L/mu/4
Nref
print "At time {0:,} generations ago, the ERY population size instantaneously increased by almost 10-fold (to {1:,}).".format(int((popt[2]+popt[3])*2*Nref), int(popt[0]*Nref)),
print "It then kept this population constant for {0:,} generations.".format(int(popt[2]*2*Nref)),
print "At time {0:,} generations in the past, the ERY population then decreased to 1.3 fold of the ancient population size or {1:,}.".format(int(popt[3]*2*Nref), int(popt[1]*Nref))
Explanation: Reasonable convergence. Divergent parameter value combinations have the same likelihood. The optimal parameter values from the optimisation runs 16 and 10 (adjacent in the table) show that quite different demographic scenarios can have almost identical likelihood. This is not unusual.
Interpretation
End of explanation
%%px --local
# set up global variables on engines required for run_dadi function call
ns = fs_ery.sample_sizes # both populations have the same sample size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_par
perturb = True
fold = 2 # perturb randomly up to `fold` times 2-fold
maxiter = 100 # run a maximum of 100 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "MODIFIED_SPECTRA/OUT_1D_models/PAR_threeEpoch" # set file name stub for opt. result files
fixed_params = None
# set starting values for perturbation
p0 = [100, 2, 1, 1]
#ar_par = lbview.map(run_dadi, repeat(p0, 20))
ar_par = []
for filename in glob("OUT_1D_models/PAR_threeEpoch*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'TB_0', 'TF_0', 'nuB_opt', 'nuF_opt', 'TB_opt', 'TF_opt', '-logL'])
df.sort_values(by='-logL', ascending=True)
Explanation: PAR
End of explanation
%%px --local
fold = 1
maxiter = 300
# set starting values for perturbation
p0 = [20, 1e-2, 0.8, 1e-3]
ar_par = lbview.map(run_dadi, repeat(p0, 10))
ar_par = []
for filename in glob("OUT_1D_models/PAR_threeEpoch*dill"):
    ar_par.append(dill.load(open(filename)))
l = 2*len(p0)+1
# show all parameter combinations
returned = [flatten(out)[:l] for out in ar_par]
df = pd.DataFrame(data=returned, \
                  columns=['nuB_0', 'nuF_0', 'TB_0', 'TF_0', 'nuB_opt', 'nuF_opt', 'TB_opt', 'TF_opt', '-logL'])
df.sort_values(by='-logL', ascending=True).head(20)
Explanation: There is no convergence. The three epoch model could not be fit to the PAR spectrum with Ludivic's correction.
End of explanation
<END_TASK>