<commit_before>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "from os import path\n",
    "import re\n",
    "from skbio import DistanceMatrix\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from kwipexpt import *\n",
    "%matplotlib inline\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "library(tidyr)\n",
    "library(dplyr, warn.conflicts=F, quietly=T)\n",
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Calculate performance of kWIP\n",
    "=============================\n",
    "\n",
    "The next bit of python code calculates the performance of kWIP against the distance between samples calulcated from the alignments of their genomes.\n",
    "\n",
    "This code caluclates spearman's $\\rho$ between the off-diagonal elements of the triagnular distance matrices.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "expts = list(map(lambda fp: path.basename(fp.rstrip('/')), glob('data/*/')))\n",
    "print(\"Number of replicate experiments:\", len(expts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def process_expt(expt):\n",
    "    expt_results = []\n",
    "    \n",
    "    def extract_info(filename):\n",
    "        return re.search(r'kwip/(\\d\\.?\\d*)x-(0\\.\\d+)-(wip|ip).dist', filename).groups()\n",
    "    \n",
    "    # dict of scale: distance matrix, populated as we go\n",
    "    truths = {}\n",
    "    \n",
    "    for distfile in glob(\"data/{}/kwip/*.dist\".format(expt)):\n",
    "        cov, scale, metric = extract_info(distfile)\n",
    "        if scale not in truths:\n",
    "            genome_dist_path = 'data/{ex}/all_genomes-{sc}.dist'.format(ex=expt, sc=scale)\n",
    "            truths[scale] = load_sample_matrix_to_runs(genome_dist_path)\n",
    "        exptmat = DistanceMatrix.read(distfile)\n",
    "        rho = distmat_corr(truths[scale], exptmat, stats.spearmanr).correlation\n",
    "        expt_results.append({\n",
    "            \"coverage\": cov,\n",
    "            \"scale\": scale,\n",
    "            \"metric\": metric,\n",
    "            \"rho\": rho,\n",
    "            \"seed\": expt,\n",
    "        })\n",
    "    return expt_results\n",
    "\n",
    "#process_expt('3662')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "results = []\n",
    "for res in map(process_expt, expts):\n",
    "    results.extend(res)\n",
    "results = pd.DataFrame(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistical analysis\n",
    "====================\n",
    "\n",
    "Is done is R, as that's easier.\n",
    "\n",
    "Below we see a summary and structure of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R -i results\n",
    "\n",
    "results$coverage = as.numeric(as.character(results$coverage))\n",
    "results$scale = as.numeric(as.character(results$scale))\n",
    "\n",
    "print(summary(results))\n",
    "str(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment design\n",
    "\n",
    "Below we see the design of the experiment in terms of the two major variables.\n",
    "\n",
    "We have a series (vertically) that, at 30x coverage, looks at the effect of genetic variation on performance. There is a second series that examines the effect of coverage at an average pairwise genetic distance of 0.001.\n",
    "\n",
    "There are 100 replicates for each data point, performed as a separate bootstrap across the random creation of the tree and sampling of reads etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "ggplot(results, aes(x=coverage, y=scale)) +\n",
    "    geom_point() +\n",
    "    scale_x_log10() +\n",
    "    scale_y_log10() +\n",
    "    theme_bw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Effect of Coverage\n",
    "------------------\n",
    "\n",
    "Here we show the spread of data across the 100 reps as boxplots per metric and covreage level.\n",
    "\n",
    "I note that the weighted product seems slightly more variable, particularly at higher coverage. Though the median is nearly always higher"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "dat = results %>%\n",
    "        filter(scale==0.001, coverage<=30) %>%\n",
    "        select(rho, metric, coverage)\n",
    "    \n",
    "dat$coverage = as.factor(dat$coverage)\n",
    "ggplot(dat, aes(x=coverage, y=rho, fill=metric)) +\n",
    "    geom_boxplot(aes(fill=metric))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# AND AGAIN WITHOUT SUBSETTING\n",
    "dat = results %>%\n",
    "        filter(scale==0.001) %>%\n",
    "        select(rho, metric, coverage)\n",
    "    \n",
    "dat$coverage = as.factor(dat$coverage)\n",
    "ggplot(dat, aes(x=coverage, y=rho, fill=metric)) +\n",
    "    geom_boxplot(aes(fill=metric)) +\n",
    "    theme_bw()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "dat = subset(results, scale==0.001, select=-scale)\n",
    "ggplot(dat, aes(x=coverage, y=rho, colour=seed, linetype=metric)) +\n",
    "    geom_line() +\n",
    "    scale_x_log10()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "summ = results %>%\n",
    "    filter(scale==0.001) %>%\n",
    "    select(-scale) %>%\n",
    "    group_by(coverage, metric) %>%\n",
    "    summarise(rho_av=mean(rho), rho_err=sd(rho))\n",
    " \n",
    "p = ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho, \" +- SD\"))) +\n",
    "    #scale_x_log10()+\n",
    "    #ggtitle(\"Performance of WIP & IP\") +\n",
    "    theme_bw()\n",
    "\n",
    "pdf(\"coverage-vs-rho_full.pdf\",width=7, height=4)\n",
    "print(p)\n",
    "dev.off()\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "summ = results %>%\n",
    "    filter(scale==0.001, coverage <= 50) %>%\n",
    "    select(-scale) %>%\n",
    "    group_by(coverage, metric) %>%\n",
    "    summarise(rho_av=mean(rho), rho_err=sd(rho))\n",
    " \n",
    "p = ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho, \" +- SD\"))) +\n",
    "    #scale_x_log10()+\n",
    "    #ggtitle(\"Performance of WIP & IP\") +\n",
    "    theme_bw()\n",
    "\n",
    "pdf(\"coverage-vs-rho_50x.pdf\",width=5, height=4)\n",
    "print(p)\n",
    "dev.off()\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "sem <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))\n",
    "summ = results %>%\n",
    "    filter(scale==0.001) %>%\n",
    "    select(-scale) %>%\n",
    "    group_by(coverage, metric) %>%\n",
    "    summarise(rho_av=mean(rho), rho_err=sem(rho))\n",
    " \n",
    "ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho))) +\n",
    "    scale_x_log10()+\n",
    "    theme_bw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "cov_diff = results %>%\n",
    "    filter(scale==0.001) %>%\n",
    "    select(rho, metric, coverage, seed) %>%\n",
    "    spread(metric, rho) %>%\n",
    "    mutate(diff=wip-ip) %>%\n",
    "    select(coverage, seed, diff) \n",
    "print(summary(cov_diff))\n",
    "\n",
    "p = ggplot(cov_diff, aes(x=coverage, y=diff, colour=seed)) +\n",
    "    geom_line() +\n",
    "    scale_x_log10() +\n",
    "    ggtitle(\"Per expt difference in performance (wip - ip)\")\n",
    "print(p)\n",
    "\n",
    "summ = cov_diff %>%\n",
    "    group_by(coverage) %>%\n",
    "    summarise(diff_av=mean(diff), diff_sd=sd(diff))\n",
    "\n",
    "ggplot(summ, aes(x=coverage, y=diff_av, ymin=diff_av-diff_sd, ymax=diff_av+diff_sd)) +\n",
    "    geom_line() + \n",
    "    geom_ribbon(alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Improvment in Spearman's \", rho, \" (wip - IP)\"))) +\n",
    "    scale_x_log10() +\n",
    "    theme_bw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "var = results %>%\n",
    "        filter(coverage == 10, scale <= 0.05) %>%\n",
    "        select(metric, rho, scale)\n",
    "var$scale = as.factor(as.character(var$scale))\n",
    "str(var)\n",
    "ggplot(var, aes(x=scale, y=rho, fill=metric)) +\n",
    "    geom_boxplot(aes(fill=metric)) +\n",
    "    xlab('Mean pairwise variation') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho))) +\n",
    "    theme_bw()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "summ = results %>%\n",
    "        filter(coverage == 10, scale <= 0.04) %>%\n",
    "        select(-coverage) %>%\n",
    "        group_by(scale, metric) %>%\n",
    "        summarise(rho_av=mean(rho), rho_sd=sd(rho))\n",
    "str(summ)\n",
    "p = ggplot(summ, aes(x=scale, y=rho_av, ymin=rho_av-rho_sd, ymax=rho_av+rho_sd, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab(expression(paste('Mean pairwise variation (', pi, ')'))) +\n",
    "    ylab(expression(paste(\"Spearman's \", rho, \" +- SD\"))) +\n",
    "    scale_x_log10()+\n",
    "    theme_bw()\n",
    "\n",
    "pdf(\"pi-vs-performance.pdf\",width=5, height=4)\n",
    "print(p)\n",
    "dev.off()\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "var_diff = results %>%\n",
    "    filter(coverage==10) %>%\n",
    "    select(rho, metric, scale, seed) %>%\n",
    "    spread(metric, rho) %>%\n",
    "    mutate(diff=wip-ip) %>%\n",
    "    select(scale, seed, diff) \n",
    "\n",
    "\n",
    "summ_var_diff = var_diff %>%\n",
    "    group_by(scale) %>%\n",
    "    summarise(diff_av=mean(diff), diff_sd=sd(diff))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "p = ggplot(var_diff, aes(x=scale, y=diff, colour=seed)) +\n",
    "    geom_line() +\n",
    "    scale_x_log10() +\n",
    "    ggtitle(\"Per expt difference in performance (wip - ip)\")\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "ggplot(summ_var_diff, aes(x=scale, y=diff_av, ymin=diff_av-diff_sd, ymax=diff_av+diff_sd)) +\n",
    "    geom_line() + \n",
    "    geom_ribbon(alpha=0.2) +\n",
    "    xlab('Average variants/site') +\n",
    "    ylab(expression(paste(\"Improvment in Spearman's \", rho, \" (wip - IP)\"))) +\n",
    "    scale_x_log10() +\n",
    "    theme_bw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2+"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}<commit_msg>add change to old notebook<commit_after>{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "from os import path\n",
    "import re\n",
    "from skbio import DistanceMatrix\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "from kwipexpt import *\n",
    "%matplotlib inline\n",
    "%load_ext rpy2.ipython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "library(tidyr)\n",
    "library(dplyr, warn.conflicts=F, quietly=T)\n",
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "Calculate performance of kWIP\n",
    "=============================\n",
    "\n",
    "The next bit of python code calculates the performance of kWIP against the distance between samples calulcated from the alignments of their genomes.\n",
    "\n",
    "This code caluclates spearman's $\\rho$ between the off-diagonal elements of the triagnular distance matrices.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "expts = list(map(lambda fp: path.basename(fp.rstrip('/')), glob('data/*/')))\n",
    "print(\"Number of replicate experiments:\", len(expts))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def process_expt(expt):\n",
    "    expt_results = []\n",
    "    \n",
    "    def extract_info(filename):\n",
    "        return re.search(r'kwip/(\\d\\.?\\d*)x-(0\\.\\d+)-(wip|ip).dist', filename).groups()\n",
    "    \n",
    "    # dict of scale: distance matrix, populated as we go\n",
    "    truths = {}\n",
    "    \n",
    "    for distfile in glob(\"data/{}/kwip/*.dist\".format(expt)):\n",
    "        cov, scale, metric = extract_info(distfile)\n",
    "        if scale not in truths:\n",
    "            genome_dist_path = 'data/{ex}/all_genomes-{sc}.dist'.format(ex=expt, sc=scale)\n",
    "            truths[scale] = load_sample_matrix_to_runs(genome_dist_path)\n",
    "        exptmat = DistanceMatrix.read(distfile)\n",
    "        rho = distmat_corr(truths[scale], exptmat, stats.spearmanr).correlation\n",
    "        expt_results.append({\n",
    "            \"coverage\": cov,\n",
    "            \"scale\": scale,\n",
    "            \"metric\": metric,\n",
    "            \"rho\": rho,\n",
    "            \"seed\": expt,\n",
    "        })\n",
    "    return expt_results\n",
    "\n",
    "#process_expt('3662')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "results = []\n",
    "for res in map(process_expt, expts):\n",
    "    results.extend(res)\n",
    "results = pd.DataFrame(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Statistical analysis\n",
    "====================\n",
    "\n",
    "Is done is R, as that's easier.\n",
    "\n",
    "Below we see a summary and structure of the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R -i results\n",
    "\n",
    "results$coverage = as.numeric(as.character(results$coverage))\n",
    "results$scale = as.numeric(as.character(results$scale))\n",
    "\n",
    "print(summary(results))\n",
    "str(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Experiment design\n",
    "\n",
    "Below we see the design of the experiment in terms of the two major variables.\n",
    "\n",
    "We have a series (vertically) that, at 30x coverage, looks at the effect of genetic variation on performance. There is a second series that examines the effect of coverage at an average pairwise genetic distance of 0.001.\n",
    "\n",
    "There are 100 replicates for each data point, performed as a separate bootstrap across the random creation of the tree and sampling of reads etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "ggplot(results, aes(x=coverage, y=scale)) +\n",
    "    geom_point() +\n",
    "    scale_x_log10() +\n",
    "    scale_y_log10() +\n",
    "    theme_bw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Effect of Coverage\n",
    "------------------\n",
    "\n",
    "Here we show the spread of data across the 100 reps as boxplots per metric and covreage level.\n",
    "\n",
    "I note that the weighted product seems slightly more variable, particularly at higher coverage. Though the median is nearly always higher"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "dat = results %>%\n",
    "        filter(scale==0.001, coverage<=30) %>%\n",
    "        select(rho, metric, coverage)\n",
    "    \n",
    "dat$coverage = as.factor(dat$coverage)\n",
    "ggplot(dat, aes(x=coverage, y=rho, fill=metric)) +\n",
    "    geom_boxplot(aes(fill=metric))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "# AND AGAIN WITHOUT SUBSETTING\n",
    "dat = results %>%\n",
    "        filter(scale==0.001) %>%\n",
    "        select(rho, metric, coverage)\n",
    "    \n",
    "dat$coverage = as.factor(dat$coverage)\n",
    "ggplot(dat, aes(x=coverage, y=rho, fill=metric)) +\n",
    "    geom_boxplot(aes(fill=metric)) +\n",
    "    theme_bw()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "\n",
    "dat = subset(results, scale==0.001, select=-scale)\n",
    "ggplot(dat, aes(x=coverage, y=rho, colour=seed, linetype=metric)) +\n",
    "    geom_line() +\n",
    "    scale_x_log10()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "summ = results %>%\n",
    "    filter(scale==0.001) %>%\n",
    "    select(-scale) %>%\n",
    "    group_by(coverage, metric) %>%\n",
    "    summarise(rho_av=mean(rho), rho_err=sd(rho))\n",
    " \n",
    "p = ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho, \" +- SD\"))) +\n",
    "    #scale_x_log10()+\n",
    "    #ggtitle(\"Performance of WIP & IP\") +\n",
    "    theme_bw()\n",
    "\n",
    "pdf(\"coverage-vs-rho_full.pdf\",width=7, height=4)\n",
    "print(p)\n",
    "dev.off()\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "summ = results %>%\n",
    "    filter(scale==0.001, coverage <= 50) %>%\n",
    "    select(-scale) %>%\n",
    "    group_by(coverage, metric) %>%\n",
    "    summarise(rho_av=mean(rho), rho_err=sd(rho))\n",
    " \n",
    "p = ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho, \" +- SD\"))) +\n",
    "    #scale_x_log10()+\n",
    "    #ggtitle(\"Performance of WIP & IP\") +\n",
    "    theme_bw()\n",
    "\n",
    "pdf(\"coverage-vs-rho_50x.pdf\",width=5, height=4)\n",
    "print(p)\n",
    "dev.off()\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "sem <- function(x) sqrt(var(x,na.rm=TRUE)/length(na.omit(x)))\n",
    "summ = results %>%\n",
    "    filter(scale==0.001) %>%\n",
    "    select(-scale) %>%\n",
    "    group_by(coverage, metric) %>%\n",
    "    summarise(rho_av=mean(rho), rho_err=sem(rho))\n",
    " \n",
    "ggplot(summ, aes(x=coverage, y=rho_av, ymin=rho_av-rho_err, ymax=rho_av+rho_err, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho))) +\n",
    "    scale_x_log10()+\n",
    "    theme_bw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "cov_diff = results %>%\n",
    "    filter(scale==0.001) %>%\n",
    "    select(rho, metric, coverage, seed) %>%\n",
    "    spread(metric, rho) %>%\n",
    "    mutate(diff=wip-ip) %>%\n",
    "    select(coverage, seed, diff) \n",
    "print(summary(cov_diff))\n",
    "\n",
    "p = ggplot(cov_diff, aes(x=coverage, y=diff, colour=seed)) +\n",
    "    geom_line() +\n",
    "    scale_x_log10() +\n",
    "    ggtitle(\"Per expt difference in performance (wip - ip)\")\n",
    "print(p)\n",
    "\n",
    "summ = cov_diff %>%\n",
    "    group_by(coverage) %>%\n",
    "    summarise(diff_av=mean(diff), diff_sd=sd(diff))\n",
    "\n",
    "ggplot(summ, aes(x=coverage, y=diff_av, ymin=diff_av-diff_sd, ymax=diff_av+diff_sd)) +\n",
    "    geom_line() + \n",
    "    geom_ribbon(alpha=0.2) +\n",
    "    xlab('Genome Coverage') +\n",
    "    ylab(expression(paste(\"Improvment in Spearman's \", rho, \" (wip - IP)\"))) +\n",
    "    scale_x_log10() +\n",
    "    theme_bw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "var = results %>%\n",
    "        filter(coverage == 10, scale <= 0.05) %>%\n",
    "        select(metric, rho, scale)\n",
    "var$scale = as.factor(as.character(var$scale))\n",
    "str(var)\n",
    "ggplot(var, aes(x=scale, y=rho, fill=metric)) +\n",
    "    geom_boxplot(aes(fill=metric)) +\n",
    "    xlab('Mean pairwise variation') +\n",
    "    ylab(expression(paste(\"Spearman's \", rho))) +\n",
    "    theme_bw()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "summ = results %>%\n",
    "        filter(coverage == 10, scale <= 0.04) %>%\n",
    "        select(-coverage) %>%\n",
    "        group_by(scale, metric) %>%\n",
    "        summarise(rho_av=mean(rho), rho_sd=sd(rho))\n",
    "str(summ)\n",
    "p = ggplot(summ, aes(x=scale, y=rho_av, ymin=rho_av-rho_sd, ymax=rho_av+rho_sd, group=metric)) +\n",
    "    geom_line(aes(linetype=metric)) +\n",
    "    geom_ribbon(aes(fill=metric), alpha=0.2) +\n",
    "    xlab(expression(paste('Mean pairwise variation (', pi, ')'))) +\n",
    "    ylab(expression(paste(\"Spearman's \", rho, \" +- SD\"))) +\n",
    "    scale_x_log10()+\n",
    "    theme_bw()\n",
    "\n",
    "pdf(\"pi-vs-performance.pdf\",width=5, height=4)\n",
    "print(p)\n",
    "dev.off()\n",
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "var_diff = results %>%\n",
    "    filter(coverage==10) %>%\n",
    "    select(rho, metric, scale, seed) %>%\n",
    "    spread(metric, rho) %>%\n",
    "    mutate(diff=wip-ip) %>%\n",
    "    select(scale, seed, diff) \n",
    "\n",
    "\n",
    "summ_var_diff = var_diff %>%\n",
    "    group_by(scale) %>%\n",
    "    summarise(diff_av=mean(diff), diff_sd=sd(diff))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "p = ggplot(var_diff, aes(x=scale, y=diff, colour=seed)) +\n",
    "    geom_line() +\n",
    "    scale_x_log10() +\n",
    "    ggtitle(\"Per expt difference in performance (wip - ip)\")\n",
    "print(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false,
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%R\n",
    "ggplot(summ_var_diff, aes(x=scale, y=diff_av, ymin=diff_av-diff_sd, ymax=diff_av+diff_sd)) +\n",
    "    geom_line() + \n",
    "    geom_ribbon(alpha=0.2) +\n",
    "    xlab('Average variants/site') +\n",
    "    ylab(expression(paste(\"Improvment in Spearman's \", rho, \" (wip - IP)\"))) +\n",
    "    scale_x_log10() +\n",
    "    theme_bw()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
{'fraction_non_alphanumeric': 0.22556629736786396, 'fraction_numerical': 0.008185731423329335, 'mean_word_length': 2.4943903341141658, 'pattern_counts': {'":': 346, '<': 13, '<?xml version=': 0, '>': 75, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '820680', 'n_tokens_mistral': 11479, 'n_tokens_neox': 11031, 'n_words': 2042}