Specifying a model in details

It is also possible to use the genetic algorithm from GADMA to a proposed model that is defined as a Python function using dadi or moments. It is the way that dadi and moments work with demographic models inference. To understand how to specify a model like that one can read manuals to these packages.

For example, consider a simple bottleneck model for one population: at time TF + TB in the past, an equilibrium population goes through a bottleneck of depth nuB, recovering to relative size nuF:

def bottleneck(params, ns, pts):
    nuB, nuF, TB, TF = params
    xx = dadi.Numerics.default_grid(pts)

    phi = dadi.PhiManip.phi_1D(xx)
    phi = dadi.Integration.one_pop(phi, xx, TB, nuB)
    phi = Integration.one_pop(phi, xx, TF, nuF)

    fs = dadi.Spectrum.from_phi(phi, ns, (xx,))
    return fs

To run optimization from GADMA one needs to run optimization function just like in dadi and moments:

# Import GADMA's optimization function:
import gadma

# Specify input data and its parameters:
data = dadi.Spectrum.from_file("fs_filename.fs")
ns = data.sample_sizes # size of AFS
pts = [40,50,60] # grid size for dadi

# Wrap our bottleneck function:
func_ex = dadi.Numerics.make_extrap_log_func(bottleneck)

# Specify upper and lower bounds for parameters:
upper_bound = [100, 100, 3, 3]
lower_bound = [1e-2, 1e-2, 0, 0]

# Run optimizations:
# Beginning GADMA optimization
popt = gadma.Inference.optimize_ga(data, func_ex, engine='dadi', args=pts_l,
                                   p_ids = ['nuB', 'nuF', 'TB', 'TF'],
                                   lower_bound=lower_bound,
                                   upper_bound=upper_bound)
# Beginning local optimization from dadi
popt = dadi.Inference.optimize_log(popt, data, func_ex, pts_l,
                                   lower_bound=lower_bound,
                                   upper_bound=upper_bound,
                                   verbose=len(p0), maxiter=3)

print('Found parameters: {0}'.format(popt))

Note

As GADMA optimization is a global search algorithm, no initial parameters p0 are set in gadma.Inference.optimize_ga function. However, it is possible to specify X_init with p0 as one of known starting points:

# Initial parameters can be set too:
p0 = [0.01, 1.5, 0.2, 0.2]

# Beginning GADMA optimization
popt = gadma.Inference.optimize_ga(data, func_ex, engine='dadi', args=pts_l, X_init=[p0],
                                   p_ids = ['nuB', 'nuF', 'TB', 'TF'],
                                   lower_bound=lower_bound,
                                   upper_bound=upper_bound)

Warning

Function gadma.Inference.optimize_ga changed in GADMA version 2. For full documentation see the API (Inference module).

If one wants to find other parameters for gadma.Inference.optimize_ga function:

>>> import gadma
>>> help(gadma.Inference.optimize_ga)