Merge branch 'devel'

Author: terhorst

.travis.yml

@@ -14,8 +14,8 @@ install:
   - conda config --set always_yes yes --set changeps1 no
   - conda update -y conda python
   - conda info -a
-  - conda install conda-build anaconda-client git ruamel_yaml pip
-  - pip install git+https://github.com/conda/constructor@6cf619d47
+  - conda install conda-build anaconda-client git pip constructor
+  # - pip install git+https://github.com/conda/constructor@6cf619d47
   - conda build -c terhorst -c bioconda conda
   - conda build -c terhorst -c bioconda conda --output | tee .packages
   - conda/template.py conda/construct.tpl $VERSION > conda/construct.yaml

README.rst

@@ -20,9 +20,9 @@ Quick start guide
    should run this once for each independent contig in your dataset,
    producing one SMC++ output file per contig.
 
-3. Fit the model using estimate_::
+3. Fit the model using cv_::
 
-     $ smc++ estimate -o analysis/ 1.25e-8 out/example.chr*.smc.gz
+     $ smc++ cv -o analysis/ 1.25e-8 out/example.chr*.smc.gz
 
    The first mandatory argument, ``1.25e-8``, is the per-generation
    mutation rate. The remaining arguments are the data files generated
@@ -250,6 +250,9 @@ is::
 
     $ smc++ estimate <mutation rate> <data file> [<data file> ...]
 
+*Please note that, in contrast to earlier versions, the recommended way to 
+estimate size history is now via the cv_ command*.
+
 Required arguments
 ^^^^^^^^^^^^^^^^^^
 
@@ -275,6 +278,20 @@ Optional arguments
 A number of other arguments concerning technical aspects of the fitting
 procedure exist. To see them, pass the ``-h`` option to ``estimate``.
 
+cv
+--
+
+This command is similar to estimate_, with the difference that it uses
+cross-validation to obtain sensible model parameters for use during estimation.
+*As of version 1.15, this is the recommended way to run SMC++*. The syntax and
+options for this command are nearly identical to estimate_:
+
+    $ smc++ cv <mutation rate> <data file> [<data file> ...]
+
+The optional `--folds` parameter can be used to specify the number of folds
+used for performing `k`-fold cross validation. The default is `2` and should be
+set higher in cases where you have more data.
+
 plot
 ----
 
@@ -400,13 +417,12 @@ good estimates*.
   may need to experiment a bit.
 
 - ``--timepoints``: This command specifies the starting and ending time points
-  of the model.  It accepts either a comma-separated list of two numbers
-  `t1,tK` specifying the starting and ending time points of the model (in
-  generations), or the default setting `h`. If `h` is specified, SMC++ will use
-  an experimental heuristic to calculate the model time points points
+  of the model.  It accepts two numbers `t1 tK` specifying the starting and
+  ending time points of the model (in generations). If not specified, SMC++
+  will use an heuristic to calculate the model time points points
   automatically.
 
-- ``--regularization-penalty``: This parameter penalizes curvature in
+- ``--regularization-penalty``, ``-rp``: This parameter penalizes curvature in
   the estimated size history. The default value of this parameter is
   ``6.0``. Lower values of the penalty shrink the estimated
   size history towards a line. If your estimates exhibit too much
@@ -423,7 +439,7 @@ good estimates*.
 
 - ``--knots``: This parameter specifies the number of spline knots 
   used in the underlying representation of the size history. The default
-  value is ``32``. Using fewer knots can lead to smoother fits, however
+  value is ``8``. Using fewer knots can lead to smoother fits, however
   underspecifying this parameter may smooth out interesting features of
   the size history.
 

bench/bench.py

@@ -7,6 +7,7 @@ build:
 
 source:
   git_url: ../
+  git_tag: devel
 
 requirements:
   build:
@@ -20,7 +21,7 @@ requirements:
     - numpy >=1.10
     - gmp >=6.1.0
     - gsl >=2.2.1
