Source code for moments.Demographics2D

"""
Two-population demographic models.
"""
import numpy

import moments


[docs]def snm(ns, pop_ids=None): """ ns = [n1, n2] Standard neutral model with a split but no divergence. :param ns: List of population sizes in first and second populations. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1]) fs = moments.Spectrum(sts) fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1]) fs.pop_ids = pop_ids return fs
[docs]def bottlegrowth(params, ns, pop_ids=None): """ params = (nuB, nuF, T) ns = [n1, n2] Instantanous size change followed by exponential growth with no population split. - nuB: Ratio of population size after instantanous change to ancient population size - nuF: Ratio of contempoary to ancient population size - T: Time in the past at which instantaneous change happened and growth began (in units of 2*Na generations) - n1, n2: Sample sizes of resulting Spectrum. :param params: List of parameters, (nuB, nuF, T). :param ns: List of population sizes in first and second populations. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") nuB, nuF, T = params return bottlegrowth_split_mig((nuB, nuF, 0, T, 0), ns, pop_ids=pop_ids)
[docs]def bottlegrowth_split(params, ns, pop_ids=None): """ params = (nuB, nuF, T, Ts) ns = [n1, n2] Instantanous size change followed by exponential growth then split. - nuB: Ratio of population size after instantanous change to ancient population size - nuF: Ratio of contempoary to ancient population size - T: Time in the past at which instantaneous change happened and growth began (in units of 2*Na generations) - Ts: Time in the past at which the two populations split. - n1, n2: Sample sizes of resulting Spectrum. :param ns: List of population sizes in first and second populations. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") nuB, nuF, T, Ts = params return bottlegrowth_split_mig((nuB, nuF, 0.0, T, Ts), ns, pop_ids=pop_ids)
[docs]def bottlegrowth_split_mig(params, ns, pop_ids=None): """ params = (nuB, nuF, m, T, Ts) ns = [n1, n2] Instantanous size change followed by exponential growth then split with migration. - nuB: Ratio of population size after instantanous change to ancient population size - nuF: Ratio of contempoary to ancient population size - m: Migration rate between the two populations (2*Na*m). - T: Time in the past at which instantaneous change happened and growth began (in units of 2*Na generations) - Ts: Time in the past at which the two populations split. - n1, n2: Sample sizes of resulting Spectrum. :param ns: List of population sizes in first and second populations. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") nuB, nuF, m, T, Ts = params nu_func = lambda t: [nuB * numpy.exp(numpy.log(nuF / nuB) * t / T)] sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1]) fs = moments.Spectrum(sts) fs.integrate(nu_func, T - Ts, dt_fac=0.01) # we split the population fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1]) nu0 = nu_func(T - Ts)[0] nu_func = lambda t: 2 * [nu0 * numpy.exp(numpy.log(nuF / nu0) * t / Ts)] fs.integrate(nu_func, Ts, m=numpy.array([[0, m], [m, 0]])) fs.pop_ids = pop_ids return fs
[docs]def split_mig(params, ns, pop_ids=None): """ Split into two populations of specifed size, with migration. params = (nu1, nu2, T, m) ns = [n1, n2] :param params: Tuple of length 4. - nu1: Size of population 1 after split. - nu2: Size of population 2 after split. - T: Time in the past of split (in units of 2*Na generations) - m: Migration rate between populations (2*Na*m) :param ns: List of length two specifying sample sizes n1 and n2. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") nu1, nu2, T, m = params sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1]) fs = moments.Spectrum(sts) fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1]) fs.integrate([nu1, nu2], T, m=numpy.array([[0, m], [m, 0]])) fs.pop_ids = pop_ids return fs
[docs]def IM(params, ns, pop_ids=None): """ Isolation-with-migration model with exponential pop growth. params = (s, nu1, nu2, T, m12, m21) ns = [n1, n2] :param params: Tuple of length 6. - s: Size of pop 1 after split. (Pop 2 has size 1-s.) - nu1: Final size of pop 1. - nu2: Final size of pop 2. - T: Time in the past of split (in units of 2*Na generations) - m12: Migration from pop 2 to pop 1 (2 * Na * m12) - m21: Migration from pop 1 to pop 2 :param ns: List of population sizes in first and second populations. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") s, nu1, nu2, T, m12, m21 = params sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1]) fs = moments.Spectrum(sts) fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1]) nu1_func = lambda t: s * (nu1 / s) ** (t / T) nu2_func = lambda t: (1 - s) * (nu2 / (1 - s)) ** (t / T) nu_func = lambda t: [nu1_func(t), nu2_func(t)] fs.integrate(nu_func, T, dt_fac=0.01, m=numpy.array([[0, m12], [m21, 0]])) fs.pop_ids = pop_ids return fs
[docs]def IM_pre(params, ns, pop_ids=None): """ params = (nuPre, TPre, s, nu1, nu2, T, m12, m21) ns = [n1, n2] Isolation-with-migration model with exponential pop growth and a size change prior to split. - nuPre: Size after first size change - TPre: Time before split of first size change. - s: Fraction of nuPre that goes to pop1. (Pop 2 has size nuPre*(1-s).) - nu1: Final size of pop 1. - nu2: Final size of pop 2. - T: Time in the past of split (in units of 2*Na generations) - m12: Migration from pop 2 to pop 1 (2*Na*m12) - m21: Migration from pop 1 to pop 2 - n1, n2: Sample sizes of resulting Spectrum. :param ns: List of population sizes in first and second populations. :param pop_ids: List of population IDs. """ if pop_ids is not None and len(pop_ids) != 2: raise ValueError("pop_ids must be a list of two population IDs") nuPre, TPre, s, nu1, nu2, T, m12, m21 = params sts = moments.LinearSystem_1D.steady_state_1D(ns[0] + ns[1]) fs = moments.Spectrum(sts) fs.integrate([nuPre], TPre, dt_fac=0.01) fs = moments.Manips.split_1D_to_2D(fs, ns[0], ns[1]) nu1_0 = nuPre * s nu2_0 = nuPre * (1 - s) nu1_func = lambda t: nu1_0 * (nu1 / nu1_0) ** (t / T) nu2_func = lambda t: nu2_0 * (nu2 / nu2_0) ** (t / T) nu_func = lambda t: [nu1_func(t), nu2_func(t)] fs.integrate(nu_func, T, dt_fac=0.01, m=numpy.array([[0, m12], [m21, 0]])) fs.pop_ids = pop_ids return fs