Table of Contents

In [2]:

    

import sys


    

In [3]:

    

x = 36


    

In [4]:

    

pts_l = [x+10, x+20, x+30]
pts_l


    

    
Out[4]:





[46, 56, 66]

In [5]:

    

?eval


    

In [6]:

    

import sys

sys.path.insert(0, '/home/claudius/Downloads/dadi')

import dadi


    

In [7]:

    

dadi_model = 'dadi.Demographics1D.bottlegrowth'


    

In [8]:

    

f = eval(dadi_model)


    

In [9]:

    

%pdoc f


    

In [17]:

    

p_init = [0] * len(pts_l)
p_init


    

    
Out[17]:





[0, 0, 0]

In [11]:

    

for i in range(len(pts_l)):
    p_init[i] = range(pts_l[i], pts_l[i]+10)


    

In [12]:

    

p_init


    

    
Out[12]:





[[46, 47, 48, 49, 50, 51, 52, 53, 54, 55],
 [56, 57, 58, 59, 60, 61, 62, 63, 64, 65],
 [66, 67, 68, 69, 70, 71, 72, 73, 74, 75]]

In [13]:

    

from itertools import product


    

In [14]:

    

?product


    

In [15]:

    

for i, comb in enumerate(product(*p_init)):
    print comb
    if i > 20: break


    

    




(46, 56, 66)
(46, 56, 67)
(46, 56, 68)
(46, 56, 69)
(46, 56, 70)
(46, 56, 71)
(46, 56, 72)
(46, 56, 73)
(46, 56, 74)
(46, 56, 75)
(46, 57, 66)
(46, 57, 67)
(46, 57, 68)
(46, 57, 69)
(46, 57, 70)
(46, 57, 71)
(46, 57, 72)
(46, 57, 73)
(46, 57, 74)
(46, 57, 75)
(46, 58, 66)
(46, 58, 67)

In [21]:

    

print "these are the starting values \'{}\'".format(p_init)


    

    




these are the starting values '[0, 0, 0]'

In [24]:

    

% ll dadiExercises/


    

    




total 33676
lrwxrwxrwx 1 claudius       53 Feb 17 15:37 ERY.FOLDED.sfs -> /data3/claudius/Big_Data/ANGSD/SFS/ERY/ERY.FOLDED.sfs
-rw-rw-r-- 1 claudius      499 Mar 24 14:04 ERY.FOLDED.sfs.dadi_format
-rw-rw-r-- 1 claudius      499 Mar 24 14:02 ERY.FOLDED.sfs.dadi_format~
lrwxrwxrwx 1 claudius       37 Feb 18 17:46 EryPar.unfolded.2dsfs -> ../../ANGSD/FST/EryPar.unfolded.2dsfs
-rw-rw-r-- 1 claudius    13051 Feb 18 19:00 EryPar.unfolded.2dsfs.dadi_format
-rw-rw-r-- 1 claudius    13051 Feb 18 18:31 EryPar.unfolded.2dsfs.dadi_format~
drwxrwxr-x 5 claudius     4096 Feb 17 13:45 examples/
-rw-rw-r-- 1 claudius   155251 Mar 22 12:37 example_YRI_CEU.ipynb
-rw-rw-r-- 1 claudius   619518 Mar 21 11:09 First_Steps_with_dadi.ipynb
-rw-rw-r-- 1 claudius     1012 Mar 16 09:54 new.bib
lrwxrwxrwx 1 claudius       53 Feb 17 15:37 PAR.FOLDED.sfs -> /data3/claudius/Big_Data/ANGSD/SFS/PAR/PAR.FOLDED.sfs
-rw-rw-r-- 1 claudius      486 Mar 24 20:08 PAR.FOLDED.sfs.dadi_format
-rw-rw-r-- 1 claudius      450 Mar 24 20:08 PAR.FOLDED.sfs.dadi_format~
-rw-rw-r-- 1 claudius       18 Mar 19 11:20 seedms
-rw-rw-r-- 1 claudius 33643874 Mar 19 11:20 test.msout

In [26]:

    

path_to_spectrum_file = "dadiExercises/ERY.FOLDED.sfs.dadi_format"
dadi_model = "dadi.Demographics1D.growth"
upper = [1e3, 4]
lower = [1e-3, 0]
p_init = [1, 1]
dadi_opt_func = "dadi.Inference.optimize_log"


    

In [30]:

    

cmd = "./run_dadi.py -p {} -m {} -u \'{}\' -l \'{}\' -i \'{}\' -d {}".format(path_to_spectrum_file, dadi_model, upper, lower, p_init, dadi_opt_func)
cmd


    

    
Out[30]:





"./run_dadi.py -p dadiExercises/ERY.FOLDED.sfs.dadi_format -m dadi.Demographics1D.growth -u '[1000.0, 4]' -l '[0.001, 0]' -i '[1, 1]' -d dadi.Inference.optimize_log"

In [31]:

    

import os


    

In [32]:

    

os.system(cmd)


    

    
Out[32]:





0

In [1]:

    

import pickle, glob


    

In [9]:

    

pickled = glob.glob('OUT_run_dadi/ERY*pickle')
sorted(pickled[:10])


    

    
Out[9]:





['OUT_run_dadi/ERY_three_epoch_NM_0.001000_0.015849_0.057708_0.057708.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.001000_0.251189_0.057708_0.057708.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.001000_0.251189_4_0.006931.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.001000_3_0.480450_0.480450.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.001000_63_0.000833_0.000833.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.001000_63_0.057708_0.000100.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.015849_0.001000_0.006931_0.000100.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.015849_3_0.000833_4.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.015849_3_4_0.480450.pickle',
 'OUT_run_dadi/ERY_three_epoch_NM_0.015849_63_0.000833_0.006931.pickle']

In [10]:

    

opt_out = []
for filename in glob.iglob('OUT_run_dadi/ERY*pickle'):
    opt_out.append( pickle.load(open(filename, "r")) )


    

In [4]:

    

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]


    

In [11]:

    

get_flag_count(opt_out, NM=True)


    

    




success 0
Maximum number of function evaluations made. 0
Maximum number of iterations reached. 414
unknown flag 0

In [ ]: