<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="lev2 toc-item"><a data-toc-modified-id="Import-dadi-on-all-engines-01" href="#Import-dadi-on-all-engines"><span class="toc-item-num">0.1&nbsp;&nbsp;</span>Import dadi on all engines</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Load-1D-SFS-02" href="#Load-1D-SFS"><span class="toc-item-num">0.2&nbsp;&nbsp;</span>Load 1D SFS</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Exponential-growth-model-03" href="#Exponential-growth-model"><span class="toc-item-num">0.3&nbsp;&nbsp;</span>Exponential growth model</a></div><div class="lev3 toc-item"><a data-toc-modified-id="create-combinations-of-starting-values-for-optimisation-031" href="#create-combinations-of-starting-values-for-optimisation"><span class="toc-item-num">0.3.1&nbsp;&nbsp;</span>create combinations of starting values for optimisation</a></div><div class="lev3 toc-item"><a data-toc-modified-id="ery-032" href="#ery"><span class="toc-item-num">0.3.2&nbsp;&nbsp;</span><em>ery</em></a></div><div class="lev4 toc-item"><a data-toc-modified-id="run-optimisation-on-1D-SFS-of-0321" href="#run-optimisation-on-1D-SFS-of"><span class="toc-item-num">0.3.2.1&nbsp;&nbsp;</span>run optimisation on 1D SFS of</a></div><div class="lev4 toc-item"><a data-toc-modified-id="calculate-log-likelihood-of-each-set-of-optimised-parameters-0322" href="#calculate-log-likelihood-of-each-set-of-optimised-parameters"><span class="toc-item-num">0.3.2.2&nbsp;&nbsp;</span>calculate log likelihood of each set of optimised parameters</a></div><div class="lev4 toc-item"><a data-toc-modified-id="residual-plots-0323" href="#residual-plots"><span class="toc-item-num">0.3.2.3&nbsp;&nbsp;</span>residual plots</a></div><div class="lev4 toc-item"><a data-toc-modified-id="different-optimisation-algorithms-0324" href="#different-optimisation-algorithms"><span class="toc-item-num">0.3.2.4&nbsp;&nbsp;</span>different optimisation algorithms</a></div><div class="lev3 toc-item"><a data-toc-modified-id="par-033" href="#par"><span class="toc-item-num">0.3.3&nbsp;&nbsp;</span><em>par</em></a></div><div class="lev4 toc-item"><a data-toc-modified-id="run-optimisation-on-1D-SFS-of-par-0331" href="#run-optimisation-on-1D-SFS-of-par"><span class="toc-item-num">0.3.3.1&nbsp;&nbsp;</span>run optimisation on 1D SFS of par</a></div><div class="lev1 toc-item"><a data-toc-modified-id="Recommended-parameter-boundaries-1" href="#Recommended-parameter-boundaries"><span class="toc-item-num">1&nbsp;&nbsp;</span>Recommended parameter boundaries</a></div><div class="lev1 toc-item"><a data-toc-modified-id="Restart-analysis-2" href="#Restart-analysis"><span class="toc-item-num">2&nbsp;&nbsp;</span>Restart analysis</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Ery-21" href="#Ery"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>Ery</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Par-22" href="#Par"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>Par</a></div><div class="lev3 toc-item"><a data-toc-modified-id="restart-221" href="#restart"><span class="toc-item-num">2.2.1&nbsp;&nbsp;</span>restart</a></div>
## Import dadi on all engines
The magic commands become automaticaly available when you create a client
Step1:  Dadi should now be imported on all remote engines as well as locally.
Step2:  The remote namespace can be checked with
Step3:  As can be seen, dadi is imported on all engines.
Load 1D SFS
Like the dadi module, I would like to load the 1D folded SFS's of par and ery into the namespaces of the local IPython session as well as on all engines.
Step4:  Exponential growth model
Step5:  This model assumes that exponential growth (or decline) started some time $T$ in the past and the current effective population size is a multiple $\nu$ of the ancient populations size, i. e. before exponential growth began. So this model just estimates two parameters. If $\nu$ is estimated to be 1, this indicates that the population size hasn't changed (although see Myers2008). If it is below one, this indicates exponential population decline (how realistic is that?). If it is above 1, this indicates exponential population growth.
Step6:  create combinations of starting values for optimisation
Step7:  These starting values cover a very wide range of possible parameter combinations. This should be better than merely randomly perturbing a single parameter combination a certain number of times (and running the optimisation from those perturbed starting values).
ery
run optimisation on 1D SFS of
Step9:  dadi also offers other otimisation algorithms
Step10:  It took 7 minutes to run these 100 optimisations.
Step11:  The tuples returned from my parallel function have been flattened into one array. That's not nice.
Step12:  Clearly, different starting parameter values lead to different optimised parameter values. There is a great variance in the returned optimised parameters. Maybe the maximum number of iterations (maxiter) is not enough?
I don't understand the the maxiter option to the optimisation function. Could it be that 100 iterations are not enough to allow some optimisations with starting values far away from the optimal values to find the optimum?
calculate log likelihood of each set of optimised parameters
Step13:  I would like to create a data frame with pandas with these arrays and then sort the table by ll.
Step14:  Given the range of starting values that result in almost exactly the same optimal parameter values, I am fairly convinced that the optimal parameters for this model are close to $\nu = 0.136$ and $T = 0.0075$. Given the model, the optimal parameters indicate that the ancestral populations to the erythropus sample started to decline exponentially (how realistic is that?) $0.0075 \times 2N_a$ generations ago to a contemporary population size of about 14% of the ancient effective population size ($N_a$).
Step16:  The parameter combination $\nu = 0.001049$ and $T = 0.000023$ has only slightly smaller log likelihood. This parameter combination corresponds to a much more recent and much more drastic population reduction. I think a likelihood ratio test would be required to see whether the data can distinguish between these two scenarios.
Maybe I should do a longer run, starting from the best-fit parameter sets found here and setting maxiter to something ridiculously high, so that it is guaranteed that the ML values will be found within gtol deviation.
Step17:  This toke only 6 seconds to compute on 11 engines.
Step18:  There are still three outlier results. Maybe I need to use a different optimsiation function.
residual plots
I would like to comparisons between model predicted and observed SFS for a couple of optimisation results.
Step19:  Compared to the standard neutral model prediction, it seems that only the first frequency class has a markedly reduced residual with this model (about 20 standard deviations).
Step20:  I cannot make out a clear difference between these this and the previous parameter set.
Step22:  Again, even this very different parameter set seems to produce a very similar SFS to the ones produced by the two more likely parameter sets.
different optimisation algorithms
Step23:  Now run dadi optimisation with Nelder-Mead algorithm.
Step25:  Half of the optimisations result in a combination of parameters that do not have the highest likelihood. So it does seem as though Nelder-Mead is more robust than BFGS (but note the random perturbation of starting values).
par
run optimisation on 1D SFS of par
Step26:  This runs definitely too long! run_dadi_on_par spends most of its time on the last three tasks. A maxiter of 100, therefore, does not seem to be limiting the search. Set to lower maxiter next time. Maybe 50.
Recommended parameter boundaries
Ryan Gutenkunst writes on the dadi forum
Step28:  The parameter combination with the highest likelihood is not found several times from different starting values. This does not give me much confidence in the maximum likelihood parameter set.
Step30:  None of these parameter combinations is the most likely determined above! I can therefore not be very confident that I have found the most likely parameter combination.
Maybe another optimisation algorithm is better at finding the global maximum.
Step31:  Again, none of these parameter combinations is close to the most likely ones determined above.
Let's see how the best-fitting model and data compare to each other.
Step32:  These logarithmic plots are not very helpfull.
Step33:  So, if I could believe the optimal parameter combination found here to be the global maximum, it would imply that the ancestral populations of the parallelus sample have undergone an exponential decline to about 1/1000th of the ancient population size, $N_{ref}$ (i. e. very strong bottleneck) and this happend in the very recent past, i. e. only $2.3 \times 10^{-5} \times 2N_{ref}$ generations ago.
Step34:  Most of the time this just runs on a single core! I don't know why. Note the difference between the wall time and the total time.
Running the simulation locally on a single core
Step35:  It really shouldn't take 2 minutes to run this 20 times with different starting values.
Step36:  Note, that p1 and popt are not in the namespace of the remote engines. They apparently have been cleared automatically.
Step38:  Jobs submitted to the load balanced view run asynchonously in non-blocking mode.
Step39:  I think it should be a good idea to use load balanced views of the engines, since I think some optimisations take longer than others.
Step40:  Look at max_iter from dadi.Inference.optimize_log.
Restart analysis
Restart the kernel from the menu to clear the namespace!
Step42:  Use unfolded ANGSD spectra and fold with dadi
Step43:  Ery
The following sets global variables on the remote engines that are called by run_dadi
Step44:  This lets the optimisations run with the faster BFGS algorithm and stop after a maximum of 10 iterations. This is intended for a broad sweep across the parameter space. Also note that this sets the optimisations to fit the ery spectrum.
Step45:  I needed to kill a single very long lasting optimisation job from the terminal (because it is not possible to do that with ipyparallel). Unfortunately that also means that the results from all 224 completed (or terminated) optimisation runs are lost as well. This is the major shortcoming of ipyparallel at the moment.
Couldn't I make ipyparallel pickle the results of the optimisations to file?
Step47:  Now, define a new run_dadi function that pickles output to file
Step48:  I killed the last optimisation.
Load pickled output
Step50:  Get flag count from optimisations
Step52:  There were 39 successfull optimisations.
Create table of results from successful optimisations
Step53:  These results are completely inconclusive. I think it would make sense to limit the search space for $T$ to something well above 0 and start a few optimisations from parameter values close to a neutral model.
Step54:  None of these optimisations was able to converge on optimal parameter values.
Step55:  The perturbation has not been very effective.
Step56:  Why are there only 19 output files? There should be 20! Two output filenames happened to be identical.
Step57:  OK. Three successful optimisations.
Step58:  OK, that looks promising. Two quite different parameter combinations have lead to the same optimal parameter combination.
Step59:  This looks like convergence.
Step60:  This is still extremely fast!
Step61:  Par
Since starting optimisations from a broad range of initial values did not lead to useful results with the ery spectrum, I am going to use the same strategy that led to a successful optimisation with ery
Step62:  I needed to kill the last two processes in the terminal because they ran more than 15 minutes. I cannot use the AsyncResult object because of that and unstead need to load the pickled output from file.
Step63:  Only 1 successful optimisation.
Step64:  That is hitting the upper limit for the time parameter.
Maybe a sweep through the parameter space could help to find an optimum.
Step65:  Unfortunately, this is too slow. I need to abort this by stopping the cluster of engines.
Step67:  After restarting the cluster, it is important to be EXTREMELY careful to load all required variables and modules on all engines again.
restart
Step69:  I am trying to reduce the initial parameter sweep.
<ASSISTANT_TASK:>
Python Code:
from ipyparallel import Client
cl = Client()
cl.ids
cl[:].targets
%%px --noblock
# run a whole cell in non-blocking mode, by default on all engines
# note: the magic has to be at the top of the cell
import time
time.sleep(1)
time.time()
%pxresult
# get the result from the AsyncResult object
%%px --local
# run whole cell on all engines a well as in the local IPython session
import sys
sys.path.insert(0, '/home/claudius/Downloads/dadi')
sys.path[0]
%%px --local
import dadi
Explanation: Table of Contents
<p><div class="lev2 toc-item"><a data-toc-modified-id="Import-dadi-on-all-engines-01" href="#Import-dadi-on-all-engines"><span class="toc-item-num">0.1&nbsp;&nbsp;</span>Import dadi on all engines</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Load-1D-SFS-02" href="#Load-1D-SFS"><span class="toc-item-num">0.2&nbsp;&nbsp;</span>Load 1D SFS</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Exponential-growth-model-03" href="#Exponential-growth-model"><span class="toc-item-num">0.3&nbsp;&nbsp;</span>Exponential growth model</a></div><div class="lev3 toc-item"><a data-toc-modified-id="create-combinations-of-starting-values-for-optimisation-031" href="#create-combinations-of-starting-values-for-optimisation"><span class="toc-item-num">0.3.1&nbsp;&nbsp;</span>create combinations of starting values for optimisation</a></div><div class="lev3 toc-item"><a data-toc-modified-id="ery-032" href="#ery"><span class="toc-item-num">0.3.2&nbsp;&nbsp;</span><em>ery</em></a></div><div class="lev4 toc-item"><a data-toc-modified-id="run-optimisation-on-1D-SFS-of-0321" href="#run-optimisation-on-1D-SFS-of"><span class="toc-item-num">0.3.2.1&nbsp;&nbsp;</span>run optimisation on 1D SFS of</a></div><div class="lev4 toc-item"><a data-toc-modified-id="calculate-log-likelihood-of-each-set-of-optimised-parameters-0322" href="#calculate-log-likelihood-of-each-set-of-optimised-parameters"><span class="toc-item-num">0.3.2.2&nbsp;&nbsp;</span>calculate log likelihood of each set of optimised parameters</a></div><div class="lev4 toc-item"><a data-toc-modified-id="residual-plots-0323" href="#residual-plots"><span class="toc-item-num">0.3.2.3&nbsp;&nbsp;</span>residual plots</a></div><div class="lev4 toc-item"><a data-toc-modified-id="different-optimisation-algorithms-0324" href="#different-optimisation-algorithms"><span class="toc-item-num">0.3.2.4&nbsp;&nbsp;</span>different optimisation algorithms</a></div><div class="lev3 toc-item"><a data-toc-modified-id="par-033" href="#par"><span class="toc-item-num">0.3.3&nbsp;&nbsp;</span><em>par</em></a></div><div class="lev4 toc-item"><a data-toc-modified-id="run-optimisation-on-1D-SFS-of-par-0331" href="#run-optimisation-on-1D-SFS-of-par"><span class="toc-item-num">0.3.3.1&nbsp;&nbsp;</span>run optimisation on 1D SFS of par</a></div><div class="lev1 toc-item"><a data-toc-modified-id="Recommended-parameter-boundaries-1" href="#Recommended-parameter-boundaries"><span class="toc-item-num">1&nbsp;&nbsp;</span>Recommended parameter boundaries</a></div><div class="lev1 toc-item"><a data-toc-modified-id="Restart-analysis-2" href="#Restart-analysis"><span class="toc-item-num">2&nbsp;&nbsp;</span>Restart analysis</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Ery-21" href="#Ery"><span class="toc-item-num">2.1&nbsp;&nbsp;</span>Ery</a></div><div class="lev2 toc-item"><a data-toc-modified-id="Par-22" href="#Par"><span class="toc-item-num">2.2&nbsp;&nbsp;</span>Par</a></div><div class="lev3 toc-item"><a data-toc-modified-id="restart-221" href="#restart"><span class="toc-item-num">2.2.1&nbsp;&nbsp;</span>restart</a></div>
## Import dadi on all engines
The magic commands become automaticaly available when you create a client:
End of explanation
%whos
Explanation: Dadi should now be imported on all remote engines as well as locally.
End of explanation
%px %who
Explanation: The remote namespace can be checked with:
End of explanation
%ll
% ll dadiExercises/
! cat dadiExercises/ERY.FOLDED.sfs.dadi_format
! cat dadiExercises/PAR.FOLDED.sfs.dadi_format
%%px --local
# import 1D spectrum of ery on all engines:
fs_ery = dadi.Spectrum.from_file('dadiExercises/ERY.FOLDED.sfs.dadi_format')
# import 1D spectrum of ery on all engines:
fs_par = dadi.Spectrum.from_file('dadiExercises/PAR.FOLDED.sfs.dadi_format')
%px %who
%who
%%px --local
ns = fs_ery.sample_sizes # both populations have the same sample size
fs_ery.pop_ids = ['ery']
fs_par.pop_ids = ['par']
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [50, 60, 70]
Explanation: As can be seen, dadi is imported on all engines.
Load 1D SFS
Like the dadi module, I would like to load the 1D folded SFS's of par and ery into the namespaces of the local IPython session as well as on all engines.
End of explanation
%psource dadi.Demographics1D.growth
Explanation: Exponential growth model
End of explanation
%%px --local
# create link to function that specifies a simple growth or decline model
func = dadi.Demographics1D.growth
# create extrapolating version of the function
func_ex = dadi.Numerics.make_extrap_log_func(func)
# set lower and upper bounds to nu and T
upper_bound = [100, 5]
lower_bound = [1e-3, 0]
%px %who
Explanation: This model assumes that exponential growth (or decline) started some time $T$ in the past and the current effective population size is a multiple $\nu$ of the ancient populations size, i. e. before exponential growth began. So this model just estimates two parameters. If $\nu$ is estimated to be 1, this indicates that the population size hasn't changed (although see Myers2008). If it is below one, this indicates exponential population decline (how realistic is that?). If it is above 1, this indicates exponential population growth.
End of explanation
from itertools import product
import numpy as np
from __future__ import print_function
# create range of starting values for nu and T, evenly distributed in log space
p0_nu = np.logspace(-3, np.log10(50), base=10.0, num=10)
p0_T = np.logspace(-3, np.log10(2), base=10, num=10)
# create all possible pairwise combinations of nu and T
i = 0
num_comb = len(p0_nu)*len(p0_T)
for p0 in product(p0_nu, p0_T):
    i+=1
    if i>20 and i<=(num_comb-20): 
        continue # continue with the next iteration of the loop
    print( i , p0 )
print("The total number of parameter combinations is ", i)
Explanation: create combinations of starting values for optimisation
End of explanation
# create load balanced view of engines
lbview = cl.load_balanced_view()
%pdoc dadi.Inference.optimize_log
Explanation: These starting values cover a very wide range of possible parameter combinations. This should be better than merely randomly perturbing a single parameter combination a certain number of times (and running the optimisation from those perturbed starting values).
ery
run optimisation on 1D SFS of
End of explanation
# create parallel function for ery spectrum
@lbview.parallel(block=False)
def run_dadi_on_ery(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    # run optimisation of paramters
    popt = dadi.Inference.optimize_log(p0=p_init, data=fs_ery, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=100, full_output=False)
    return (p_init, popt) # returns a tuple with starting values and optimal nu and T
# run dadi in parallel on all combinations of starting values
ar = run_dadi_on_ery( product(p0_nu, p0_T) )
# wait for results to come in and print progress report
#ar.wait_interactive()
ar.progress
# run after optimisations have finished
ar.wall_time/60
Explanation: dadi also offers other otimisation algorithms: optimize_log_fmin and optimize_grid.
End of explanation
# save optimsation results to file
# can be loaded back into another IPython session with: ar = pickle.load(filehandle)
# unfortunately, I cannot save the AsyncResult object with pickle (dill can save it, but emits a 
# "maximum recursion depth exceeded" error when I try to load the object from file again).
# I need to save the result as a plain list
import dill
with open("exp_growth_optim_res_ERY.pickle", "w") as fhandle:
    dill.dump(list(ar.get()), fhandle)
# this can now be loaded back from file into memory:
with open("exp_growth_optim_res_ERY.pickle", "r") as fhandle:
    ra = dill.load(fhandle)
?pickle
ar.get()[:10]
Explanation: It took 7 minutes to run these 100 optimisations.
End of explanation
?print # print_function from __future__
# print result of optimisations
for i, p in enumerate(ar):
    if (i+1) % 2:
        #print(i, p)
        print("starting values: nu = %.4f, T = %.6f" % (p[0][0], p[0][1]), end=" ==> ")
    else:
        #print(i, p)
        print("optim. values: nu = %.4f, T = %.6f" % (p[0], p[1]))    
Explanation: The tuples returned from my parallel function have been flattened into one array. That's not nice.
End of explanation
# get the optimized parameter values
popt = [p for i, p in enumerate(ar) if not (i+1) % 2]
print( popt[:10] )
len(popt)
# get the starting values
pstart = [p for i, p in enumerate(ar) if (i+1) % 2]
pstart[:10]
# get log likelihood values for parameter combinations
@lbview.parallel(block=True)
def get_ll(p):
    # calculate the best-fit model SFS given paramter combination
    expected_sfs = func_ex(p, ns, pts_l)
    # calculate the log likelihood of the expected SFS given the observed SFS
    return dadi.Inference.ll_multinom(expected_sfs, fs_ery)
# run get_ll on all sets of parameters in parallel
ll = get_ll.map(popt)
sorted(ll, reverse=True)[:20]
# get arrays of individual parameters from each optimisation run
nu = [p_set[0] for p_set in popt]
T = [p_set[1] for p_set in popt]
nu_start = [p_set[0][0] for p_set in pstart]
T_start = [p_set[0][1] for p_set in pstart]
print(len(ll), len(nu), len(T), len(nu_start), len(T_start))
Explanation: Clearly, different starting parameter values lead to different optimised parameter values. There is a great variance in the returned optimised parameters. Maybe the maximum number of iterations (maxiter) is not enough?
I don't understand the the maxiter option to the optimisation function. Could it be that 100 iterations are not enough to allow some optimisations with starting values far away from the optimal values to find the optimum?
calculate log likelihood of each set of optimised parameters
End of explanation
import pandas as pd
# create pandas data frame
df = pd.DataFrame({
        'a_nu_start' : nu_start,
        'b_T_start' : T_start,
        'c_nu' : nu,
        'd_T' : T,
        'e_log ll' : ll
    })
df.head()
df.dtypes
# get statistics for each column
df.describe()
df.sort_values(by='e_log ll', ascending=False).head(50)
Explanation: I would like to create a data frame with pandas with these arrays and then sort the table by ll.
End of explanation
%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams['font.size'] = 20.0
plt.figure(figsize=[12.0, 10.0])
plt.scatter(np.log(nu), np.log(T), c=ll, s=75, alpha=0.5)
plt.xlabel(r'log $\nu$')
plt.ylabel("log T")
plt.title("log likelihoods for parameter combinations")
plt.grid()
plt.colorbar()
Explanation: Given the range of starting values that result in almost exactly the same optimal parameter values, I am fairly convinced that the optimal parameters for this model are close to $\nu = 0.136$ and $T = 0.0075$. Given the model, the optimal parameters indicate that the ancestral populations to the erythropus sample started to decline exponentially (how realistic is that?) $0.0075 \times 2N_a$ generations ago to a contemporary population size of about 14% of the ancient effective population size ($N_a$).
End of explanation
p0 = df.sort_values(by='e_log ll', ascending=False).iloc[0, 2:4]
p_init = np.array(p0)
p_init
# create direct view of engines
dview = cl[:]
# define function to run dadi with perturbed starting values
# note the high maxiter
def run_dadi_on_ery(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    # perturb starting values by up to a factor of 2 times 'fold'
    p0 = dadi.Misc.perturb_params(p_init, fold=2, upper_bound=upper_bound, lower_bound=lower_bound)
    # run optimisation of paramters
    popt = dadi.Inference.optimize_log(p0=p0, data=fs_ery, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=1e9, full_output=False)
    return popt # returns a tuple with starting values and optimal nu and T
%pdoc run_dadi_on_ery
# run optimisation with perturbed initial values on all 11 engines
fine_opt_res = dview.apply_async( run_dadi_on_ery, p_init )
# get running time
fine_opt_res.wall_time
Explanation: The parameter combination $\nu = 0.001049$ and $T = 0.000023$ has only slightly smaller log likelihood. This parameter combination corresponds to a much more recent and much more drastic population reduction. I think a likelihood ratio test would be required to see whether the data can distinguish between these two scenarios.
Maybe I should do a longer run, starting from the best-fit parameter sets found here and setting maxiter to something ridiculously high, so that it is guaranteed that the ML values will be found within gtol deviation.
End of explanation
# get optimised parameters
fine_opt_res.get()
# run get_ll on all sets of parameters in parallel
ll = get_ll.map(fine_opt_res.get())
ll
nu_opt = [p[0] for p in fine_opt_res.get()]
T_opt = [p[1] for p in fine_opt_res.get()]
df = pd.DataFrame({
        'a_nu' : nu_opt,
        'b_T' : T_opt,
        'c_ll' : ll
    })
df
Explanation: This toke only 6 seconds to compute on 11 engines.
End of explanation
p_init
# calculate the best-fit model SFS given the paramter combination with the highest likelihood
expected_sfs = func_ex(p_init, ns, pts_l)
expected_sfs.fold()
plt.rcParams['figure.figsize'] = [12.0, 10.0]
print(p_init)
dadi.Plotting.plot_1d_comp_multinom(expected_sfs.fold()[:19], fs_ery[:19], residual='linear')
Explanation: There are still three outlier results. Maybe I need to use a different optimsiation function.
residual plots
I would like to comparisons between model predicted and observed SFS for a couple of optimisation results.
End of explanation
p = np.array(df.iloc[0, 0:2])
print(p)
expected_sfs = func_ex(p_init, ns, pts_l)
dadi.Plotting.plot_1d_comp_multinom(expected_sfs.fold()[:19], fs_ery[:19], residual='linear')
Explanation: Compared to the standard neutral model prediction, it seems that only the first frequency class has a markedly reduced residual with this model (about 20 standard deviations).
End of explanation
df
p = np.array(df.iloc[6, 0:2])
print(p)
expected_sfs = func_ex(p_init, ns, pts_l)
dadi.Plotting.plot_1d_comp_multinom(expected_sfs.fold()[:19], fs_ery[:19], residual='linear')
Explanation: I cannot make out a clear difference between these this and the previous parameter set.
End of explanation
# define function to run dadi with perturbed starting values
def run_dadi_on_ery(p_init, dadi_opt_func, maxiter=100): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to perturb and then run optimisation from
    dadi_opt_func:
        + optimize_log (BFGS based)
        + optimize_log_fmin (Nelder-Mead)
    maxiter: maximum number of iterations
    
    # perturb starting values by up to a factor of 2 times 'fold'
    p0 = dadi.Misc.perturb_params(p_init, fold=2, upper_bound=upper_bound, lower_bound=lower_bound)
    # run optimisation of paramters
    popt = dadi_opt_func(p0=p0, data=fs_ery, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=maxiter, full_output=False)
    return popt # returns a tuple with starting values and optimal nu and T
%pdoc run_dadi_on_ery
Explanation: Again, even this very different parameter set seems to produce a very similar SFS to the ones produced by the two more likely parameter sets.
different optimisation algorithms
End of explanation
# run optimisation with perturbed initial values on all 11 engines
fine_opt_res = dview.apply_async( run_dadi_on_ery, p_init, dadi.Inference.optimize_log_fmin, maxiter=1e9 )
fine_opt_res.wall_time
ll = get_ll.map( fine_opt_res.get() )
nu = [p[0] for p in fine_opt_res.get()]
T = [p[1] for p in fine_opt_res.get()]
df = pd.DataFrame({
        'a_nu' : nu,
        'b_T' : T,
        'c_ll' : ll
    })
df
Explanation: Now run dadi optimisation with Nelder-Mead algorithm.
End of explanation
# create parallel function for par spectrum
@lbview.parallel(block=False)
def run_dadi_on_par(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    # run optimisation of paramters
    popt = dadi.Inference.optimize_log(p0=p_init, data=fs_par, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=100, full_output=False)
    return (p_init, popt) # returns a tuple with starting values and optimal nu and T
# run dadi in parallel on all combinations of starting values
#ar = run_dadi_on_par( product(p0_nu, p0_T) )
# wait for results to come in and print progress report
#ar.wait_interactive()
ar.wall_time/60
Explanation: Half of the optimisations result in a combination of parameters that do not have the highest likelihood. So it does seem as though Nelder-Mead is more robust than BFGS (but note the random perturbation of starting values).
par
run optimisation on 1D SFS of par
End of explanation
%psource dadi.Inference.optimize_log
# save optimsation results to file
# can be loaded back into another IPython session with: ar = pickle.load(filehandle)
# unfortunately, I cannot save the AsyncResult object with pickle (dill can save it, but emits a 
# "maximum recursion depth exceeded" error when I try to load the object from file again).
# I need to save the result as a plain list
import dill
with open("exp_growth_optim_res_PAR.pickle", "w") as fhandle:
    dill.dump(list(ar.get()), fhandle)
# get the optimized parameter values
popt = [p for i, p in enumerate(ar) if not (i+1) % 2]
# get the starting values
pstart = [p for i, p in enumerate(ar) if (i+1) % 2]
# get arrays of individual parameters from each optimisation run
nu = [p_set[0] for p_set in popt]
T = [p_set[1] for p_set in popt]
nu_start = [p_set[0][0] for p_set in pstart]
T_start = [p_set[0][1] for p_set in pstart]
# run get_ll on all sets of parameters in parallel
ll = get_ll.map(popt)
# create pandas data frame
df = pd.DataFrame({
        'a_nu_start' : nu_start,
        'b_T_start' : T_start,
        'c_nu' : nu,
        'd_T' : T,
        'e_log ll' : ll
    })
df.head()
df.sort_values(by='e_log ll', ascending=False).head(30)
Explanation: This runs definitely too long! run_dadi_on_par spends most of its time on the last three tasks. A maxiter of 100, therefore, does not seem to be limiting the search. Set to lower maxiter next time. Maybe 50.
Recommended parameter boundaries
Ryan Gutenkunst writes on the dadi forum:
The parameter bounds are primarily about avoiding really slow fits. Calculation with high migration rates and times or small population sizes is much slower, and the optimizer can explore extreme values before settling down. So typically we set migrations rates to be bounded [0, 20ish], times to be [0, 5ish], population sizes to be [1e-3, 1e6ish].
Maybe I should set the time parameter to something above 0.
End of explanation
plt.rcParams['font.size'] = 20.0
plt.figure(figsize=[12.0, 10.0])
plt.scatter(np.log(nu), np.log(T), c=ll, s=75, alpha=0.5)
plt.xlabel(r'log $\nu$')
plt.ylabel("log T")
plt.title("log likelihoods for parameter combinations")
plt.grid()
plt.colorbar()
# define function to run dadi with perturbed starting values
# note the high maxiter
def run_dadi_on_par(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    # perturb starting values by up to a factor of 2 times 'fold'
    p0 = dadi.Misc.perturb_params(p_init, fold=2, upper_bound=upper_bound, lower_bound=lower_bound)
    # run optimisation of paramters
    popt = dadi.Inference.optimize_log(p0=p0, data=fs_par, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=1e9, full_output=False)
    return popt # returns optimal nu and T
p_init
# get parameter combination with highest likelihood as starting values
p_init = np.array( df.sort_values(by='e_log ll', ascending=False).head().iloc[0, 2:4] )
# run optimisation with perturbed initial values on all 11 engines
fine_opt_res = dview.apply_async( run_dadi_on_par, p_init )
fine_opt_res.wall_time
ll = get_ll.map(fine_opt_res.get())
nu_opt = [p[0] for p in fine_opt_res.get()]
T_opt = [p[1] for p in fine_opt_res.get()]
df = pd.DataFrame({
        'a_nu' : nu_opt,
        'b_T' : T_opt,
        'c_ll' : ll
    })
df
Explanation: The parameter combination with the highest likelihood is not found several times from different starting values. This does not give me much confidence in the maximum likelihood parameter set.
End of explanation
# define function to run dadi with perturbed starting values
def run_dadi_optimisation(p_init, dadi_opt_func, sfs, maxiter=100): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to perturb and then run optimisation from
    dadi_opt_func:
        + optimize_log (BFGS based)
        + optimize_log_fmin (Nelder-Mead)
    sfs: ipyparallel Reference to dadi Spectrum object (previously imported into all engines)
    maxiter: maximum number of iterations
    
    # perturb starting values by up to a factor of 2 times 'fold'
    p0 = dadi.Misc.perturb_params(p_init, fold=2, upper_bound=upper_bound, lower_bound=lower_bound)
    # run optimisation of paramters
    popt = dadi_opt_func(p0=p0, data=sfs, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=maxiter, full_output=False)
    return popt # returns a tuple with starting values and optimal nu and T
from ipyparallel import Reference
sfs = Reference('fs_par')
# run optimisation with perturbed initial values on all 11 engines
# use Nelder-Mead algorithm this time
fine_opt_res = dview.apply_async( run_dadi_optimisation, \
                                 p_init, dadi.Inference.optimize_log_fmin, sfs, maxiter=100 )
fine_opt_res.wall_time
ll = get_ll.map(fine_opt_res.get())
nu_opt = [p[0] for p in fine_opt_res.get()]
T_opt = [p[1] for p in fine_opt_res.get()]
df = pd.DataFrame({
        'a_nu' : nu_opt,
        'b_T' : T_opt,
        'c_ll' : ll
    })
df
Explanation: None of these parameter combinations is the most likely determined above! I can therefore not be very confident that I have found the most likely parameter combination.
Maybe another optimisation algorithm is better at finding the global maximum.
End of explanation
# calculate the best-fit model SFS given the paramter combination with the highest likelihood
expected_sfs = func_ex(p_init, ns, pts_l)
print(p_init)
dadi.Plotting.plot_1d_comp_multinom(expected_sfs.fold()[:19], fs_par[:19], residual='linear')
Explanation: Again, none of these parameter combinations is close to the most likely ones determined above.
Let's see how the best-fitting model and data compare to each other.
End of explanation
# get optimally scaled (by theta) expected SFS for given data
model = dadi.Inference.optimally_scaled_sfs(expected_sfs.fold(), fs_par)
print(p_init)
plt.plot(model, 'ro-', label="model")
plt.plot(fs_par[:19], 'bo-', label="data")
plt.legend()
plt.xlabel('minor allele frequency')
plt.ylabel('SNP count')
Explanation: These logarithmic plots are not very helpfull.
End of explanation
%time run_dadi.map(range(20))
Explanation: So, if I could believe the optimal parameter combination found here to be the global maximum, it would imply that the ancestral populations of the parallelus sample have undergone an exponential decline to about 1/1000th of the ancient population size, $N_{ref}$ (i. e. very strong bottleneck) and this happend in the very recent past, i. e. only $2.3 \times 10^{-5} \times 2N_{ref}$ generations ago.
End of explanation
%time
# import 1D spectrum of ery on all engines:
fs_ery = dadi.Spectrum.from_file('dadiExercises/ERY.FOLDED.sfs.dadi_format')
fs_ery.pop_ids = ['ery']
ns = fs_ery.sample_sizes
# create link to function that specifies a simple growth or decline model
func = dadi.Demographics1D.growth
# create extrapolating version of the function
func_ex = dadi.Numerics.make_extrap_log_func(func)
# set lower and upper bounds to nu and T
upper_bound = [100, 3]
lower_bound = [1e-2, 0]
# set starting value
p0 = [1, 1] # corresponds to constant population size
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
# perturb starting values by up to a factor of 2 times fold:
p1 = dadi.Misc.perturb_params(p0, fold=1.5, upper_bound=upper_bound, lower_bound=lower_bound)
# run optimisation of paramters
popt = dadi.Inference.optimize_log(p0=p1, data=fs_ery, model_func=func_ex, pts=pts_l, \
                               lower_bound=lower_bound, upper_bound=upper_bound, \
                               verbose=0, maxiter=100, full_output=False)
popt
%time
# perturb starting values by up to a factor of 2 times fold:
p1 = dadi.Misc.perturb_params(p0, fold=1.5, upper_bound=upper_bound, lower_bound=lower_bound)
# run optimisation of paramters
popt = dadi.Inference.optimize_log(p0=p1, data=fs_ery, model_func=func_ex, pts=pts_l, \
                               lower_bound=lower_bound, upper_bound=upper_bound, \
                               verbose=0, maxiter=100, full_output=False)
popt
Explanation: Most of the time this just runs on a single core! I don't know why. Note the difference between the wall time and the total time.
Running the simulation locally on a single core:
End of explanation
%px %who
Explanation: It really shouldn't take 2 minutes to run this 20 times with different starting values.
End of explanation
?run_dadi.map
lbview.block
Explanation: Note, that p1 and popt are not in the namespace of the remote engines. They apparently have been cleared automatically.
End of explanation
cl.metadata
# create serial function
def serial_run_dadi(x): # for the function to be called with map, it needs to have one input variable
    # perturb starting values by up to a factor of 2 times fold:
    p1 = dadi.Misc.perturb_params(p0, fold=1.5, upper_bound=upper_bound, lower_bound=lower_bound)
    # run optimisation of paramters
    popt = dadi.Inference.optimize_log(p0=p1, data=fs_ery, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=100, full_output=False)
    return popt # returns a numpy array of optimal nu and T
# run dadi optimisation in non blocking mode
ar = lbview.map(serial_run_dadi, range(20))
ar.wait_interactive()
# the sum of the computation time of all of the tasks:
ar.serial_time
# the time between the first task submitted and last result received. 
# This is the actual cost of computation, including IPython overhead
ar.wall_time
parallel_speedup = ar.serial_time / ar.wall_time
parallel_speedup
ar.get()
for i, r in enumerate(ar):
    print "Optimsation number %i: %.3f (nu), %.3f (T)" % (i, r[0], r[1])
# create serial function
def serial_run_dadi_on_data(sfs):
    
    sfs: Spectrum object, provided to 'data' option of dadi.Inference.optimize_log
    
    # perturb starting values by up to a factor of 2 times fold:
    p1 = dadi.Misc.perturb_params(p0, fold=1.5, upper_bound=upper_bound, lower_bound=lower_bound)
    # run optimisation of paramters
    popt = dadi.Inference.optimize_log(p0=p1, data=sfs, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=0, maxiter=100, full_output=False)
    return popt # returns a numpy array of optimal nu and T
from ipyparallel import Reference
# create a Reference to `fs_ery` on the engines.
data = Reference('fs_ery')
ar = lbview.map(serial_run_dadi_on_data, [data] * 20)
for i, r in enumerate(ar):
    print "Optimsation number %i: %.3f (nu), %.3f (T)" % (i, r[0], r[1])
# run parallel optimisation on 1D SFS of par
ar = lbview.map(serial_run_dadi_on_data, [data] * 20)
# print result
for i, r in enumerate(ar):
    print "Optimsation number %i: %.3f (nu), %.3f (T)" % (i, r[0], r[1])
# run parallel optimisation on 1D SFS of par
# return results as they become ready (unordered)
# 3 calls of the function per task (chunksize)
ar = lbview.map(serial_run_dadi_on_data, [data] * 50, ordered=False, chunksize=3)
# print result
for i, r in enumerate(ar):
    print "Optimsation number %i: %.3f (nu), %.3f (T)" % (i, r[0], r[1])
Explanation: Jobs submitted to the load balanced view run asynchonously in non-blocking mode.
End of explanation
numpy.mean(ar[0])
Explanation: I think it should be a good idea to use load balanced views of the engines, since I think some optimisations take longer than others.
End of explanation
%who
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
Explanation: Look at max_iter from dadi.Inference.optimize_log.
Restart analysis
Restart the kernel from the menu to clear the namespace!
End of explanation
%%px --local
# import 1D spectrum of ery on all engines:
fs_ery = dadi.Spectrum.from_file('ERY.unfolded.sfs').fold()
# import 1D spectrum of ery on all engines:
fs_par = dadi.Spectrum.from_file('PAR.unfolded.sfs').fold()
cl[0]['fs_ery']
%%px --local
ns = fs_ery.sample_sizes # both populations have the same sample size
fs_ery.pop_ids = ['ery']
fs_par.pop_ids = ['par']
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
%%px --local
# create link to function that specifies the model
func = dadi.Demographics1D.growth
# create extrapolating version of the model function
func_ex = dadi.Numerics.make_extrap_log_func(func)
%pdoc dadi.Demographics1D.growth
# create load balanced view of engines
lbview = cl.load_balanced_view()
def run_dadi(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    if perturb == True:
        p_init = dadi.Misc.perturb_params(p_init, fold=fold, 
                                      upper_bound=upper_bound, lower_bound=lower_bound)
        # note upper_bound and lower_bound variables are expected to be in the namespace of each engine
    # run optimisation of paramters
    popt = dadi_opt_func(p0=p_init, data=sfs, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=verbose, maxiter=maxiter, full_output=full_output)
    return p_init, popt
Explanation: Use unfolded ANGSD spectra and fold with dadi:
End of explanation
%%px
# set up global variables on engines required for run_dadi function call
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
sfs = fs_ery # use ERY spectrum
perturb = False
fold = 1
maxiter = 10 # run a maximum of 20 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
Explanation: Ery
The following sets global variables on the remote engines that are called by run_dadi:
End of explanation
%%px
# set lower and upper bounds to nu and T
upper_bound = [1e4, 4]
lower_bound = [1e-4, 0]
p0_nu = np.logspace(-3, 3, base=10.0, num=15)
p0_T = np.logspace(-3, np.log10(3), base=10, num=15)
from itertools import product
# DO NOT RUN
#ar_ery = lbview.map(run_dadi, product(p0_nu, p0_T), block=False, order=True)
ar_ery.progress
ar_ery.get()[0]
Explanation: This lets the optimisations run with the faster BFGS algorithm and stop after a maximum of 10 iterations. This is intended for a broad sweep across the parameter space. Also note that this sets the optimisations to fit the ery spectrum.
End of explanation
%%px --local
import dill # better version of pickle
%%px --local
# set outname stub on all engines
outname = "OUT_exp_growth_model/ERY"
print outname
i=0
for p_init in product(p0_nu, p0_T):
    i+=1
    if i > 3: break
    name = outname[:]
    for p in p_init:
        if p < 1:
            name += "_%.4f" % (p)
        else:
            name += "_%4d" % (p)
    print name
Explanation: I needed to kill a single very long lasting optimisation job from the terminal (because it is not possible to do that with ipyparallel). Unfortunately that also means that the results from all 224 completed (or terminated) optimisation runs are lost as well. This is the major shortcoming of ipyparallel at the moment.
Couldn't I make ipyparallel pickle the results of the optimisations to file?
End of explanation
def run_dadi(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    if perturb == True:
        p_init = dadi.Misc.perturb_params(p_init, fold=fold, 
                                      upper_bound=upper_bound, lower_bound=lower_bound)
        # note upper_bound and lower_bound variables are expected to be in the namespace of each engine
    # run optimisation of paramters
    popt = dadi_opt_func(p0=p_init, data=sfs, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=verbose, maxiter=maxiter, full_output=full_output)
    # pickle to file
    import dill
    name = outname[:] # make copy of file name stub!
    for p in p_init:
        name += "_%.4f" % (p)
    with open(name + ".dill", "w") as fh:
        dill.dump((p_init, popt), fh)
    
    return p_init, popt
%ll
ar_ery = lbview.map(run_dadi, product(p0_nu, p0_T), block=False, order=True)
ar_ery.progress
Explanation: Now, define a new run_dadi function that pickles output to file:
End of explanation
import glob
len( glob.glob('OUT_exp_growth_model/*dill') )
ar_ery = []
for filename in glob.glob('OUT_exp_growth_model/*dill'):
    ar_ery.append( dill.load(open(filename)) )
ar_ery[1]
Explanation: I killed the last optimisation.
Load pickled output:
End of explanation
def get_flag_count(out, NM=True):
    
    out: list of tuples, each containing p_init and popt + additional info, including warnflags
    as produced by run_dadi.py
    
    from collections import defaultdict
    
    if NM: # if ar from Nelder-Mead
        i = 4 # the warnflag is reported at index position 4 in the output array
    else: # ar from BFGS optimisation
        i = 6
    
    warnflag = defaultdict(int)
    for res in out:
        if res[1][i] == 1: # notice the change in indexing
            warnflag[1] +=1
        elif res[1][i] == 2:
            warnflag[2] += 1
        elif res[1][i] == 0:
            warnflag[0] += 1
        else:
            warnflag[999] +=1
    if NM:
        print "success", warnflag[0]
        print "Maximum number of function evaluations made.", warnflag[1]
        print "Maximum number of iterations reached.", warnflag[2]
        print "unknown flag", warnflag[999]
    else:
        print "success", warnflag[0]
        print "Maximum number of iterations exceeded.", warnflag[1]
        print "Gradient and/or function calls not changing.", warnflag[2]
        print "unknown flag", warnflag[999]
get_flag_count(ar_ery, NM=False)
Explanation: Get flag count from optimisations:
End of explanation
def flatten(array):
    
        Returns a list of flattened elements of every inner lists (or tuples)
        ****RECURSIVE****
    
    import numpy
    res = []
    for el in array:
        if isinstance(el, (list, tuple, numpy.ndarray)):
            res.extend(flatten(el))
            continue
        res.append(el)
    return list( res )
import pandas as pd
i = 6 # where to find flag, 6 for BFGS, 4 for Nelder-Mead
successfull_popt_ery = [flatten(out)[:5] for out in ar_ery if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt_ery, \
                  columns=['nu_0', 'T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True) # smaller is better
Explanation: There were 39 successfull optimisations.
Create table of results from successful optimisations:
End of explanation
%%px
lower_bound = [1e-4, 1e-6]
cl[0]['lower_bound']
%%px
# set up global variables on engines required for run_dadi function call
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
sfs = fs_ery # use ERY spectrum
perturb = True
fold = 1
maxiter = 10 # run a maximum of 20 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "OUT_exp_growth_model/ERY_perturb"
p_init = [1, 1] # this should a neutral model, i. e. constant population size
from itertools import repeat
ar_ery = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=True)
ar_ery.progress
get_flag_count(ar_ery, NM=False)
Explanation: These results are completely inconclusive. I think it would make sense to limit the search space for $T$ to something well above 0 and start a few optimisations from parameter values close to a neutral model.
End of explanation
%ll OUT_exp_growth_model/ERY_perturb*
Explanation: None of these optimisations was able to converge on optimal parameter values.
End of explanation
! ls -l OUT_exp_growth_model/ERY_perturb* | wc -l
Explanation: The perturbation has not been very effective.
End of explanation
%%px
# use Nelder-Mead algorithm
dadi_opt_func = dadi.Inference.optimize_log_fmin 
! rm -f OUT_exp_growth_model/ERY_perturb*
! ls -l OUT_exp_growth_model/ERY_perturb* | wc -l
ar_ery = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=False)
ar_ery.elapsed
get_flag_count(ar_ery, NM=True)
%%px
# increase perturbation
fold = 3
ar_ery = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=False)
ar_ery.progress
ar_ery.elapsed
get_flag_count(ar_ery, NM=True)
! rm -f OUT_exp_growth_model/ERY_perturb*
%%px
maxiter = 50
ar_ery = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=False)
ar_ery.progress
ar_ery.elapsed
get_flag_count(ar_ery, NM=True)
Explanation: Why are there only 19 output files? There should be 20! Two output filenames happened to be identical.
End of explanation
i = 4 # where to find flag, 6 for BFGS, 4 for Nelder-Mead
successfull_popt_ery = [flatten(out)[:5] for out in ar_ery if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt_ery, \
                  columns=['nu_0', 'T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True) # smaller is better
Explanation: OK. Three successful optimisations.
End of explanation
p_init = np.array( df.sort_values(by='-logL').iloc[0, 2:4] )
p_init
ar_ery = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=False)
ar_ery.elapsed
get_flag_count(ar_ery, NM=True)
i = 4 # where to find flag, 6 for BFGS, 4 for Nelder-Mead
successfull_popt_ery = [flatten(out)[:5] for out in ar_ery if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt_ery, \
                  columns=['nu_0', 'T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True) # smaller is better
Explanation: OK, that looks promising. Two quite different parameter combinations have lead to the same optimal parameter combination.
End of explanation
%%px
maxiter = 100
ar_ery = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=False)
ar_ery.elapsed
Explanation: This looks like convergence.
End of explanation
get_flag_count(ar_ery, NM=True)
i = 4 # where to find flag, 6 for BFGS, 4 for Nelder-Mead
successfull_popt_ery = [flatten(out)[:5] for out in ar_ery if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt_ery, \
                  columns=['nu_0', 'T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True) # smaller is better
# safe optimisation result to file
# dill can dump but not load an AsyncResult object
dill.dump(list(ar_ery.get()), open("OUT_exp_growth_model/ERY_perturb_ar_ery.dill", "w"))
test = dill.load(open("OUT_exp_growth_model/ERY_perturb_ar_ery.dill"))
get_flag_count(test, NM=True)
Explanation: This is still extremely fast!
End of explanation
%%px
# set up global variables on engines required for run_dadi function call
dadi_opt_func = dadi.Inference.optimize_log_fmin # uses Nelder-Mead algorithm
sfs = fs_par # use PAR spectrum
perturb = True
fold = 3
maxiter = 50 # run a maximum of 50 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "OUT_exp_growth_model/PAR_perturb"
p_init = [1, 1]
ar_par = lbview.map(run_dadi, repeat(p_init, 20), block=False, order=False)
ar_par.elapsed/60
ar_par.progress
Explanation: Par
Since starting optimisations from a broad range of initial values did not lead to useful results with the ery spectrum, I am going to use the same strategy that led to a successful optimisation with ery: creating initial parameter values from a perturbation of parameters for a neutral model.
End of explanation
ar_par = []
for filename in glob.glob('OUT_exp_growth_model/PAR_perturb*'):
    ar_par.append( dill.load(open(filename)) )
get_flag_count(ar_par, NM=True)
Explanation: I needed to kill the last two processes in the terminal because they ran more than 15 minutes. I cannot use the AsyncResult object because of that and unstead need to load the pickled output from file.
End of explanation
i = 4 # where to find flag, 6 for BFGS, 4 for Nelder-Mead
successfull_popt_par = [flatten(out)[:5] for out in ar_par if out[1][i] == 0]
df = pd.DataFrame(data=successfull_popt_par, \
                  columns=['nu_0', 'T_0', 'nu_opt', 'T_opt', '-logL'])
df.sort_values(by='-logL', ascending=True) # smaller is better
Explanation: Only 1 successful optimisation.
End of explanation
%%px
# set up global variables on engines required for run_dadi function call
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
sfs = fs_par # use PAR spectrum
perturb = False
fold = 3
maxiter = 10 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "OUT_exp_growth_model/PAR"
ar_par = lbview.map(run_dadi, product(p0_nu, p0_T), block=False, order=False)
ar_par.progress
ar_par.elapsed/60
Explanation: That is hitting the upper limit for the time parameter.
Maybe a sweep through the parameter space could help to find an optimum.
End of explanation
import glob, dill
ar_par = []
for filename in glob.glob('OUT_exp_growth_model/PAR_[0123456789]*'):
    try:
        ar_par.append( dill.load( open(filename) ) )
    except:
        pass
get_flag_count(ar_par, NM=False)
Explanation: Unfortunately, this is too slow. I need to abort this by stopping the cluster of engines.
End of explanation
%who
! rm -f OUT_exp_growth_model/PAR_[0123456789]*
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
%%px --local
# import 1D spectrum of ery on all engines:
fs_ery = dadi.Spectrum.from_file('ERY.unfolded.sfs').fold()
# import 1D spectrum of ery on all engines:
fs_par = dadi.Spectrum.from_file('PAR.unfolded.sfs').fold()
%%px --local
ns = fs_ery.sample_sizes # both populations have the same sample size
fs_ery.pop_ids = ['ery']
fs_par.pop_ids = ['par']
# setting the smallest grid size slightly larger than the largest population sample size (36)
pts_l = [40, 50, 60]
%%px --local
# create link to function that specifies the model
func = dadi.Demographics1D.growth
# create extrapolating version of the model function
func_ex = dadi.Numerics.make_extrap_log_func(func)
# create load balanced view of engines
lbview = cl.load_balanced_view()
%%px --local
import dill # better version of pickle
from itertools import product
def run_dadi(p_init): # for the function to be called with map, it needs to have one input variable
    
    p_init: initial parameter values to run optimisation from
    
    if perturb == True:
        p_init = dadi.Misc.perturb_params(p_init, fold=fold, 
                                      upper_bound=upper_bound, lower_bound=lower_bound)
        # note upper_bound and lower_bound variables are expected to be in the namespace of each engine
    # run optimisation of paramters
    popt = dadi_opt_func(p0=p_init, data=sfs, model_func=func_ex, pts=pts_l, \
                                   lower_bound=lower_bound, upper_bound=upper_bound, \
                                   verbose=verbose, maxiter=maxiter, full_output=full_output)
    import dill
    # pickle to file
    name = outname[:] # make copy of file name stub!
    for p in p_init:
        name += "_%.4f" % (p)
    with open(name + ".dill", "w") as fh:
        dill.dump((p_init, popt), fh)
    
    return p_init, popt
%%px
# set lower and upper bounds to nu and T
upper_bound = [1e4, 4]
lower_bound = [1e-4, 1e-6]
Explanation: After restarting the cluster, it is important to be EXTREMELY careful to load all required variables and modules on all engines again.
restart
End of explanation
p0_nu = np.logspace(-2, 2, base=10.0, num=10)
p0_T = np.logspace(-2, np.log10(2), base=10, num=10)
%%px
# set up global variables on engines required for run_dadi function call
dadi_opt_func = dadi.Inference.optimize_log # uses BFGS algorithm
sfs = fs_par # use PAR spectrum
perturb = False
fold = 3
maxiter = 10 # run a maximum of 10 iterations
verbose = 0
full_output = True # need to have full output to get the warnflags (see below)
outname = "OUT_exp_growth_model/PAR"
ar_par = lbview.map(run_dadi, product(p0_nu, p0_T), block=False, order=False)
def get_flag_count(out, NM=True):
    
    out: list of tuples, each containing p_init and popt + additional info, including warnflags
    as produced by run_dadi.py
    
    from collections import defaultdict
    
    if NM: # if ar from Nelder-Mead
        i = 4 # the warnflag is reported at index position 4 in the output array
    else: # ar from BFGS optimisation
        i = 6
    
    warnflag = defaultdict(int)
    for res in out:
        if res[1][i] == 1: # notice the change in indexing
            warnflag[1] +=1
        elif res[1][i] == 2:
            warnflag[2] += 1
        elif res[1][i] == 0:
            warnflag[0] += 1
        else:
            warnflag[999] +=1
    if NM:
        print "success", warnflag[0]
        print "Maximum number of function evaluations made.", warnflag[1]
        print "Maximum number of iterations reached.", warnflag[2]
        print "unknown flag", warnflag[999]
    else:
        print "success", warnflag[0]
        print "Maximum number of iterations exceeded.", warnflag[1]
        print "Gradient and/or function calls not changing.", warnflag[2]
        print "unknown flag", warnflag[999]
ar_par = []
for filename in glob.glob('OUT_exp_growth_model/PAR_[0123456789]*'):
    try:
        ar_par.append( dill.load( open(filename) ) )
    except:
        pass
get_flag_count(ar_par, NM=False)
Explanation: I am trying to reduce the initial parameter sweep.
End of explanation
<END_TASK>