tests/test_pairgbs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ipyrad testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.1.26\n"
     ]
    }
   ],
   "source": [
    "import ipyrad as ip      ## for RADseq assembly\n",
    "print ip.__version__     ## print version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## clear test directory if it already exists\n",
    "import shutil\n",
    "import os\n",
    "if os.path.exists(\"./test_pairgbs\"):\n",
    "    shutil.rmtree(\"./test_pairgbs\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting started -- Assembly objects\n",
    "The first step is to create an Assembly object. It takes an optional argument that provides it with an internal name. We could imagine that we planned to assemble and later combine data from multiple sequencing runs, but before combining them each group of samples has to be analyzed under a different set of parameters. As an example, we could call this data set \"2014_data_set\" and another \"2015_data_set\". "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  New Assembly: test_pairgbs\n"
     ]
    }
   ],
   "source": [
    "## create an Assembly object called data1. \n",
    "## It takes an 'test'\n",
    "data1 = ip.Assembly('test_pairgbs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data1.set_params(\"project_dir\", \"./test_pairgbs\")\n",
    "data1.set_params(\"raw_fastq_path\", \"./data/sim_pairgbs_*.gz\")\n",
    "data1.set_params(\"barcodes_path\", \"./data/sim_pairgbs_barcodes.txt\")\n",
    "data1.set_params(\"datatype\", \"pairgbs\")\n",
    "#data1.set_params(17, 0)\n",
    "#data1.set_params(19, 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1   project_dir           ./test_pairgbs                               \n",
      "  2   raw_fastq_path              ./data/sim_pairgbs_*.gz                      \n",
      "  3   barcodes_path               ./data/sim_pairgbs_barcodes.txt              \n",
      "  4   sorted_fastq_path                                                        \n",
      "  5   assembly_method             denovo                                       \n",
      "  6   reference_sequence                                                       \n",
      "  7   datatype                    pairgbs                                      \n",
      "  8   restriction_overhang        ('TGCAG', '')                                \n",
      "  9   max_low_qual_bases          5                                            \n",
      "  10  phred_Qscore_offset         33                                           \n",
      "  11  mindepth_statistical        6                                            \n",
      "  12  mindepth_majrule            6                                            \n",
      "  13  maxdepth                    1000                                         \n",
      "  14  clust_threshold             0.85                                         \n",
      "  15  max_barcode_mismatch        1                                            \n",
      "  16  filter_adapters             0                                            \n",
      "  17  filter_min_trim_len         35                                           \n",
      "  18  max_alleles_consens         2                                            \n",
      "  19  max_Ns_consens              (5, 5)                                       \n",
      "  20  max_Hs_consens              (8, 8)                                       \n",
      "  21  min_samples_locus           4                                            \n",
      "  22  max_SNPs_locus              (100, 100)                                   \n",
      "  23  max_Indels_locus            (5, 99)                                      \n",
      "  24  max_shared_Hs_locus         0.25                                         \n",
      "  25  edit_cutsites               (0, 0)                                       \n",
      "  26  trim_overhang               (1, 2, 2, 1)                                 \n",
      "  27  output_formats              *                                            \n",
      "  28  pop_assign_file                                                          \n",
      "  29  excludes                                                                 \n",
      "  30  outgroups                                                                \n"
     ]
    }
   ],
   "source": [
    "data1.get_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 1: Demultiplex the raw data files\n",
    "This demultiplexes, and links the new fastq data files as Samples to the Assembly object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     state  reads_raw\n",
      "1A0      1       4000\n",
      "1B0      1       4000\n",
      "1C0      1       4000\n",
      "1D0      1       4000\n",
      "2E0      1       4000\n",
      "2F0      1       4000\n",
      "2G0      1       4000\n",
      "2H0      1       4000\n",
      "3I0      1       4000\n",
      "3J0      1       4000\n",
      "3K0      1       4000\n",
      "3L0      1       4000\n"
     ]
    }
   ],
   "source": [
    "## demultiplex the raw_data files\n",
    "## set step1 to only go if no samples are present...\n",
    "data1.step1() #append=True)\n",
    "\n",
    "print data1.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#data1.set_params(4, \"./test_pairgbs/fastq/*\")\n",
    "#data1.link_fastqs()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2: Filter reads \n",
    "If for some reason we wanted to execute on just a subsample of our data, we could do this by selecting only certain samples to call the `step2` function on. Because `step2` is a function of `data`, it will always execute with the parameters that are linked to `data`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     state  reads_raw  reads_filtered\n",
      "1A0      2       4000            4000\n",
      "1B0      2       4000            4000\n",
      "1C0      2       4000            4000\n",
      "1D0      2       4000            4000\n",
      "2E0      2       4000            4000\n",
      "2F0      2       4000            4000\n",
      "2G0      2       4000            4000\n",
      "2H0      2       4000            4000\n",
      "3I0      2       4000            4000\n",
      "3J0      2       4000            4000\n",
      "3K0      2       4000            4000\n",
      "3L0      2       4000            4000\n"
     ]
    }
   ],
   "source": [
    "data1.step2(force=True) #[\"1A0\",\"1B0\"])#\n",
    "\n",
    "print data1.stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the `name` and `fname` of the `Sample` objects and edit them as desired without affecting the original data files. The `name` identifier is equal to the filename (`fname`) by default, but is the name used in the final output files, and thus it may be desirable to reduce it to something more readable, like below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  loading Assembly: test_pairgbs [test_pairgbs/test_pairgbs.assembly]\n"
     ]
    }
   ],
   "source": [
    "import ipyrad as ip\n",
    "data1 = ip.load.load_assembly(\"test_pairgbs/test_pairgbs.assembly\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "data1.step3()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     state  reads_raw  reads_filtered  reads_merged  clusters_total  \\\n",
      "1A0      3       4000            4000          4000             200   \n",
      "1B0      3       4000            4000          4000             200   \n",
      "1C0      3       4000            4000          4000             200   \n",
      "1D0      3       4000            4000          4000             200   \n",
      "2E0      3       4000            4000          4000             200   \n",
      "2F0      3       4000            4000          4000             200   \n",
      "2G0      3       4000            4000          4000             200   \n",
      "2H0      3       4000            4000          4000             200   \n",
      "3I0      3       4000            4000          4000             200   \n",
      "3J0      3       4000            4000          4000             200   \n",
      "3K0      3       4000            4000          4000             200   \n",
      "3L0      3       4000            4000          4000             200   \n",
      "\n",
      "     clusters_hidepth  \n",
      "1A0               200  \n",
      "1B0               200  \n",
      "1C0               200  \n",
      "1D0               200  \n",
      "2E0               200  \n",
      "2F0               200  \n",
      "2G0               200  \n",
      "2H0               200  \n",
      "3I0               200  \n",
      "3J0               200  \n",
      "3K0               200  \n",
      "3L0               200  \n"
     ]
    }
   ],
   "source": [
    "print data1.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "CompositeError",
     "evalue": "one or more exceptions from call to method: optim\n[3:apply]: ValueError: 'axis' entry is out of bounds\n[0:apply]: ValueError: 'axis' entry is out of bounds\n[1:apply]: ValueError: 'axis' entry is out of bounds\n[2:apply]: ValueError: 'axis' entry is out of bounds\n.... 8 more exceptions ...",
     "output_type": "error",
     "traceback": [
      "[3:apply]: ",
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)\u001b[1;32m<string>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36moptim\u001b[1;34m(args)\u001b[0m",
      "\u001b[0;32m    255\u001b[0m                                             \u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcounts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    256\u001b[0m                                              \u001b[0mdisp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 257\u001b[1;33m                                              full_output=False)",
      "\u001b[0m\u001b[0;32m    258\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msample\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhetero\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    259\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfmin\u001b[1;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback)\u001b[0m",
      "\u001b[0;32m    375\u001b[0m             'return_all': retall}",
      "\u001b[0;32m    376\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 377\u001b[1;33m     \u001b[0mres\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcallback\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mopts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    378\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    379\u001b[0m         \u001b[0mretlist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'x'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'fun'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nit'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nfev'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'status'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[1;34m(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)\u001b[0m",
      "\u001b[0;32m    433\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mretall\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    434\u001b[0m         \u001b[0mallvecs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 435\u001b[1;33m     \u001b[0mfsim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    436\u001b[0m     \u001b[0mnonzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.05\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    437\u001b[0m     \u001b[0mzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.00025\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[1;34m(*wrapper_args)\u001b[0m",
      "\u001b[0;32m    283\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    284\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 285\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    286\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    287\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mget_diploid_lik\u001b[1;34m(pstart, bfreqs, ustacks, counts)\u001b[0m",
      "\u001b[0;32m     94\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     95\u001b[0m         \u001b[1;31m## get likelihood for all sites\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 96\u001b[1;33m         \u001b[0mlik1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood1\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     97\u001b[0m         \u001b[0mlik2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     98\u001b[0m         \u001b[0mliks\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlik1\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mlik2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mlikelihood1\u001b[1;34m(errors, bfreqs, ustacks)\u001b[0m",
      "\u001b[0;32m     46\u001b[0m     \u001b[1;31m## make sure base_frequencies are in the right order\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     47\u001b[0m     \u001b[1;31m#print uniqstackl.sum()-uniqstack, uniqstackl.sum(), 0.001\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 48\u001b[1;33m     \u001b[0mtotals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     49\u001b[0m     \u001b[0mprob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbinom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtotals\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtotals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     50\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mprob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/numpy/core/_methods.pyc\u001b[0m in \u001b[0;36m_sum\u001b[1;34m(a, axis, dtype, out, keepdims)\u001b[0m",
      "\u001b[0;32m     30\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     31\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 32\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     33\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     34\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: 'axis' entry is out of bounds",
      "",
      "[0:apply]: ",
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)\u001b[1;32m<string>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36moptim\u001b[1;34m(args)\u001b[0m",
      "\u001b[0;32m    255\u001b[0m                                             \u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcounts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    256\u001b[0m                                              \u001b[0mdisp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 257\u001b[1;33m                                              full_output=False)",
      "\u001b[0m\u001b[0;32m    258\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msample\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhetero\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    259\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfmin\u001b[1;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback)\u001b[0m",
      "\u001b[0;32m    375\u001b[0m             'return_all': retall}",
      "\u001b[0;32m    376\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 377\u001b[1;33m     \u001b[0mres\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcallback\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mopts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    378\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    379\u001b[0m         \u001b[0mretlist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'x'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'fun'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nit'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nfev'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'status'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[1;34m(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)\u001b[0m",
      "\u001b[0;32m    433\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mretall\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    434\u001b[0m         \u001b[0mallvecs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 435\u001b[1;33m     \u001b[0mfsim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    436\u001b[0m     \u001b[0mnonzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.05\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    437\u001b[0m     \u001b[0mzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.00025\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[1;34m(*wrapper_args)\u001b[0m",
      "\u001b[0;32m    283\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    284\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 285\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    286\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    287\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mget_diploid_lik\u001b[1;34m(pstart, bfreqs, ustacks, counts)\u001b[0m",
      "\u001b[0;32m     94\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     95\u001b[0m         \u001b[1;31m## get likelihood for all sites\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 96\u001b[1;33m         \u001b[0mlik1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood1\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     97\u001b[0m         \u001b[0mlik2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     98\u001b[0m         \u001b[0mliks\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlik1\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mlik2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mlikelihood1\u001b[1;34m(errors, bfreqs, ustacks)\u001b[0m",
      "\u001b[0;32m     46\u001b[0m     \u001b[1;31m## make sure base_frequencies are in the right order\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     47\u001b[0m     \u001b[1;31m#print uniqstackl.sum()-uniqstack, uniqstackl.sum(), 0.001\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 48\u001b[1;33m     \u001b[0mtotals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     49\u001b[0m     \u001b[0mprob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbinom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtotals\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtotals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     50\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mprob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/numpy/core/_methods.pyc\u001b[0m in \u001b[0;36m_sum\u001b[1;34m(a, axis, dtype, out, keepdims)\u001b[0m",
      "\u001b[0;32m     30\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     31\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 32\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     33\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     34\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: 'axis' entry is out of bounds",
      "",
      "[1:apply]: ",
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)\u001b[1;32m<string>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36moptim\u001b[1;34m(args)\u001b[0m",
      "\u001b[0;32m    255\u001b[0m                                             \u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcounts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    256\u001b[0m                                              \u001b[0mdisp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 257\u001b[1;33m                                              full_output=False)",
      "\u001b[0m\u001b[0;32m    258\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msample\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhetero\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    259\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfmin\u001b[1;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback)\u001b[0m",
      "\u001b[0;32m    375\u001b[0m             'return_all': retall}",
      "\u001b[0;32m    376\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 377\u001b[1;33m     \u001b[0mres\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcallback\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mopts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    378\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    379\u001b[0m         \u001b[0mretlist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'x'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'fun'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nit'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nfev'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'status'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[1;34m(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)\u001b[0m",
      "\u001b[0;32m    433\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mretall\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    434\u001b[0m         \u001b[0mallvecs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 435\u001b[1;33m     \u001b[0mfsim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    436\u001b[0m     \u001b[0mnonzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.05\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    437\u001b[0m     \u001b[0mzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.00025\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[1;34m(*wrapper_args)\u001b[0m",
      "\u001b[0;32m    283\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    284\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 285\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    286\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    287\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mget_diploid_lik\u001b[1;34m(pstart, bfreqs, ustacks, counts)\u001b[0m",
      "\u001b[0;32m     94\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     95\u001b[0m         \u001b[1;31m## get likelihood for all sites\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 96\u001b[1;33m         \u001b[0mlik1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood1\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     97\u001b[0m         \u001b[0mlik2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     98\u001b[0m         \u001b[0mliks\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlik1\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mlik2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mlikelihood1\u001b[1;34m(errors, bfreqs, ustacks)\u001b[0m",
      "\u001b[0;32m     46\u001b[0m     \u001b[1;31m## make sure base_frequencies are in the right order\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     47\u001b[0m     \u001b[1;31m#print uniqstackl.sum()-uniqstack, uniqstackl.sum(), 0.001\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 48\u001b[1;33m     \u001b[0mtotals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     49\u001b[0m     \u001b[0mprob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbinom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtotals\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtotals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     50\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mprob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/numpy/core/_methods.pyc\u001b[0m in \u001b[0;36m_sum\u001b[1;34m(a, axis, dtype, out, keepdims)\u001b[0m",
      "\u001b[0;32m     30\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     31\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 32\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     33\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     34\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: 'axis' entry is out of bounds",
      "",
      "[2:apply]: ",
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m\u001b[1;31mValueError\u001b[0m                                Traceback (most recent call last)\u001b[1;32m<string>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m()\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36moptim\u001b[1;34m(args)\u001b[0m",
      "\u001b[0;32m    255\u001b[0m                                             \u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcounts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    256\u001b[0m                                              \u001b[0mdisp\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 257\u001b[1;33m                                              full_output=False)",
      "\u001b[0m\u001b[0;32m    258\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msample\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mhetero\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    259\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfmin\u001b[1;34m(func, x0, args, xtol, ftol, maxiter, maxfun, full_output, disp, retall, callback)\u001b[0m",
      "\u001b[0;32m    375\u001b[0m             'return_all': retall}",
      "\u001b[0;32m    376\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 377\u001b[1;33m     \u001b[0mres\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_minimize_neldermead\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcallback\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mcallback\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mopts\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    378\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mfull_output\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    379\u001b[0m         \u001b[0mretlist\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'x'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'fun'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nit'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'nfev'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mres\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;34m'status'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36m_minimize_neldermead\u001b[1;34m(func, x0, args, callback, xtol, ftol, maxiter, maxfev, disp, return_all, **unknown_options)\u001b[0m",
      "\u001b[0;32m    433\u001b[0m     \u001b[1;32mif\u001b[0m \u001b[0mretall\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    434\u001b[0m         \u001b[0mallvecs\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0msim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 435\u001b[1;33m     \u001b[0mfsim\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    436\u001b[0m     \u001b[0mnonzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.05\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    437\u001b[0m     \u001b[0mzdelt\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;36m0.00025\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/scipy/optimize/optimize.pyc\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[1;34m(*wrapper_args)\u001b[0m",
      "\u001b[0;32m    283\u001b[0m     \u001b[1;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    284\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m0\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m--> 285\u001b[1;33m         \u001b[1;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[1;33m+\u001b[0m \u001b[0margs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m    286\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m    287\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mget_diploid_lik\u001b[1;34m(pstart, bfreqs, ustacks, counts)\u001b[0m",
      "\u001b[0;32m     94\u001b[0m     \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     95\u001b[0m         \u001b[1;31m## get likelihood for all sites\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 96\u001b[1;33m         \u001b[0mlik1\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood1\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     97\u001b[0m         \u001b[0mlik2\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mhetero\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mlikelihood2\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merrors\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mbfreqs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mustacks\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     98\u001b[0m         \u001b[0mliks\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mlik1\u001b[0m\u001b[1;33m+\u001b[0m\u001b[0mlik2\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/Documents/ipyrad/ipyrad/assemble/jointestimate.pyc\u001b[0m in \u001b[0;36mlikelihood1\u001b[1;34m(errors, bfreqs, ustacks)\u001b[0m",
      "\u001b[0;32m     46\u001b[0m     \u001b[1;31m## make sure base_frequencies are in the right order\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     47\u001b[0m     \u001b[1;31m#print uniqstackl.sum()-uniqstack, uniqstackl.sum(), 0.001\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 48\u001b[1;33m     \u001b[0mtotals\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0marray\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m4\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     49\u001b[0m     \u001b[0mprob\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mscipy\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mstats\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbinom\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mpmf\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mtotals\u001b[0m\u001b[1;33m-\u001b[0m\u001b[0mustacks\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtotals\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0merrors\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     50\u001b[0m     \u001b[1;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mbfreqs\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mprob\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m/home/deren/anaconda/lib/python2.7/site-packages/numpy/core/_methods.pyc\u001b[0m in \u001b[0;36m_sum\u001b[1;34m(a, axis, dtype, out, keepdims)\u001b[0m",
      "\u001b[0;32m     30\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     31\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;32m---> 32\u001b[1;33m     \u001b[1;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0m\u001b[0;32m     33\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[0;32m     34\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mout\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mNone\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m",
      "\u001b[1;31mValueError\u001b[0m: 'axis' entry is out of bounds",
      "",
      "... 8 more exceptions ..."
     ]
    }
   ],
   "source": [
    "data1.step4()\n",
    "print data1.stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data1.step5([\"1A0\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "#data1.step5([\"1A0\"])  ## better filters for -N-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Quick parameter explanations are always on-hand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ip.get_params_info(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Log history \n",
    "A common problem after struggling through an analysis is that you find you've completely forgotten what parameters you used at what point, and when you changed them. The log history time stamps all calls to `set_params()`, as well as calls to `step` methods. It also records copies/branching of data objects.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "for i in data.log:\n",
    "    print i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Saving Assembly objects\n",
    "Assembly objects can be saved and loaded so that interactive analyses can be started, stopped, and returned to quite easily. The format of these saved files is a serialized 'pickle' object used by Python. Individual Sample objects are saved within Assembly objects. While it is important to remember that some of the information in Assembly objects is in their links to data files, most of the useful information that we would want to analyze post assembly is stored in the object itself. Thus these objects will be useful for making plots and tables of assembly statistics later. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## save assembly object\n",
    "#ip.save_assembly(\"data1.p\")\n",
    "\n",
    "## load assembly object\n",
    "#data = ip.load_assembly(\"data1.p\")\n",
    "#print data.name"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### figureing out insert"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "indels = []\n",
    "catg = []\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arcatg = np.array([[5, 0, 0, 0],\n",
    "                   [0, 5, 0, 0],\n",
    "                   [0, 0, 5, 0],\n",
    "                   [0, 0, 0, 5]],\n",
    "                  [[5, 0, 0, 0],\n",
    "                   [0, 5, 0, 0],\n",
    "                   [0, 0, 5, 0],\n",
    "                   [0, 0, 0, 5]], dtype=uint32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "arcatg"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}