{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cbcaae60",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import seaborn as sns\n",
    "import matplotlib as mpl\n",
    "from matplotlib import pyplot as plt\n",
    "import statistics\n",
    "sns.reset_defaults()\n",
    "sns.set_context(\"paper\")\n",
    "tex_fonts = {\n",
    "    # Use LaTeX to write all text\n",
    "    \"text.usetex\": True,\n",
    "    \"font.family\": \"sans-serif\",\n",
    "    \"font.sans-serif\": \"Helvetica\",\n",
    "    # Use 10pt font in plots, to match 10pt font in document\n",
    "    \"axes.labelsize\": 14,\n",
    "    \"axes.titlesize\": 14,\n",
    "    \"font.size\": 14,\n",
    "    # Make the legend/label fonts a little smaller\n",
    "    \"legend.fontsize\": 12,\n",
    "    \"xtick.labelsize\": 12,\n",
    "    \"ytick.labelsize\": 12\n",
    "}\n",
    "\n",
    "plt.rcParams.update(tex_fonts)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7ffe3990",
   "metadata": {},
   "outputs": [],
   "source": [
    "# accuracy by k value\n",
    "# total number of reads: 1490410\n",
    "vals_k = list(range(2, 13))\n",
    "data_k = []\n",
    "\n",
    "for k in vals_k:\n",
    "    path = \"mapquik-\" + str(k) + \"-31.stats\"\n",
    "    mapped_q60 = 0\n",
    "    wrong_q60 = 0\n",
    "    mapped_q0 = 0\n",
    "    wrong_q0 = 0\n",
    "    with open(path) as fm:\n",
    "        for line in fm:\n",
    "            spl = line.strip().split('\\t')\n",
    "            #print(spl)\n",
    "            if mapped_q60 == 0:\n",
    "                mapped_q60 = int(spl[0])\n",
    "                wrong_q60 = int(spl[1])\n",
    "            else:\n",
    "                mapped_q0 = int(spl[0])\n",
    "                wrong_q0 = int(spl[1])\n",
    "    data_k += [(100.0*((mapped_q60 + mapped_q0) / 1490410), 100.0*(mapped_q60 / (mapped_q60 + mapped_q0)), 100.0*(1.0 - ((wrong_q60 / mapped_q60))))]\n",
    "\n",
    "k_df = pd.DataFrame([x for x in data_k],columns=(\"Reads mapped\",\"Reads mapped at Q$60$\",\"Q$60$ accuracy\"),index=vals_k)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "44f380ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "# total number of reads: 1490410\n",
    "vals_l = list(range(10, 32))\n",
    "data_l = []\n",
    "\n",
    "for l in vals_l:\n",
    "    path = \"mapquik-5-\" + str(l) + \".stats\"\n",
    "    mapped_q60 = 0\n",
    "    wrong_q60 = 0\n",
    "    mapped_q0 = 0\n",
    "    wrong_q0 = 0\n",
    "    with open(path) as fm:\n",
    "        for line in fm:\n",
    "            spl = line.strip().split('\\t')\n",
    "            #print(spl)\n",
    "            if mapped_q60 == 0:\n",
    "                mapped_q60 = int(spl[0])\n",
    "                wrong_q60 = int(spl[1])\n",
    "            else:\n",
    "                mapped_q0 = int(spl[0])\n",
    "                wrong_q0 = int(spl[1])\n",
    "    data_l += [(100.0*((mapped_q60 + mapped_q0) / 1490410), 100.0*(mapped_q60 / (mapped_q60 + mapped_q0)), 100.0*(1.0 - ((wrong_q60 / mapped_q60))))]\n",
    "\n",
    "l_df = pd.DataFrame([x for x in data_l],columns=(\"Reads mapped\",\"Reads mapped at Q$60$\",\"Q$60$ accuracy\"),index=vals_l)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "ecfc9677",
   "metadata": {},
   "outputs": [],
   "source": [
    "# total number of reads: 1490410\n",
    "vals_d = [x/1000.0 for x in list(range(5, 21))]\n",
    "data_d = []\n",
    "for d in vals_d:\n",
    "    path = \"mapquik-5-31-\" + str(round(d, 3)) + \".stats\"\n",
    "    mapped_q60 = 0\n",
    "    wrong_q60 = 0\n",
    "    mapped_q0 = 0\n",
    "    wrong_q0 = 0\n",
    "    with open(path) as fm:\n",
    "        for line in fm:\n",
    "            spl = line.strip().split('\\t')\n",
    "            #print(spl)\n",
    "            if mapped_q60 == 0:\n",
    "                mapped_q60 = int(spl[0])\n",
    "                wrong_q60 = int(spl[1])\n",
    "            else:\n",
    "                mapped_q0 = int(spl[0])\n",
    "                wrong_q0 = int(spl[1])\n",
    "    data_d += [(100.0*((mapped_q60 + mapped_q0) / 1490410), 100.0*(mapped_q60 / (mapped_q60 + mapped_q0)), 100.0*(1.0 - ((wrong_q60 / mapped_q60))))]\n",
    "\n",
    "d_df = pd.DataFrame([x for x in data_d],columns=(\"Reads mapped\",\"Reads mapped at Q$60$\",\"Q$60$ accuracy\"),index=vals_d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f7ab64a0",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def adjust(df, a, xlab):\n",
    "    p = df.plot(ax=a)\n",
    "    p.set(xlabel=xlab, ylabel='\\%')\n",
    "    p.set(ylim=(92.4, 100.4))\n",
    "    p.xaxis.set_label_coords(0.5, -0.1)\n",
    "\n",
    "def adjust_r(df, a):\n",
    "    ax = df.plot(ax = a, x=df.columns[0], y=df.columns[1], legend=False)\n",
    "    ax2 = ax.twinx()\n",
    "    df.plot(ax = ax2, x=df.columns[0], y=df.columns[2], legend=False, color='r')\n",
    "    ax.set(xlabel = \"\", ylim=(5, 35))\n",
    "    ax2.set(xlabel = \"\", ylim=(7, 13))\n",
    "    ax.xaxis.set_label_coords(0.5, -0.1)\n",
    "    handles, labels = ax2.get_legend_handles_labels()\n",
    "    handles_r, labels_r = ax.get_legend_handles_labels()\n",
    "    handles.extend(handles_r)\n",
    "    labels.extend(labels_r)\n",
    "    ax.legend(handles, labels, loc='lower left')\n",
    "    return ax, ax2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7fe60a43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 83540 characters]",
      "text/plain": [
       "<Figure size 1300x800 with 9 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(2, 3, figsize=(13, 8), sharey = \"row\")\n",
    "\n",
    "r_k_df = pd.read_csv('k_perf.csv')\n",
    "r_l_df = pd.read_csv('l_perf.csv')\n",
    "r_d_df = pd.read_csv('d_perf.csv')\n",
    "\n",
    "k_ax, k_ax2 = adjust_r(r_k_df, axes[0][0])\n",
    "l_ax, l_ax2 = adjust_r(r_l_df, axes[0][1])\n",
    "d_ax, d_ax2 = adjust_r(r_d_df, axes[0][2])\n",
    "\n",
    "adjust(k_df, axes[1][0], r'$k$-min-mer length ($k$)')\n",
    "adjust(l_df, axes[1][1], r'Minimizer length ($\\ell$)')\n",
    "adjust(d_df, axes[1][2], r'Minimizer density ($\\delta$)')\n",
    "\n",
    "k_ax.set(ylabel='Runtime (s)')\n",
    "d_ax2.set(ylabel='Memory (GB)')\n",
    "\n",
    "fig.align_xlabels()\n",
    "plt.show()\n",
    "fig.savefig(\"kld_analysis_all.pdf\", bbox_inches='tight', dpi = 300)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4e387bde",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
