Commit e82a42cb authored by Ryan Gutenkunst's avatar Ryan Gutenkunst
Browse files

Add 5D admixture methods and tests

parent ec968f8e
......@@ -317,24 +317,44 @@ def _four_pop_admixture_intermediates(phi_4D, f1,f2,f3, xx,yy,zz,aa,bb):
+ (1-f1-f2-f3)*aa[nuax,nuax,nuax,:]
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _admixture_intermediates(phi_4D, ad_w, aa)
= _admixture_intermediates(phi_4D, ad_w, bb)
return lower_w_index, upper_w_index, frac_lower, frac_upper, norm
def phi_2D_to_3D_admix(phi, f, xx,yy,zz):
def _five_pop_admixture_intermediates(phi_5D, f1,f2,f3,f4, xx,yy,zz,aa,bb,cc):
"""
Intermediate results used when splitting a new population out from a 4D phi.
ww: frequency mapping for the new population.
"""
# For each point x,y,z,a,b in phi, this is the corresponding frequency c that
# SNPs with frequency x,y,z,a,b in populations 1,2,3,4,5 would map to.
if f1 + f2 + f3 + f4 > 1:
raise ValueError('Admixture proportions (f1=%f, f2 = %f, f3=%f, f4=%f) are '
'non-sensible.' % (f1, f2, f3, f4))
ad_w = f1*xx[:,nuax,nuax,nuax,nuax] + f2*yy[nuax,:,nuax,nuax,nuax] + f3*zz[nuax,nuax,:,nuax,nuax]\
+ f4*aa[nuax,nuax,nuax,:,nuax] + (1-f1-f2-f3-f4)*bb[nuax,nuax,nuax,nuax,:]
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _admixture_intermediates(phi_5D, ad_w, cc)
return lower_w_index, upper_w_index, frac_lower, frac_upper, norm
def phi_2D_to_3D_admix(phi, f1, xx,yy,zz):
"""
Create population 3 admixed from populations 1 and 2.
Returns a 3D sfs of shape (len(xx),len(yy),len(zz))
phi: phi corresponding to original 2 populations
f: Fraction of population 3 derived from population 1. (A fraction 1-f
f1: Fraction of population 3 derived from population 1. (A fraction 1-f1
will be derived from population 2.)
xx,yy: Mapping of points in phi to frequencies in populations 1 and 2.
zz: Frequency mapping that will be used along population 3 axis.
"""
lower_z_index, upper_z_index, frac_lower, frac_upper, norm \
= _two_pop_admixture_intermediates(phi, f, xx,yy,zz)
= _two_pop_admixture_intermediates(phi, f1, xx,yy,zz)
# Assemble our result.
......@@ -693,6 +713,166 @@ def phi_4D_admix_into_2(phi, f1,f3,f4, xx,yy,zz,aa):
return phi
def phi_5D_admix_into_1(phi, f2,f3,f4,f5, xx,yy,zz,aa,bb):
"""
Admix populations 2, 3, 4, and 5 into population 1.
Alters phi in place and returns the new version.
phi: phi corresponding to original 5 populations.
f2: Fraction of updated population 1 to be derived from population 2.
f3: Fraction of updated population 1 to be derived from population 3.
f4: Fraction of updated population 1 to be derived from population 4.
f5: Fraction of updated population 1 to be derived from population 5.
A fraction (1-f2-f3-f4-f5) will be derived from the original pop 1.
xx,yy,zz,aa,bb: Mapping of points in phi to frequencies in populations 1,2,3,4, and 5.
"""
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _five_pop_admixture_intermediates(phi, 1-f2-f3-f4-f5,f2,f3,f4, xx,yy,zz,aa,bb, xx)
lower_cont = frac_lower * norm
upper_cont = frac_upper * norm
idx_i = numpy.arange(phi.shape[0])
for jj in range(phi.shape[1]):
for kk in range(phi.shape[2]):
for ll in range(phi.shape[3]):
for mm in range(phi.shape[4]):
phi_int = numpy.zeros((phi.shape[0], phi.shape[0]))
phi_int[idx_i, lower_w_index[:,jj,kk,ll,mm]] = lower_cont[:,jj,kk,ll,mm]
phi_int[idx_i, upper_w_index[:,jj,kk,ll,mm]] = upper_cont[:,jj,kk,ll,mm]
phi[:,jj,kk,ll,mm] = Numerics.trapz(phi_int, xx, axis=0)
return phi
def phi_5D_admix_into_2(phi, f1,f3,f4,f5, xx,yy,zz,aa,bb):
"""
Admix populations 1, 3, 4, and 5 into population 2.
Alters phi in place and returns the new version.
phi: phi corresponding to original 5 populations.
f1: Fraction of updated population 2 to be derived from population 1.
f3: Fraction of updated population 2 to be derived from population 3.
f4: Fraction of updated population 2 to be derived from population 4.
f5: Fraction of updated population 2 to be derived from population 5.
A fraction (1-f1-f3-f4-f5) will be derived from the original pop 2.
xx,yy,zz,aa,bb: Mapping of points in phi to frequencies in populations 1,2,3,4, and 5.
"""
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _five_pop_admixture_intermediates(phi, f1, 1-f1-f3-f4-f5,f3,f4, xx,yy,zz,aa,bb, xx)
lower_cont = frac_lower * norm
upper_cont = frac_upper * norm
idx_j = numpy.arange(phi.shape[1])
for ii in range(phi.shape[0]):
for kk in range(phi.shape[2]):
for ll in range(phi.shape[3]):
for mm in range(phi.shape[4]):
phi_int = numpy.zeros((phi.shape[1], phi.shape[1]))
phi_int[idx_j, lower_w_index[ii,:,kk,ll,mm]] = lower_cont[ii,:,kk,ll,mm]
phi_int[idx_j, upper_w_index[ii,:,kk,ll,mm]] = upper_cont[ii,:,kk,ll,mm]
phi[ii,:,kk,ll,mm] = Numerics.trapz(phi_int, xx, axis=0)
return phi
def phi_5D_admix_into_3(phi, f1,f2,f4,f5, xx,yy,zz,aa,bb):
"""
Admix populations 1, 2, 4, and 5 into population 3.
Alters phi in place and returns the new version.
phi: phi corresponding to original 5 populations.
f1: Fraction of updated population 3 to be derived from population 1.
f2: Fraction of updated population 3 to be derived from population 2.
f4: Fraction of updated population 3 to be derived from population 4.
f5: Fraction of updated population 3 to be derived from population 5.
A fraction (1-f1-f2-f4-f5) will be derived from the original pop 3.
xx,yy,zz,aa,bb: Mapping of points in phi to frequencies in populations 1,2,3,4, and 5.
"""
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _five_pop_admixture_intermediates(phi, f1, f2, 1-f1-f2-f4-f5,f4, xx,yy,zz,aa,bb, xx)
lower_cont = frac_lower * norm
upper_cont = frac_upper * norm
idx_k = numpy.arange(phi.shape[2])
for ii in range(phi.shape[0]):
for jj in range(phi.shape[1]):
for ll in range(phi.shape[3]):
for mm in range(phi.shape[4]):
phi_int = numpy.zeros((phi.shape[2], phi.shape[2]))
phi_int[idx_k, lower_w_index[ii,jj,:,ll,mm]] = lower_cont[ii,jj,:,ll,mm]
phi_int[idx_k, upper_w_index[ii,jj,:,ll,mm]] = upper_cont[ii,jj,:,ll,mm]
phi[ii,jj,:,ll,mm] = Numerics.trapz(phi_int, xx, axis=0)
return phi
def phi_5D_admix_into_4(phi, f1,f2,f3,f5, xx,yy,zz,aa,bb):
"""
Admix populations 1, 2, 3, and 5 into population 4.
Alters phi in place and returns the new version.
phi: phi corresponding to original 5 populations.
f1: Fraction of updated population 4 to be derived from population 1.
f2: Fraction of updated population 4 to be derived from population 2.
f3: Fraction of updated population 4 to be derived from population 3.
f5: Fraction of updated population 4 to be derived from population 5.
A fraction (1-f1-f2-f3-f5) will be derived from the original pop 4.
xx,yy,zz,aa,bb: Mapping of points in phi to frequencies in populations 1,2,3,4, and 5.
"""
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _five_pop_admixture_intermediates(phi, f1, f2, f3, 1-f1-f2-f3-f5, xx,yy,zz,aa,bb, xx)
lower_cont = frac_lower * norm
upper_cont = frac_upper * norm
idx_l = numpy.arange(phi.shape[3])
for ii in range(phi.shape[0]):
for jj in range(phi.shape[1]):
for kk in range(phi.shape[2]):
for mm in range(phi.shape[4]):
phi_int = numpy.zeros((phi.shape[3], phi.shape[3]))
phi_int[idx_l, lower_w_index[ii,jj,kk,:,mm]] = lower_cont[ii,jj,kk,:,mm]
phi_int[idx_l, upper_w_index[ii,jj,kk,:,mm]] = upper_cont[ii,jj,kk,:,mm]
phi[ii,jj,kk,:,mm] = Numerics.trapz(phi_int, xx, axis=0)
return phi
def phi_5D_admix_into_5(phi, f1,f2,f3,f4, xx,yy,zz,aa,bb):
"""
Admix populations 1, 2, 3, and 4 into population 5.
Alters phi in place and returns the new version.
phi: phi corresponding to original 5 populations.
f1: Fraction of updated population 5 to be derived from population 1.
f2: Fraction of updated population 5 to be derived from population 2.
f3: Fraction of updated population 5 to be derived from population 3.
f4: Fraction of updated population 5 to be derived from population 3.
A fraction (1-f1-f2-f3-f4) will be derived from the original pop 5.
xx,yy,zz,aa,bb: Mapping of points in phi to frequencies in populations 1,2,3,4, and 5.
"""
lower_w_index, upper_w_index, frac_lower, frac_upper, norm \
= _five_pop_admixture_intermediates(phi, f1, f2, f3, f4, xx,yy,zz,aa,bb, xx)
lower_cont = frac_lower * norm
upper_cont = frac_upper * norm
idx_m = numpy.arange(phi.shape[4])
for ii in range(phi.shape[0]):
for jj in range(phi.shape[1]):
for kk in range(phi.shape[2]):
for ll in range(phi.shape[3]):
phi_int = numpy.zeros((phi.shape[4], phi.shape[4]))
phi_int[idx_m, lower_w_index[ii,jj,kk,ll,:]] = lower_cont[ii,jj,kk,ll,:]
phi_int[idx_m, upper_w_index[ii,jj,kk,ll,:]] = upper_cont[ii,jj,kk,ll,:]
phi[ii,jj,kk,ll,:] = Numerics.trapz(phi_int, xx, axis=0)
return phi
def remove_pop(phi, xx, popnum):
"""
Remove a population from phi.
......
......@@ -31,26 +31,26 @@ class Int5DTestCase(unittest.TestCase):
"""
Test simple SNM integration.
"""
pts = 30
pts = 5
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
fs1 = dadi.Spectrum.from_phi(phi, [5], (xx,))
ref = dadi.Integration.one_pop(phi, xx, T=0.1)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.PhiManip.phi_2D_to_3D_admix(phi, 0, xx,xx,xx)
phi = dadi.PhiManip.phi_3D_to_4D(phi, 0,0, xx,xx,xx,xx)
phi = dadi.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
# No demography, so all pops should still be SNM
phi = dadi.Integration.five_pops(phi, xx, T=0.1)
fs5 = dadi.Spectrum.from_phi(phi, [5,5,5,5,5], (xx,xx,xx,xx,xx))
# Loose tolerance of this comparison is okay. We're looking
# for gross violation due to bugs.
# Check all marginal phi, which should be identical except for 0,1 entries.
for ii in range(5):
tomarg = list(range(5))
tomarg.remove(ii)
fsm = fs5.marginalize(tomarg)
self.assertTrue(np.allclose(fs1, fsm, rtol=5e-2, atol=5e-2))
toremove = list(range(5))
toremove.remove(ii)
# Remove all the populations except the one we care about
test = phi
for pop in toremove[::-1]:
test = dadi.PhiManip.remove_pop(test, xx, pop)
print(np.allclose(ref[1:-1],test[1:-1]))
def test_integration_nomig(self):
nu1 = lambda t: 0.5 + 50*t
......@@ -149,86 +149,101 @@ class Int5DTestCase(unittest.TestCase):
test = test_func_marg(pin, None, 2, [16, 18, 20])
self.assertTrue(np.allclose(ref, test, atol=1e-3))
#def test_admix_props(self):
# """
# Test admix_props in from_phi
# """
# pts = 20
# xx = dadi.Numerics.default_grid(pts)
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
# fs2 = dadi.Spectrum.from_phi(phi, [5,5], (xx,xx), admix_props=((0.3,0.7),(0.9,0.1)))
#
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
# phi = dadi.PhiManip.phi_3D_to_4D(phi, 0, 0, xx,xx,xx,xx)
# phi = dadi.Integration.four_pops(phi, xx, T=0.1, nu3=0.5, nu4=10, m34=2, m43=0.5, gamma3=-1, gamma4=1, h3=0.2, h4=0.7)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx),
# admix_props=((1,0,0,0),(0,1,0,0),(0,0,0.3,0.7),(0,0,0.9,0.1)))
# fsm = fs4.marginalize((0,1))
# self.assertTrue(np.allclose(fs2, fsm, rtol=1e-2, atol=1e-2))
#def test_admixture(self):
# """
# Test phi_4D_admix_into_4
# """
# pts = 20
# xx = dadi.Numerics.default_grid(pts)
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
# phi = dadi.PhiManip.phi_2D_admix_1_into_2(phi,0.8,xx,xx)
# fs2 = dadi.Spectrum.from_phi(phi, [5,5], (xx,xx))
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
# phi = dadi.PhiManip.phi_3D_to_4D(phi, 0, 0, xx,xx,xx,xx)
# phi = dadi.Integration.four_pops(phi, xx, T=0.1, nu3=0.5, nu4=10, m34=2, m43=0.5, gamma3=-1, gamma4=1, h3=0.2, h4=0.7)
# phi = dadi.PhiManip.phi_4D_admix_into_4(phi,0,0,0.8,xx,xx,xx,xx)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((0,1))
# self.assertTrue(np.allclose(fs2, fsm, rtol=1e-2, atol=1e-2))
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
# phi = dadi.PhiManip.phi_3D_to_4D(phi, 0, 0, xx,xx,xx,xx)
# phi = dadi.Integration.four_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.7)
# phi = dadi.PhiManip.phi_4D_admix_into_2(phi,0.8,0,0,xx,xx,xx,xx)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((2,3))
# self.assertTrue(np.allclose(fs2, fsm, rtol=1e-2, atol=1e-2))
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
# phi = dadi.PhiManip.phi_3D_to_4D(phi, 0, 0, xx,xx,xx,xx)
# phi = dadi.Integration.four_pops(phi, xx, T=0.1, nu1=0.5, nu3=10, m13=2, m31=0.5, gamma1=-1, gamma3=1, h1=0.2, h3=0.7)
# phi = dadi.PhiManip.phi_4D_admix_into_3(phi,0.8,0,0,xx,xx,xx,xx)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((1,3))
# self.assertTrue(np.allclose(fs2, fsm, rtol=1e-2, atol=1e-2))
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
# phi = dadi.PhiManip.phi_2D_admix_2_into_1(phi,0.8,xx,xx)
# fs2 = dadi.Spectrum.from_phi(phi, [5,5], (xx,xx))
# phi = dadi.PhiManip.phi_1D(xx)
# phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
# phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
# phi = dadi.PhiManip.phi_3D_to_4D(phi, 0, 0, xx,xx,xx,xx)
# phi = dadi.Integration.four_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.7)
# phi = dadi.PhiManip.phi_4D_admix_into_1(phi,0.8,0,0,xx,xx,xx,xx)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((2,3))
# self.assertTrue(np.allclose(fs2, fsm, rtol=1e-2, atol=1e-2))
suite = unittest.TestLoader().loadTestsFromTestCase(Int5DTestCase)
def test_admix_into(self):
"""
Test phi_5D_admix_into methods
"""
pts = 5
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
ref = dadi.PhiManip.phi_2D_admix_1_into_2(phi,0.8,xx,xx)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
phi = dadi.PhiManip.phi_3D_to_4D(phi, 0,0, xx,xx,xx,xx)
phi = dadi.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
phi = dadi.PhiManip.phi_5D_admix_into_2(phi,0.8,0,0,0,xx,xx,xx,xx,xx)
test = dadi.PhiManip.remove_pop(phi,xx,5)
test = dadi.PhiManip.remove_pop(test,xx,4)
test = dadi.PhiManip.remove_pop(test,xx,3)
self.assertTrue(np.allclose(ref, test))
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
phi = dadi.PhiManip.phi_3D_to_4D(phi, 0,0, xx,xx,xx,xx)
phi = dadi.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
phi = dadi.PhiManip.phi_5D_admix_into_3(phi,0.8,0,0,0,xx,xx,xx,xx,xx)
test = dadi.PhiManip.remove_pop(phi,xx,5)
test = dadi.PhiManip.remove_pop(test,xx,4)
test = dadi.PhiManip.remove_pop(test,xx,2)
self.assertTrue(np.allclose(ref, test))
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
phi = dadi.PhiManip.phi_3D_to_4D(phi, 0,0, xx,xx,xx,xx)
phi = dadi.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
phi = dadi.PhiManip.phi_5D_admix_into_4(phi,0.8,0,0,0,xx,xx,xx,xx,xx)
test = dadi.PhiManip.remove_pop(phi,xx,5)
test = dadi.PhiManip.remove_pop(test,xx,3)
test = dadi.PhiManip.remove_pop(test,xx,2)
self.assertTrue(np.allclose(ref, test))
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=0.1, nu1=0.5, nu2=10, m12=2, m21=0.5, gamma1=-1, gamma2=1, h1=0.2, h2=0.9)
phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
phi = dadi.PhiManip.phi_3D_to_4D(phi, 0,0, xx,xx,xx,xx)
phi = dadi.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
phi = dadi.PhiManip.phi_5D_admix_into_5(phi,0.8,0,0,0,xx,xx,xx,xx,xx)
test = dadi.PhiManip.remove_pop(phi,xx,4)
test = dadi.PhiManip.remove_pop(test,xx,3)
test = dadi.PhiManip.remove_pop(test,xx,2)
self.assertTrue(np.allclose(ref, test))
def test_4D_to_5D(self):
"""
Test splitting with admixture
"""
# Not a comprhensive test, since it only considers a single scenario
pts = 4
xx = dadi.Numerics.default_grid(pts)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=1, nu1=10, nu2=0.1)
# Create pop 3 as a mixture of 1&2
ref = dadi.PhiManip.phi_2D_to_3D(phi, 0.3, xx,xx,xx)
phi = dadi.PhiManip.phi_1D(xx)
phi = dadi.PhiManip.phi_1D_to_2D(xx, phi)
phi = dadi.Integration.two_pops(phi, xx, T=1, nu1=10, nu2=0.1)
phi = dadi.PhiManip.phi_2D_to_3D(phi, 0, xx,xx,xx)
phi = dadi.PhiManip.phi_3D_to_4D(phi, 0, 0, xx,xx,xx,xx)
# Create pop 5 as a mixture of 1&2
phi = dadi.PhiManip.phi_4D_to_5D(phi, 0.3, 0.7, 0, xx,xx,xx,xx,xx)
# Remove pops 4 and 3
test = dadi.PhiManip.remove_pop(phi,xx,4)
test = dadi.PhiManip.remove_pop(test,xx,3)
self.assertTrue(np.allclose(ref, test))
# Don't run these tests by default, because they are very slow
#suite = unittest.TestLoader().loadTestsFromTestCase(Int5DTestCase)
if __name__ == '__main__':
unittest.main()
\ No newline at end of file
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment