<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Table of Contents
<p><div class="lev1 toc-item"><a href="#Preparation" data-toc-modified-id="Preparation-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Preparation</a></div><div class="lev1 toc-item"><a href="#Model-definition" data-toc-modified-id="Model-definition-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Model definition</a></div><div class="lev1 toc-item"><a href="#LRT" data-toc-modified-id="LRT-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>LRT</a></div><div class="lev2 toc-item"><a href="#get-optimal-parameter-values" data-toc-modified-id="get-optimal-parameter-values-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>get optimal parameter values</a></div><div class="lev2 toc-item"><a href="#get-bootstrap-replicates" data-toc-modified-id="get-bootstrap-replicates-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>get bootstrap replicates</a></div><div class="lev2 toc-item"><a href="#calculate-adjustment-for-D" data-toc-modified-id="calculate-adjustment-for-D-33"><span class="toc-item-num">3.3&nbsp;&nbsp;</span>calculate adjustment for D</a></div>
# Preparation
Step2:  Model definition
Step3:  LRT
get optimal parameter values
Step4:  This two-epoch model can be reduced to one-epoch model by either setting $Nery_2 = Nery_1$ and $Npar_2 = Npar_1$ or by setting $T_2 = 0$.
Step5:  get bootstrap replicates
Step6:  calculate adjustment for D
Step7:  From Coffman2016, suppl. mat.
<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
from glob import glob
import dill
import pandas as pd
# turn on floating point division by default, old behaviour via '//'
from __future__ import division
from itertools import repeat
def flatten(array):
    import numpy as np
    res  = []
    for el in array:
        if isinstance(el, (list, tuple, np.ndarray)):
            res.extend(flatten(el))
            continue
        res.append(el)
    return list(res)
%matplotlib inline
import pylab
pylab.rcParams['figure.figsize'] = [10, 8]
pylab.rcParams['font.size'] = 12
%%px --local
# load spectrum
sfs2d = dadi.Spectrum.from_file("EryPar.unfolded.sfs.dadi")
sfs2d = sfs2d.transpose()
sfs2d.pop_ids = ['ery', 'par']
sfs2d = sfs2d.fold()
ns = sfs2d.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.Plotting.plot_single_2d_sfs(sfs2d, vmin=1, cmap='jet')
pylab.savefig("2DSFS_folded.png")
# get number of segregating sites from SFS
sfs2d.S()
Explanation: Table of Contents
<p><div class="lev1 toc-item"><a href="#Preparation" data-toc-modified-id="Preparation-1"><span class="toc-item-num">1&nbsp;&nbsp;</span>Preparation</a></div><div class="lev1 toc-item"><a href="#Model-definition" data-toc-modified-id="Model-definition-2"><span class="toc-item-num">2&nbsp;&nbsp;</span>Model definition</a></div><div class="lev1 toc-item"><a href="#LRT" data-toc-modified-id="LRT-3"><span class="toc-item-num">3&nbsp;&nbsp;</span>LRT</a></div><div class="lev2 toc-item"><a href="#get-optimal-parameter-values" data-toc-modified-id="get-optimal-parameter-values-31"><span class="toc-item-num">3.1&nbsp;&nbsp;</span>get optimal parameter values</a></div><div class="lev2 toc-item"><a href="#get-bootstrap-replicates" data-toc-modified-id="get-bootstrap-replicates-32"><span class="toc-item-num">3.2&nbsp;&nbsp;</span>get bootstrap replicates</a></div><div class="lev2 toc-item"><a href="#calculate-adjustment-for-D" data-toc-modified-id="calculate-adjustment-for-D-33"><span class="toc-item-num">3.3&nbsp;&nbsp;</span>calculate adjustment for D</a></div>
# Preparation
End of explanation
def split_asym_mig_2epoch(params, ns, pts):
    
    params = (nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2)
    ns = (n1,n2)
    Split into two populations of specified size, with potentially asymmetric migration.
    The split coincides with a stepwise size change in the daughter populations. Then,
    have a second stepwise size change at some point in time after the split. This is
    enforced to happen at the same time for both populations. Migration is assumed to
    be the same during both epochs.
    nu1_1: pop size ratio of pop 1 after split (with respect to Na)
    nu2_1: pop size ratio of pop 2 after split (with respect to Na)
    T1: Time from split to second size change (in units of 2*Na generations)
    nu1_2: pop size ratio of pop 1 after second size change (with respect to Na)
    nu2_2: pop size ratio of pop 2 after second size change (with respect to Na)
    T2: time in past of second size change (in units of 2*Na generations)
    m1: Migration rate from ery into par (in units of 2*Na ind per generation)
    m2: Migration rate from par into ery (in units of 2*Na ind per generation)
    n1,n2: Sample sizes of resulting Spectrum
    pts: Number of grid points to use in integration.
    
    nu1_1,nu2_1,T1,nu1_2,nu2_2,T2,m1,m2 = params
    xx = dadi.Numerics.default_grid(pts)
    phi = dadi.PhiManip.phi_1D(xx)
    
    # split
    phi = dadi.PhiManip.phi_1D_to_2D(xx, phi) 
    
    # divergence with potentially asymmetric migration for time T1
    phi = dadi.Integration.two_pops(phi, xx, T1, nu1_1, nu2_1, m12=m2, m21=m1)
    
    # divergence with potentially asymmetric migration and different pop size for time T2
    phi = dadi.Integration.two_pops(phi, xx, T2, nu1_2, nu2_2, m12=m2, m21=m1)
    fs = dadi.Spectrum.from_phi(phi, ns, (xx,xx))
    return fs
cl[:].push(dict(split_asym_mig_2epoch=split_asym_mig_2epoch))
%%px --local
func_ex = dadi.Numerics.make_extrap_log_func(split_asym_mig_2epoch)
Explanation: Model definition
End of explanation
ar_split_asym_mig_2epoch = []
for filename in glob("OUT_2D_models/split_asym_mig_2epoch_[0-9]*dill"):
    ar_split_asym_mig_2epoch.append(dill.load(open(filename)))
l = 2*8+1
returned = [flatten(out)[:l] for out in ar_split_asym_mig_2epoch]
df = pd.DataFrame(data=returned, \
                 columns=['ery_1_0','par_1_0','T1_0','ery_2_0','par_2_0','T2_0','m1_0','m2_0', 'Nery_1_opt','Npar_1_opt','T1_opt','Nery_2_opt','Npar_2_opt','T2_opt','m1_opt','m2_opt','-logL'])
df.sort_values(by='-logL', ascending=True).iloc[:10,8:17]
# optimal parameter values for complex model
popt_c = np.array(df.sort_values(by='-logL', ascending=True).iloc[0, 8:16]) # take the 8th best parameter combination
popt_c
Explanation: LRT
get optimal parameter values
End of explanation
# optimal paramter values for simple model (1 epoch)
# note: Nery_2=Nery_1, Npar_2=Npar_1 and T2=0
popt_s = [1.24966921,  3.19164623,  1.42043464,  1.24966921,  3.19164623, 0.0, 0.08489757,  0.39827944]
Explanation: This two-epoch model can be reduced to one-epoch model by either setting $Nery_2 = Nery_1$ and $Npar_2 = Npar_1$ or by setting $T_2 = 0$.
End of explanation
# load bootstrapped 2D SFS
all_boot = [dadi.Spectrum.from_file("../SFS/bootstrap/2DSFS/{0:03d}.unfolded.2dsfs.dadi".format(i)).fold() for i in range(200)]
Explanation: get bootstrap replicates
End of explanation
# calculate adjustment for D evaluating at the *simple* model parameterisation
# specifying only T2 as fixed
adj_s = dadi.Godambe.LRT_adjust(func_ex, pts_l, all_boot, popt_s, sfs2d, nested_indices=[5], multinom=True)
adj_s
# calculate adjustment for D evaluating at the *complex* model parameterisation
# specifying only T2 as fixed
adj_c = dadi.Godambe.LRT_adjust(func_ex, pts_l, all_boot, popt_c, sfs2d, nested_indices=[5], multinom=True)
adj_c
Explanation: calculate adjustment for D
End of explanation
# calculate adjustment for D evaluating at the *simple* model parameterisation
# treating Nery_2, Npar_2 and T2 as nested
adj_s = dadi.Godambe.LRT_adjust(func_ex, pts_l, all_boot, popt_s, sfs2d, nested_indices=[3,4,5], multinom=True)
adj_s
# calculate adjustment for D evaluating at the *complex* model parameterisation
# treating Nery_2, Npar_2 and T2 as nested
adj_c = dadi.Godambe.LRT_adjust(func_ex, pts_l, all_boot, popt_c, sfs2d, nested_indices=[3,4,5], multinom=True)
adj_c
Explanation: From Coffman2016, suppl. mat.:
The two-epoch model can be marginalized down to the SNM model for an LRT by either setting η = 1 or T = 0. We found that the LRT adjustment performed well when treating both parameters as nested, so μ(θ) was evaluated with T = 0 and η = 1.
End of explanation
<END_TASK>