delta_c14_plot.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "chinese-dodge",
   "metadata": {},
   "outputs": [],
   "source": [
    "import datetime\n",
    "from IPython.display import display\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import warnings\n",
    "import pickle\n",
    "%matplotlib inline\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "external-teach",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "import c14\n",
    "import c14.models.liver_new_mitosis as lm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "electric-ethiopia",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "presidential-alpha",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = pd.read_csv('C14data_liver_samples_20211020.csv')\n",
    "data = data.groupby(['type', 'sample', 'ploidy', 'pathology']).mean().dropna(how='all').reset_index()\n",
    "#data['Dbirth'] -= 0.5\n",
    "data['age'] = data['Dcoll'] - data['Dbirth']\n",
    "data = data.query('type == \"hepatocyte\" and pathology != [\"Y\", \"C\"]' )\n",
    "\n",
    "m2n = data['ploidy']==\"2n\"\n",
    "m4n = data['ploidy']==\"4n\"\n",
    "m2n4n = data['ploidy']==\"2n4n\"\n",
    "mC2n = data['ploidy']==\"C2n\"\n",
    "mC4n = data['ploidy']==\"C4n\"\n",
    "\n",
    "with pd.HDFStore('point_est.pandas') as st:\n",
    "    point_est = st['F_IS_spline_newploidy'] \n",
    "    \n",
    "    \n",
    "paraP =dict(point_est.loc['POP1spline','median'].dropna())\n",
    "mP = c14.models.liver_new_mitosis.POP1(dnatotal=lm.dnatotal_spline)\n",
    "mP.set_parameters_phy(paraP,mode='bayes')\n",
    "opP=c14.optimize(mP,data=c14.exp_data(data),step_size=0.1)\n",
    "simdataP=opP.calc_sim_data_dict(\\\n",
    "    mP.transform_physical_to_fit(paraP,mode='bayes'),mode='bayes')    \n",
    "    \n",
    "paraRm =dict(point_est.loc['Rmspline','median'].dropna())\n",
    "mRm = c14.models.liver_new_mitosis.Rm(dnatotal=lm.dnatotal_spline)\n",
    "mRm.set_parameters_phy(paraRm,mode='bayes')\n",
    "opRm=c14.optimize(mRm,data=c14.exp_data(data),step_size=0.1)\n",
    "simdataRm=opRm.calc_sim_data_dict(\\\n",
    "    mRm.transform_physical_to_fit(paraRm,mode='bayes'),mode='bayes')\n",
    "\n",
    "paraPOP3p_2x2n =dict(point_est.loc['POP3p_2x2n_NPspline','median'].dropna())\n",
    "mPOP3p_2x2n = c14.models.liver_new_mitosis.POP3p_2x2n_NP(dnatotal=lm.dnatotal_spline)\n",
    "mPOP3p_2x2n.set_parameters_phy(paraPOP3p_2x2n,mode='bayes')\n",
    "opPOP3p_2x2n=c14.optimize(mPOP3p_2x2n ,data=c14.exp_data(data),step_size=0.1)\n",
    "simdataPOP3p_2x2n=opPOP3p_2x2n.calc_sim_data_dict(\\\n",
    "    mPOP3p_2x2n.transform_physical_to_fit(paraPOP3p_2x2n,mode='bayes'),mode='bayes')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "random-treasurer",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.lines import Line2D\n",
    "from matplotlib import colors\n",
    "plt.figure(11,figsize=[6.4*3/1.5, 2*4.8/1.5])\n",
    "sh=(2,3)\n",
    "ax1 = plt.subplot2grid(sh,(0,0))\n",
    "ax2 = plt.subplot2grid(sh,(0,1))\n",
    "ax3 = plt.subplot2grid(sh,(0,2))\n",
    "ax4 = plt.subplot2grid(sh,(1,0))\n",
    "ax5 = plt.subplot2grid(sh,(1,1))\n",
    "plt.subplot2grid(sh,(1,2)).remove()\n",
    "\n",
    "\n",
    "\n",
    "cmap = colors.ListedColormap(['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown'])\n",
    "bounds=[2004,2012,2016,2017,2018,2019,2020]\n",
    "norm = colors.BoundaryNorm(bounds, cmap.N)\n",
    "sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
    "sm.set_array([])\n",
    "\n",
    "#divider = make_axes_locatable(ax3)\n",
    "#cax = divider.append_axes(position='right', size='2%', pad=0.1)\n",
    "#plt.colorbar(sm, cax=cax)\n",
    "#cax.set_ylabel('Collection date')\n",
    "marker1='o'\n",
    "marker2='x'\n",
    "size1=40\n",
    "size2=30\n",
    "\n",
    "ax1.set_title('2n')\n",
    "datam = data[m2n]\n",
    "simd = simdataP[m2n]\n",
    "ax=ax1\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax2.set_title('4n')   \n",
    "datam = data[m4n]\n",
    "simd = simdataP[m4n]\n",
    "ax=ax2\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax3.set_title('2n4n')    \n",
    "datam = data[m2n4n]\n",
    "simd = simdataP[m2n4n]\n",
    "ax=ax3\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "    \n",
    "ax4.set_title('C2n')   \n",
    "datam = data[mC2n]\n",
    "simd = simdataP[mC2n]\n",
    "ax=ax4\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax5.set_title('C4n')    \n",
    "datam = data[mC4n]\n",
    "simd = simdataP[mC4n]\n",
    "ax=ax5\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "    \n",
    "\n",
    "lims = np.array([ax1.get_ylim(),ax2.get_ylim(),ax3.get_ylim(),(-0.05,0.3)])\n",
    "yl = (lims.min(),lims.max())\n",
    "\n",
    "a = np.linspace(1930,2020,1000)\n",
    "ax1.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax2.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax3.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax4.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax5.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "#yl = (m.Catm.lin(a).min(),m.Catm.lin(a).max()*1.01)\n",
    "\n",
    "ax4.set_xlabel(r'Birth date')\n",
    "ax5.set_xlabel(r'Birth date')\n",
    "ax3.set_xlabel(r'Birth date')\n",
    "ax1.set_ylabel(r'$\\Delta14C$')\n",
    "ax4.set_ylabel(r'$\\Delta14C$')\n",
    "\n",
    "ax2.set_ylabel('')\n",
    "ax3.set_ylabel('')\n",
    "ax5.set_ylabel('')\n",
    "ax2.set_yticklabels('')\n",
    "ax3.set_yticklabels('')\n",
    "ax5.set_yticklabels('')\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout(w_pad=0)\n",
    "plt.subplots_adjust(wspace=0.018)\n",
    "ax1.set_ylim(yl)\n",
    "ax2.set_ylim(yl)\n",
    "ax3.set_ylim(yl)\n",
    "ax4.set_ylim(yl)\n",
    "ax5.set_ylim(yl)\n",
    "\n",
    "\n",
    "ax1.set_xlim(data.Dbirth.min()-5,2020)#data.Dbirth.max()+5+15)\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax3.set_xlim(ax1.get_xlim())\n",
    "ax4.set_xlim(ax1.get_xlim())\n",
    "ax5.set_xlim(ax1.get_xlim())\n",
    "custom_lines = [Line2D([0], [0], color='black',marker=marker1,lw=0 ),\n",
    "                Line2D([0], [0], color='black',marker=marker2,lw=0 ),\n",
    "                Line2D([0], [0], color='black',lw=0,alpha=1 ),\n",
    "                Line2D([0], [0], color=sm.to_rgba(2004),marker='s', lw=0),\n",
    "                Line2D([0], [0], color=sm.to_rgba(2012),marker='s', lw=0)]+\\\n",
    "                [Line2D([0], [0], color=sm.to_rgba(y),marker='s', lw=0) for y in range(2016,2020)]\n",
    "l=ax5.legend(custom_lines, ['Data','Model','Collection year','2004', '2012']+list(range(2016,2020)),\\\n",
    "          bbox_to_anchor=(1.1, 0.9),frameon=True,title='')\n",
    "t =l.get_texts()[2]\n",
    "t.set_position((-28,0))\n",
    "#plt.tight_layout()\n",
    "plt.savefig('model_data_comp_POP1_newploidy.pdf', bbox_inches='tight')\n",
    "plt.savefig('model_data_comp_POP1_newploidy.png', dpi=600,bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "entertaining-hygiene",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "civil-entity",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.lines import Line2D\n",
    "from matplotlib import colors\n",
    "plt.figure(11,figsize=[6.4*3/1.5, 2*4.8/1.5])\n",
    "sh=(2,3)\n",
    "ax1 = plt.subplot2grid(sh,(0,0))\n",
    "ax2 = plt.subplot2grid(sh,(0,1))\n",
    "ax3 = plt.subplot2grid(sh,(0,2))\n",
    "ax4 = plt.subplot2grid(sh,(1,0))\n",
    "ax5 = plt.subplot2grid(sh,(1,1))\n",
    "plt.subplot2grid(sh,(1,2)).remove()\n",
    "\n",
    "\n",
    "\n",
    "cmap = colors.ListedColormap(['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown'])\n",
    "bounds=[2004,2012,2016,2017,2018,2019,2020]\n",
    "norm = colors.BoundaryNorm(bounds, cmap.N)\n",
    "sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
    "sm.set_array([])\n",
    "\n",
    "#divider = make_axes_locatable(ax3)\n",
    "#cax = divider.append_axes(position='right', size='2%', pad=0.1)\n",
    "#plt.colorbar(sm, cax=cax)\n",
    "#cax.set_ylabel('Collection date')\n",
    "marker1='o'\n",
    "marker2='x'\n",
    "size1=40\n",
    "size2=30\n",
    "\n",
    "ax1.set_title('2n')\n",
    "datam = data[m2n]\n",
    "simd = simdataRm[m2n]\n",
    "ax=ax1\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax2.set_title('4n')   \n",
    "datam = data[m4n]\n",
    "simd = simdataRm[m4n]\n",
    "ax=ax2\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax3.set_title('2n4n')    \n",
    "datam = data[m2n4n]\n",
    "simd = simdataRm[m2n4n]\n",
    "ax=ax3\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "    \n",
    "ax4.set_title('C2n')   \n",
    "datam = data[mC2n]\n",
    "simd = simdataRm[mC2n]\n",
    "ax=ax4\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax5.set_title('C4n')    \n",
    "datam = data[mC4n]\n",
    "simd = simdataRm[mC4n]\n",
    "ax=ax5\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "    \n",
    "\n",
    "lims = np.array([ax1.get_ylim(),ax2.get_ylim(),ax3.get_ylim(),(-0.05,0.3)])\n",
    "yl = (lims.min(),lims.max())\n",
    "\n",
    "a = np.linspace(1930,2020,1000)\n",
    "ax1.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax2.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax3.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax4.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax5.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "#yl = (m.Catm.lin(a).min(),m.Catm.lin(a).max()*1.01)\n",
    "\n",
    "ax4.set_xlabel(r'Birth date')\n",
    "ax5.set_xlabel(r'Birth date')\n",
    "ax3.set_xlabel(r'Birth date')\n",
    "ax1.set_ylabel(r'$\\Delta14C$')\n",
    "ax4.set_ylabel(r'$\\Delta14C$')\n",
    "\n",
    "ax2.set_ylabel('')\n",
    "ax3.set_ylabel('')\n",
    "ax5.set_ylabel('')\n",
    "ax2.set_yticklabels('')\n",
    "ax3.set_yticklabels('')\n",
    "ax5.set_yticklabels('')\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout(w_pad=0)\n",
    "plt.subplots_adjust(wspace=0.018)\n",
    "ax1.set_ylim(yl)\n",
    "ax2.set_ylim(yl)\n",
    "ax3.set_ylim(yl)\n",
    "ax4.set_ylim(yl)\n",
    "ax5.set_ylim(yl)\n",
    "\n",
    "\n",
    "ax1.set_xlim(data.Dbirth.min()-5,2020)#data.Dbirth.max()+5+15)\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax3.set_xlim(ax1.get_xlim())\n",
    "ax4.set_xlim(ax1.get_xlim())\n",
    "ax5.set_xlim(ax1.get_xlim())\n",
    "custom_lines = [Line2D([0], [0], color='black',marker=marker1,lw=0 ),\n",
    "                Line2D([0], [0], color='black',marker=marker2,lw=0 ),\n",
    "                Line2D([0], [0], color='black',lw=0,alpha=1 ),\n",
    "                Line2D([0], [0], color=sm.to_rgba(2004),marker='s', lw=0),\n",
    "                Line2D([0], [0], color=sm.to_rgba(2012),marker='s', lw=0)]+\\\n",
    "                [Line2D([0], [0], color=sm.to_rgba(y),marker='s', lw=0) for y in range(2016,2020)]\n",
    "l=ax5.legend(custom_lines, ['Data','Model','Collection year','2004', '2012']+list(range(2016,2020)),\\\n",
    "          bbox_to_anchor=(1.1, 0.9),frameon=True,title='')\n",
    "t =l.get_texts()[2]\n",
    "t.set_position((-28,0))\n",
    "#plt.tight_layout()\n",
    "plt.savefig('model_data_comp_Rm_newploidy.pdf', bbox_inches='tight')\n",
    "plt.savefig('model_data_comp_Rm_newploidy.png', dpi=600,bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cubic-cookie",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.lines import Line2D\n",
    "from matplotlib import colors\n",
    "plt.figure(11,figsize=[6.4*3/1.5, 2*4.8/1.5])\n",
    "sh=(2,3)\n",
    "ax1 = plt.subplot2grid(sh,(0,0))\n",
    "ax2 = plt.subplot2grid(sh,(0,1))\n",
    "ax3 = plt.subplot2grid(sh,(0,2))\n",
    "ax4 = plt.subplot2grid(sh,(1,0))\n",
    "ax5 = plt.subplot2grid(sh,(1,1))\n",
    "plt.subplot2grid(sh,(1,2)).remove()\n",
    "\n",
    "\n",
    "\n",
    "cmap = colors.ListedColormap(['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown'])\n",
    "bounds=[2004,2012,2016,2017,2018,2019,2020]\n",
    "norm = colors.BoundaryNorm(bounds, cmap.N)\n",
    "sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)\n",
    "sm.set_array([])\n",
    "\n",
    "#divider = make_axes_locatable(ax3)\n",
    "#cax = divider.append_axes(position='right', size='2%', pad=0.1)\n",
    "#plt.colorbar(sm, cax=cax)\n",
    "#cax.set_ylabel('Collection date')\n",
    "marker1='o'\n",
    "marker2='x'\n",
    "size1=40\n",
    "size2=30\n",
    "\n",
    "ax1.set_title('2n')\n",
    "datam = data[m2n]\n",
    "simd = simdataPOP3p_2x2n[m2n]\n",
    "ax=ax1\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax2.set_title('4n')   \n",
    "datam = data[m4n]\n",
    "simd = simdataPOP3p_2x2n[m4n]\n",
    "ax=ax2\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax3.set_title('2n4n')    \n",
    "datam = data[m2n4n]\n",
    "simd = simdataPOP3p_2x2n[m2n4n]\n",
    "ax=ax3\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "    \n",
    "ax4.set_title('C2n')   \n",
    "datam = data[mC2n]\n",
    "simd = simdataPOP3p_2x2n[mC2n]\n",
    "ax=ax4\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "ax5.set_title('C4n')    \n",
    "datam = data[mC4n]\n",
    "simd = simdataPOP3p_2x2n[mC4n]\n",
    "ax=ax5\n",
    "ax.scatter(datam.Dbirth, datam.d14C, c=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker1,zorder=3,s=size1)\n",
    "ax.scatter(datam.Dbirth,simd , color=sm.to_rgba(datam.Dcoll), facecolor=sm.to_rgba(datam.Dcoll),\n",
    "             marker=marker2,zorder=2,s=size2)\n",
    "for index,row in datam.iterrows():\n",
    "    ax.plot([row.Dbirth,row.Dbirth],[row.d14C,simd[index]],':', color=sm.to_rgba(row.Dcoll),alpha=1,zorder=1)\n",
    "\n",
    "    \n",
    "\n",
    "lims = np.array([ax1.get_ylim(),ax2.get_ylim(),ax3.get_ylim(),(-0.05,0.3)])\n",
    "yl = (lims.min(),lims.max())\n",
    "\n",
    "a = np.linspace(1930,2020,1000)\n",
    "ax1.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax2.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax3.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax4.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "ax5.plot(a,c14.Catm().lin(a),':',color='gray')\n",
    "#yl = (m.Catm.lin(a).min(),m.Catm.lin(a).max()*1.01)\n",
    "\n",
    "\n",
    "\n",
    "ax2.set_ylabel('')\n",
    "ax3.set_ylabel('')\n",
    "ax5.set_ylabel('')\n",
    "ax2.set_yticklabels('')\n",
    "ax3.set_yticklabels('')\n",
    "ax5.set_yticklabels('')\n",
    "\n",
    "\n",
    "\n",
    "plt.tight_layout(w_pad=0)\n",
    "plt.subplots_adjust(wspace=0.018)\n",
    "ax1.set_ylim(yl)\n",
    "ax2.set_ylim(yl)\n",
    "ax3.set_ylim(yl)\n",
    "ax4.set_ylim(yl)\n",
    "ax5.set_ylim(yl)\n",
    "\n",
    "ax4.set_xlabel(r'Birth date')\n",
    "ax5.set_xlabel(r'Birth date')\n",
    "ax3.set_xlabel(r'Birth date')\n",
    "ax1.set_ylabel(r'$\\Delta14C$')\n",
    "ax4.set_ylabel(r'$\\Delta14C$')\n",
    "\n",
    "ax1.set_xlim(data.Dbirth.min()-5,2020)#data.Dbirth.max()+5+15)\n",
    "ax2.set_xlim(ax1.get_xlim())\n",
    "ax3.set_xlim(ax1.get_xlim())\n",
    "ax4.set_xlim(ax1.get_xlim())\n",
    "ax5.set_xlim(ax1.get_xlim())\n",
    "custom_lines = [Line2D([0], [0], color='black',marker=marker1,lw=0 ),\n",
    "                Line2D([0], [0], color='black',marker=marker2,lw=0 ),\n",
    "                Line2D([0], [0], color='black',lw=0,alpha=1 ),\n",
    "                Line2D([0], [0], color=sm.to_rgba(2004),marker='s', lw=0),\n",
    "                Line2D([0], [0], color=sm.to_rgba(2012),marker='s', lw=0)]+\\\n",
    "                [Line2D([0], [0], color=sm.to_rgba(y),marker='s', lw=0) for y in range(2016,2020)]\n",
    "l=ax5.legend(custom_lines, ['Data','Model','Collection year','2004', '2012']+list(range(2016,2020)),\\\n",
    "          bbox_to_anchor=(1.1, 0.9),frameon=True,title='')\n",
    "t =l.get_texts()[2]\n",
    "t.set_position((-28,0))\n",
    "#plt.tight_layout()\n",
    "plt.savefig('model_data_comp_POP3p_2x2n_newploidy.pdf', bbox_inches='tight')\n",
    "plt.savefig('model_data_comp_POP3p_2x2n_newploidy.png', dpi=600,bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "offshore-vermont",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "c14",
   "language": "python",
   "name": "c14"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}