I am able to run the program when polarized = False, but when setting polarized = True I get the following error.
/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/dadi-1.7.0-py2.7-macosx-10.9-x86_64.egg/dadi/Numerics.py:171: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use
arr[tuple(seq)]instead of
arr[seq]. In the future this will be interpreted as an array index,
arr[np.array(seq)]`, which will result either in an error or a different result.
Traceback (most recent call last):
File "dadi_Run_2D_Set.py", line 249, in
Optimize_Functions.Optimize_Routine(fs, pts, prefix, "asym_mig_twoepoch", Models_2D.asym_mig_twoepoch, rounds, 8, reps=reps, maxiters=maxiters, folds=folds, param_labels = "nu1, nu2, m12a, m21a, m12b, m21b, T1, T2")
File "/Users/peterstokes/Documents/Science/Blackman_Lab/dadi/dadi_data/Optimize_Functions.py", line 237, in Optimize_Routine
rep_results = collect_results(fs, sim_model, params_opt, roundrep)
File "/Users/peterstokes/Documents/Science/Blackman_Lab/dadi/dadi_data/Optimize_Functions.py", line 121, in collect_results
chi2 = numpy.sum((folded_sim_model - fs)**2/folded_sim_model)
File "", line 3, in sub
File "build/bdist.macosx-10.9-x86_64/egg/dadi/Spectrum_mod.py", line 1878, in _check_other_folding
ValueError: Cannot operate with a folded Spectrum and an unfolded one.
Do you see anything that would cause errors in the first bit of my script?
#===========================================================================
Import data to create joint-site frequency spectrum
#===========================================================================
#**************
snps = "/Users/peterstokes/Documents/Science/Blackman_Lab/dadi/dadi_data/python_scripts/angsd_joined_LRW_60_REF_ANC_included_ready.txt"
#Create python dictionary from snps file
dd = dadi.Misc.make_data_dict(snps)
#**************
#pop_ids is a list which should match the populations headers of your SNPs file columns
pop_ids=["WILD", "LANDRACE"]
#**************
#projection sizes, in ALLELES not individuals
proj = [14,14]
#Convert this dictionary into folded AFS object
#[polarized = False] creates folded spectrum object
fs = dadi.Spectrum.from_data_dict(dd, pop_ids=pop_ids, projections = proj, polarized = True)
#print some useful information about the afs or jsfs
print "\n\n============================================================================\nData for site frequency spectrum\n============================================================================\n"
print "projection", proj
print "sample sizes", fs.sample_sizes
sfs_sum = numpy.around(fs.S(), 2)
print "Sum of SFS = ", sfs_sum, '\n', '\n'
#================================================================================
Calling external 2D models from the Models_2D.py script
#================================================================================
'''
We will use a function from the Optimize_Functions.py script for our optimization routines:
Optimize_Routine(fs, pts, outfile, model_name, func, rounds, param_number, reps=None, maxiters=None, folds=None, in_params=None, in_upper=None, in_lower=None, param_labels=" ")
Mandatory Arguments =
fs: spectrum object name
pts: grid size for extrapolation, list of three values
outfile: prefix for output naming
model_name: a label to help label the output files; ex. "no_mig"
func: access the model function from within 'moments_Run_Optimizations.py' or from a separate python model script, ex. after importing Models_2D, calling Models_2D.no_mig
rounds: number of optimization rounds to perform
param_number: number of parameters in the model selected (can count in params line for the model)
Optional Arguments =
reps: a list of integers controlling the number of replicates in each of the optimization rounds
maxiters: a list of integers controlling the maxiter argument in each of the optimization rounds
folds: a list of integers controlling the fold argument when perturbing input parameter values
in_params: a list of parameter values
in_upper: a list of upper bound values
in_lower: a list of lower bound values
param_labels: list of labels for parameters that will be written to the output file to keep track of their order
Below, I give all the necessary information to call each model available in the
Models_2D.py script. I have set the optimization routine to be the same for each
model using the optional lists below, which are included as optional arguments for
each model. This particular configuration will run 4 rounds as follows:
Round1 - 10 replicates, maxiter = 3, fold = 3
Round2 - 20 replicates, maxiter = 5, fold = 2
Round3 - 30 replicates, maxiter = 10, fold = 2
Round4 - 40 replicates, maxiter = 15, fold = 1
If this script was run as is, each model would be called and optimized sequentially;
this could take a very long time. For your actual analyses, I strongly recommend
creating multiple scripts with only a few models each and running them
independently. It is also not a good idea to mix models from the Diversification Set
and the Island Set, as each was meant to be mutually exclusive.
'''
#create a prefix based on the population names to label the output files
#ex. Pop1_Pop2
prefix = "_".join(pop_ids)
#**************
#make sure to define your extrapolation grid size (based on your projections)
pts = [50,60,70]
#**************
#Set the number of rounds here
rounds = 4
#define the lists for optional arguments
#you can change these to alter the settings of the optimization routine
reps = [10,20,30,40]
maxiters = [3,5,10,15]
folds = [3,2,2,1]`
Originally posted by @pstokespmb in #3 (comment)