-    - mpfr >=3.1.5
+    - mpfr >=4.0.0
     - setuptools_scm >=1.15
   run:
     - python
@@ -45,6 +46,7 @@ requirements:
     - ad >=1.3.2
     - readline >=6.2
     - attrs >=17.2.0
+    - msprime >=0.5.0
 
 about:
     home: https://github.com/popgenmethods/smcpp

bench/large.txt.gz

@@ -1,31 +1,47 @@
 #!/bin/bash -x
 SMC=$(which smc++)
+TMP=$(mktemp -d)
 set -e
-$SMC vcf2smc -v example/example.vcf.gz /tmp/example.1.smc.gz 1 msp1:msp_0
-$SMC vcf2smc -d msp_0 msp_0 example/example.vcf.gz /tmp/example.2.smc.gz 1 msp2:msp_0,msp_3,msp_4
-$SMC vcf2smc -d msp_1 msp_1 example/example.vcf.gz /tmp/example.12.smc.gz 1 msp1:msp_1,msp_2 msp2:msp_3,msp_4,msp_0
-$SMC estimate -o /tmp/out/1 --unfold --knots 5 --em-iterations 1 1.25e-8 /tmp/example.1.smc.gz
-$SMC estimate -o /tmp/out/1 --unfold --knots 5 --timepoints 33,1000 --em-iterations 1 1.25e-8 /tmp/example.1.smc.gz
-$SMC estimate -p 0.01 -r 1e-8 -o /tmp/out/2 --knots 5 --em-iterations 1 1.25e-8 /tmp/example.2.smc.gz
-$SMC split -o /tmp/out/split --em-iterations 1 \
-    /tmp/out/1/model.final.json \
-    /tmp/out/2/model.final.json \
-    /tmp/example.*.smc.gz
-$SMC split --polarization-error .02 -o /tmp/out/split --em-iterations 1 \
-    /tmp/out/1/model.final.json \
-    /tmp/out/2/model.final.json \
-    /tmp/example.*.smc.gz
-$SMC posterior /tmp/out/1/model.final.json \
-    /tmp/matrix.npz /tmp/example.1.smc.gz /tmp/example.1.smc.gz
+$SMC vcf2smc -v example/example.vcf.gz $TMP/example.1.smc.gz 1 msp1:msp_0
+$SMC vcf2smc -v example/example.vcf.gz $TMP/example.11.smc.gz 1 msp1:msp_1
+$SMC vcf2smc -d msp_0 msp_0 example/example.vcf.gz $TMP/example.2.smc.gz 1 msp2:msp_0,msp_3,msp_4
+$SMC vcf2smc -d msp_1 msp_1 example/example.vcf.gz $TMP/example.12.smc.gz 1 msp1:msp_1,msp_2 msp2:msp_3,msp_4,msp_0
+$SMC chunk 10 200000 $TMP/chunk.1. $TMP/example.1.smc.gz
+$SMC chunk 10 200000 $TMP/chunk.2. $TMP/example.2.smc.gz
+$SMC chunk 10 200000 $TMP/chunk.12. $TMP/example.12.smc.gz
+$SMC estimate --em-iterations 1 -o $TMP/out/1 --unfold --knots 5 --em-iterations 1 1.25e-8 $TMP/example.1.smc.gz
+$SMC estimate --em-iterations 1 -o $TMP/out/1 --unfold --knots 5 --timepoints 33 1000 --em-iterations 1 1.25e-8 $TMP/example.1.smc.gz
+$SMC estimate --em-iterations 1 -p 0.01 -r 1e-8 -o $TMP/out/2 --knots 5 --em-iterations 1 1.25e-8 $TMP/example.2.smc.gz
+$SMC cv --em-iterations 1 --folds 2 -o $TMP/out/cv --fold 0 1e-8 $TMP/example.1.smc.gz $TMP/example.11.smc.gz
+$SMC cv --em-iterations 1 --folds 2 -o $TMP/out/cv --fold 1 1e-8 $TMP/example.1.smc.gz $TMP/example.11.smc.gz
+$SMC cv --em-iterations 1 --folds 2 -o $TMP/out/cv 1e-8 $TMP/example.1.smc.gz $TMP/example.11.smc.gz
+$SMC cv --em-iterations 1 --folds 2 -o $TMP/out/cv1 1e-8 $TMP/example.1.smc.gz $TMP/example.11.smc.gz
+$SMC split -o $TMP/out/split --em-iterations 1 \
+    $TMP/out/1/model.final.json \
+    $TMP/out/2/model.final.json \
+    $TMP/example.*.smc.gz
+$SMC split -o $TMP/out/split --em-iterations 1 \
+    $TMP/out/1/model.final.json \
+    $TMP/out/2/model.final.json \
+    $TMP/chunk*
+$SMC split --polarization-error .02 -o $TMP/out/split --em-iterations 1 \
+    $TMP/out/1/model.final.json \
+    $TMP/out/2/model.final.json \
+    $TMP/example.*.smc.gz
+$SMC posterior $TMP/out/1/model.final.json \
+    $TMP/matrix.npz $TMP/example.1.smc.gz $TMP/example.1.smc.gz
+$SMC simulate $TMP/out/1/model.final.json 2 .01 $TMP/out/1/sim.vcf
+$SMC simulate $TMP/out/split/model.final.json 2 .01 $TMP/out/split/sim.vcf
 $SMC posterior --colorbar -v \
-    --heatmap /tmp/plot.png \
-    /tmp/out/1/model.final.json \
-    /tmp/matrix.npz /tmp/example.1.smc.gz
+    --heatmap $TMP/plot.png \
+    $TMP/out/1/model.final.json \
+    $TMP/matrix.npz $TMP/example.1.smc.gz
 $SMC posterior -v \
-    /tmp/out/split/model.final.json \
-    /tmp/matrix.npz \
-    /tmp/example.12.smc.gz
-$SMC plot -c -g 29 /tmp/1.png /tmp/out/1/model.final.json
-$SMC plot /tmp/2.pdf /tmp/out/2/model.final.json
-$SMC plot -c /tmp/12.png /tmp/out/split/model.final.json
-$SMC plot -c /tmp/all.pdf /tmp/out/*/model.final.json
+    $TMP/out/split/model.final.json \
+    $TMP/matrix.npz \
+    $TMP/example.12.smc.gz
+$SMC plot -c -g 29 $TMP/1.png $TMP/out/1/model.final.json
+$SMC plot $TMP/2.pdf $TMP/out/2/model.final.json
+$SMC plot -c $TMP/12.png $TMP/out/split/model.final.json
+$SMC plot -c $TMP/all.pdf $TMP/out/*/model.final.json
+

bench/small.txt.gz

@@ -53,6 +53,7 @@ std::set<block_key> bin_key<1>::run(
     else
     {
         ret.emplace(key);
+        /*
         if (true) // or undistinguish)
             for (int aa = 0; aa <= na(0); ++aa)
             {
@@ -63,7 +64,6 @@ std::set<block_key> bin_key<1>::run(
                         (bb <= key(2)))
                     ret.emplace(tmp);
             }
-        /*
         if (nb > 0 and ((double)b / (double)nb > cutoff))
             for (int bb = (int)(cutoff * nb); bb <= nb; ++bb)
             {

bench/tiny.txt.gz

@@ -62,7 +62,8 @@
             "matplotlib>=1.5",
             "pysam>=0.9",
             "attrs>=17.2.0",
-            "pandas"],
+            "pandas",
+            "msprime>=0.5.0"],
         entry_points = {
             'console_scripts': ['smc++ = smcpp.frontend.console:main'],
             }