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

5D integration tests

parent d8a33c53
......@@ -41,8 +41,7 @@ class Int5DTestCase(unittest.TestCase):
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,
frozen1=True, frozen2=True, frozen3=True, frozen4=True, frozen5=False)
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
......@@ -53,56 +52,102 @@ class Int5DTestCase(unittest.TestCase):
fsm = fs5.marginalize(tomarg)
self.assertTrue(np.allclose(fs1, fsm, rtol=5e-2, atol=5e-2))
#def test_integration_2Dcomp(self):
# """
# Integration tested by comparison to 2D integrations
# """
# pts = 20
# nu1 = lambda t: 0.5 + 5*t
# nu2 = lambda t: 10-50*t
# m12 = lambda t: 2-t
# m21 = lambda t: 0.5+3*t
# gamma1 = lambda t: -2*t
# gamma2 = lambda t: 3*t
# h1 = lambda t: 0.2+t
# h2 = lambda t: 0.9-t
# 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=nu1, nu2=nu2,
# m12=m12, m21=m21, gamma1=gamma1, gamma2=gamma2, h1=h1, h2=h2)
# 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=nu1, nu4=nu2,
# m14=m12, m41=m21, gamma1=gamma1, gamma4=gamma2, h1=h1, h4=h2)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((1,2))
# 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, nu2=nu1, nu4=nu2,
# m24=m12, m42=m21, gamma2=gamma1, gamma4=gamma2, h2=h1, h4=h2)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((0,2))
# 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, nu3=nu1, nu4=nu2,
# m34=m12, m43=m21, gamma3=gamma1, gamma4=gamma2, h3=h1, h4=h2)
# 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))
def test_integration_nomig(self):
nu1 = lambda t: 0.5 + 50*t
nu2 = lambda t: 10-50*t
gamma1 = lambda t: -30*t
gamma2 = lambda t: 30*t
h1 = lambda t: 0.2+5*t
h2 = lambda t: 0.9-5*t
T = 0.1
pin = (T, nu1, nu2, gamma1, gamma2, h1, h2)
@dadi.Numerics.make_extrap_func
def ref_func(params, ns, pts):
T, nu1, nu2, gamma1, gamma2, h1, h2 = params
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=T, nu1=nu1, gamma1=gamma1, h1=h1,
nu2=nu2, gamma2=gamma2, h2=h2)
return dadi.Spectrum.from_phi(phi, [5,5], (xx,xx))
@dadi.Numerics.make_extrap_func
def test_func_all(params, ns, pts):
T, nu1, nu2, gamma1, gamma2, h1, h2 = params
xx = dadi.Numerics.default_grid(pts)
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.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
phi = dadi.Integration.five_pops(phi, xx, T=T, nu1=nu1, gamma1=gamma1, h1=h1,
nu2=nu2, gamma2=gamma2, h2=h2,
nu3=nu2, gamma3=gamma2, h3=h2,
nu4=nu2, gamma4=gamma2, h4=h2,
nu5=nu2, gamma5=gamma2, h5=h2)
fs_all = dadi.Spectrum.from_phi(phi, [5,5,5,5,5], (xx,xx,xx,xx,xx))
return fs_all
ref = ref_func(pin, None, [16,18,20])
fs_all = test_func_all(pin, None, [16,18,20])
fs2 = fs_all.marginalize([2,3,4])
fs3 = fs_all.marginalize([1,3,4])
fs4 = fs_all.marginalize([1,2,4])
fs5 = fs_all.marginalize([1,2,3])
self.assertTrue(np.allclose(fs2, ref, atol=1e-2))
self.assertTrue(np.allclose(fs3, ref, atol=1e-2))
self.assertTrue(np.allclose(fs4, ref, atol=1e-2))
self.assertTrue(np.allclose(fs5, ref, atol=1e-2))
def test_integration_mig(self):
"""
Integration tested by comparison to 2D integrations
"""
m12 = lambda t: 2-19*t
m21 = lambda t: 0.5+30*t
T = 0.1
pin = (T,m12,m21)
@dadi.Numerics.make_extrap_func
def ref_func(params, ns, pts):
T, m12, m21 = params
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=T, m12=m12, m21=m21)
return dadi.Spectrum.from_phi(phi, [5,5], (xx,xx))
@dadi.Numerics.make_extrap_func
def test_func_marg(params, ns, popkept, pts):
T, m12, m21 = params
xx = dadi.Numerics.default_grid(pts)
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.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
kwargs = {'m1{0}'.format(popkept):m12,
'm{0}1'.format(popkept):m21}
phi = dadi.Integration.five_pops(phi, xx, T=T, **kwargs)
fs_all = dadi.Spectrum.from_phi(phi, [5,5,5,5,5], (xx,xx,xx,xx,xx))
tomarg = [1,2,3,4]
tomarg.remove(popkept-1)
fs = fs_all.marginalize(tomarg)
return fs
ref = ref_func(pin, None, [16,18,20])
for ii in range(2,6):
test = test_func_marg(pin, None, 2, [16, 18, 20])
self.assertTrue(np.allclose(ref, test, atol=1e-3))
#def test_admix_props(self):
# """
......@@ -186,128 +231,4 @@ class Int5DTestCase(unittest.TestCase):
suite = unittest.TestLoader().loadTestsFromTestCase(Int5DTestCase)
if __name__ == '__main__':
#cm = Int5DTestCase()
#cm.test_integration_SNM()
#unittest.main()
pts = 100
nu1 = lambda t: 0.5 + 50*t
nu2 = lambda t: 10-50*t
gamma1 = lambda t: -30*t
gamma2 = lambda t: 30*t
h1 = lambda t: 0.2+5*t
h2 = lambda t: 0.9-5*t
T = 0.1
nu2, gamma2, h2 = 1, 0, 0.5
pin = (T, nu1, nu2, gamma1, gamma2, h1, h2)
@dadi.Numerics.make_extrap_func
def ref_func(params, ns, pts):
T, nu1, nu2, gamma1, gamma2, h1, h2 = params
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=T, nu1=nu1, gamma1=gamma1, h1=h1,
nu2=nu2, gamma2=gamma2, h2=h2)
return dadi.Spectrum.from_phi(phi, [5,5], (xx,xx))
@dadi.Numerics.make_extrap_func
def test_func_all(params, ns, pts):
T, nu1, nu2, gamma1, gamma2, h1, h2 = params
xx = dadi.Numerics.default_grid(pts)
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.PhiManip.phi_4D_to_5D(phi, 0,0,0, xx,xx,xx,xx,xx)
phi = dadi.Integration.five_pops(phi, xx, T=T, nu1=nu1, gamma1=gamma1, h1=h1,
nu2=nu2, gamma2=gamma2, h2=h2,
nu3=nu2, gamma3=gamma2, h3=h2,
nu4=nu2, gamma4=gamma2, h4=h2,
nu5=nu2, gamma5=gamma2, h5=h2,
frozen2=False, frozen3=False, frozen4=False, frozen5=False)
fs_all = dadi.Spectrum.from_phi(phi, [5,5,5,5,5], (xx,xx,xx,xx,xx))
return fs_all
@dadi.Numerics.make_extrap_func
def test_func_four(params, ns, pts):
T, nu1, nu2, gamma1, gamma2, h1, h2 = params
xx = dadi.Numerics.default_grid(pts)
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=T, nu1=nu1, gamma1=gamma1, h1=h1,
nu2=nu2, gamma2=gamma2, h2=h2,
nu3=nu2, gamma3=gamma2, h3=h2,
nu4=nu2, gamma4=gamma2, h4=h2)
fs_all = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
return fs_all
dadi.cuda_enabled(False)
ref = ref_func(pin, None, [100,120,140])
fs_all = test_func_all(pin, None, [16,18,20])
fs2 = fs_all.marginalize([2,3,4])
fs3 = fs_all.marginalize([1,3,4])
fs4 = fs_all.marginalize([1,2,4])
fs5 = fs_all.marginalize([1,2,3])
import matplotlib.pyplot as plt
plt.close('all')
for q in [fs2,fs3,fs4,fs5]:
#for q in [fs2,fs3,fs4]:
plt.figure()
dadi.Plotting.plot_2d_comp_Poisson(ref, q)
if dadi.cuda_enabled(True):
fs_all_cuda = test_func_all(pin, None, [16,18,20])
fs2 = fs_all_cuda.marginalize([2,3,4])
fs3 = fs_all_cuda.marginalize([1,3,4])
fs4 = fs_all_cuda.marginalize([1,2,4])
fs5 = fs_all_cuda.marginalize([1,2,3])
for q in [fs2,fs3,fs4,fs5]:
plt.figure()
dadi.Plotting.plot_2d_comp_Poisson(ref, q)
#fs4_all = test_func_four(pin, None, [16,18,20])
#fs2_4 = fs4_all.marginalize([2,3])
#fs3_4 = fs4_all.marginalize([1,3])
#fs4_4 = fs4_all.marginalize([1,2])
#for q in [fs2_4,fs3_4,fs4_4]:
# plt.figure()
# dadi.Plotting.plot_2d_comp_Poisson(ref, q)
#fsm = fs5.marginalize((1,2,3))
#print(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, nu2=nu1, nu4=nu2,
# m24=m12, m42=m21, gamma2=gamma1, gamma4=gamma2, h2=h1, h4=h2)
# fs4 = dadi.Spectrum.from_phi(phi, [5,5,5,5], (xx,xx,xx,xx))
# fsm = fs4.marginalize((0,2))
# 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, nu3=nu1, nu4=nu2,
# m34=m12, m43=m21, gamma3=gamma1, gamma4=gamma2, h3=h1, h4=h2)
# 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))
#
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