generate_plots.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:43:54.497062Z",
     "start_time": "2021-07-31T13:43:53.802370Z"
    }
   },
   "outputs": [],
   "source": [
    "import warnings\n",
    "\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "from math import sqrt,log\n",
    "from iminuit import Minuit\n",
    "import pandas as pd\n",
    "idx = pd.IndexSlice\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1.inset_locator import InsetPosition,BboxConnector\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "from matplotlib.transforms import blended_transform_factory,TransformedBbox\n",
    "import matplotlib as mpl\n",
    "from matplotlib.ticker import MultipleLocator\n",
    "%matplotlib inline\n",
    "import seaborn as sb\n",
    "\n",
    "import clapy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:43:54.517770Z",
     "start_time": "2021-07-31T13:43:54.498632Z"
    }
   },
   "outputs": [],
   "source": [
    "def fitfunc( t, Tc, r, GF):\n",
    "    if t<(Tc-Tc*r):\n",
    "        res = GF/Tc*(t+Tc*r)\n",
    "    else:\n",
    "        res =  GF\n",
    "    return res \n",
    "\n",
    "def mypdfln( x,  mu,  std):\n",
    "    #cpdef means generate both c function and python function\n",
    "    norm = (sqrt(2*np.pi*std*std))\n",
    "    ret = (-1*(x-mu)*(x-mu)/(2.*std*std))-log(norm)\n",
    "    return ret\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "class Nowakowski_LH:  \n",
    "    def __init__(self,data,t):\n",
    "        self.data = data\n",
    "        self.t = t\n",
    "        self.ndata = len(data)\n",
    "    \n",
    "    def compute(self,  Tc, r, GF, s):\n",
    "        #this line is a cast not a copy. Let cython knows mydata will spit out double\n",
    "        mydata = self.data\n",
    "        myt = self.t\n",
    "        loglh = 0.\n",
    "        lmu = 0.\n",
    "        ler = 0.\n",
    "        for i in range(self.ndata):\n",
    "            lmu = fitfunc(myt[i],Tc,r,GF)\n",
    "            loglh -= mypdfln(mydata[i],lmu,s)\n",
    "        return loglh\n",
    "\n",
    "    \n",
    "def calc_statistic(dataframe,parameter):\n",
    "    tmp = dataframe.copy()\n",
    "    try:\n",
    "        tmp.rename({'GF_error' : 'GFf_error'}, axis=\"columns\",inplace=True)\n",
    "    except:\n",
    "        pass\n",
    "    for i in ['Tc','r','ss','sc','GFf']:\n",
    "        try:\n",
    "            tmp[i+'_2std'] = tmp[i+'_error'].mul(2)\n",
    "            tmp[i+'_mse'] = tmp[i].sub(1).pow(2)\n",
    "            tmp[i+'_d'] = tmp[i].sub(parameter[i])\n",
    "            tmp[i+'_p'] = tmp[i+'_d'].abs() < tmp[i+'_2std']\n",
    "            tmp[i+'_p2'] = np.logical_and( tmp[i+'_hpd46_l'] < parameter[i] ,parameter[i] < tmp[i+'_hpd46_u']    )\n",
    "        except:\n",
    "            print(i+' not found')\n",
    "    sdata = tmp.mean(level=['model','sPopulation','sSample','M','GF'])\n",
    "    for i in ['Tc','r','ss','sc','GFf']:\n",
    "        try:\n",
    "            sdata[i+'_mse'] = sdata[i+'_mse'].pow(0.5)\n",
    "        except:\n",
    "            print(i+' not found')\n",
    "    return tmp,sdata\n",
    "\n",
    "def ecdf(x):\n",
    "    ''' calc emerical cdf '''\n",
    "    xs = np.sort(x)\n",
    "    ys = np.arange(1, len(xs)+1)/float(len(xs))\n",
    "    return xs, ys\n",
    "\n",
    "def idx_name_levels(df,name):\n",
    "    ''' get levels from level 'name' '''\n",
    "    pos = df.index.names.index(name)\n",
    "    return df.index.levels[pos]\n",
    "\n",
    "def idx_name_labels(df,name):\n",
    "    ''' get labels from level 'name' '''\n",
    "    pos = df.index.names.index(name)\n",
    "    return df.index.labels[pos]\n",
    "\n",
    "def idx_name(df,name):\n",
    "    ''' get name-level values from df '''\n",
    "    return idx_name_levels(df,name)[idx_name_labels(df,name)]\n",
    "\n",
    "\n",
    "\n",
    "def plot2d_pandas(df,column,levelx,levely,return_data=False,sym_limits=None,sym_scale=1,axis=None,**kwargs):\n",
    "    ''' plot 2d slice from 2 indexes of df '''\n",
    "    xindx = idx_name_levels(df,levelx)\n",
    "    yindx = idx_name_levels(df,levely)\n",
    "    datag = df[column].groupby(level=[levely,levelx])\n",
    "    if np.sum(datag.count() == 1) != len(yindx)*len(xindx):\n",
    "        warnings.warn(\"DataFrame has more than 2 dimension. Average is plotted\")\n",
    "    data = datag.mean().sort_index().values.reshape((len(yindx),len(xindx)))\n",
    "    dx = (xindx[1] - xindx[0])\n",
    "    dy = (yindx[1] - yindx[0])\n",
    "    ext = [xindx[0] - dx/2,xindx[-1] + dx/2,yindx[0] - dy/2,yindx[-1] + dy/2]\n",
    "    args = dict()\n",
    "    if sym_limits is not None:\n",
    "        limit = np.absolute((data - sym_limits)).max()*sym_scale\n",
    "        args.update({'vmin':sym_limits-limit,'vmax':sym_limits+limit})\n",
    "    if axis is not None:\n",
    "        im = axis.imshow(data,origin='lower',extent=ext,**args,**kwargs)\n",
    "    else:\n",
    "        im = plt.imshow(data,origin='lower',extent=ext,**args,**kwargs)\n",
    "\n",
    "    if return_data:\n",
    "        return im,data\n",
    "    return im\n",
    "\n",
    "def zoom_effect02(ax1, ax2,loc1a=3, loc2a=2, loc1b=4, loc2b=1, **kwargs):\n",
    "    \"\"\"\n",
    "    ax1 : the main axes\n",
    "    ax1 : the zoomed axes\n",
    "\n",
    "    Similar to zoom_effect01.  The xmin & xmax will be taken from the\n",
    "    ax1.viewLim.\n",
    "    \"\"\"\n",
    "\n",
    "    tt = ax1.transScale + (ax1.transLimits + ax2.transAxes)\n",
    "    trans = blended_transform_factory(ax2.transData, tt)\n",
    "\n",
    "    bbox1 = ax1.bbox\n",
    "    bbox2 = TransformedBbox(ax1.viewLim, trans)\n",
    "    c1 = BboxConnector(bbox1, bbox2, loc1=loc1a, loc2=loc2a, **kwargs)\n",
    "    c1.set_clip_on(False)\n",
    "    c2 = BboxConnector(bbox1, bbox2, loc1=loc1b, loc2=loc2b, **kwargs)\n",
    "    c2.set_clip_on(False)\n",
    "\n",
    "    ax2.add_patch(c1)\n",
    "    ax2.add_patch(c2)\n",
    "\n",
    "    return c1, c2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:43:54.532188Z",
     "start_time": "2021-07-31T13:43:54.520577Z"
    }
   },
   "outputs": [],
   "source": [
    "#do not run if latex renderer for matplotlib doesn't work\n",
    "from matplotlib import rc\n",
    "rc('text', usetex=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:43:54.541458Z",
     "start_time": "2021-07-31T13:43:54.533739Z"
    }
   },
   "outputs": [],
   "source": [
    "#some settings\n",
    "sizeB = 3.25*1.1\n",
    "sizeH = 2.5*1.1\n",
    "fonts = 0.9\n",
    "color_pr = np.array([0.5647058823529412, 0.9333333333333333, 0.5647058823529412]) - np.array([0.4,0.1,0.4])\n",
    "color_no = '#CC79A7'\n",
    "\n",
    "font ='DejaVu Serif'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:43:54.550007Z",
     "start_time": "2021-07-31T13:43:54.543250Z"
    }
   },
   "outputs": [],
   "source": [
    "## default parameters\n",
    "\n",
    "\n",
    "asym_dist = clapy.dist()\n",
    "\n",
    "dargs = {\n",
    "        'samples': 10000,\n",
    "        'nCells': 100,\n",
    "        'mCells': 100,\n",
    "        'GF': 0.95,\n",
    "        'G1': 0.5,\n",
    "        'S': 0.3,\n",
    "        'G2M': 0.2,\n",
    "        'sCells' : 0.3,\n",
    "        'sSamples' : 0.2\n",
    "}\n",
    "\n",
    "dTC = dargs['G1']+dargs['G2M']+dargs['S']\n",
    "dFS = dargs['S']/dTC\n",
    "X = np.arange(0,dargs['nCells']+1)\n",
    "time_points = np.linspace(0.01,1.965625,22)\n",
    "measure_times = np.ravel(np.array(time_points)[:,np.newaxis]*np.ones(dargs['samples']))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load data for figuers 2 and 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:43:54.858271Z",
     "start_time": "2021-07-31T13:43:54.551925Z"
    }
   },
   "outputs": [
    {
     "ename": "OSError",
     "evalue": "``data.pandas`` does not exist",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mOSError\u001b[0m                                   Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-bec1974f4196>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;31m#load data\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0;32mwith\u001b[0m \u001b[0mpd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mHDFStore\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'data.pandas'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'r'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcomplib\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'zlib'\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcomplevel\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m9\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      3\u001b[0m     \u001b[0mdata_asy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data_asy'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m     \u001b[0mdata_asy_d\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data_asy_d'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0mdata_sym\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mst\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'data_sym'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib64/python3.9/site-packages/pandas/io/pytables.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, path, mode, complevel, complib, fletcher32, **kwargs)\u001b[0m\n\u001b[1;32m    535\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fletcher32\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfletcher32\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    536\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_filters\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 537\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    538\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    539\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m__fspath__\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib64/python3.9/site-packages/pandas/io/pytables.py\u001b[0m in \u001b[0;36mopen\u001b[0;34m(self, mode, **kwargs)\u001b[0m\n\u001b[1;32m    653\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    654\u001b[0m         \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 655\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_handle\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtables\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen_file\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_path\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_mode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    656\u001b[0m         \u001b[0;32mexcept\u001b[0m \u001b[0mIOError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0merr\u001b[0m\u001b[0;34m:\u001b[0m  \u001b[0;31m# pragma: no cover\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    657\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;34m\"can not be written\"\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib64/python3.9/site-packages/tables/file.py\u001b[0m in \u001b[0;36mopen_file\u001b[0;34m(filename, mode, title, root_uep, filters, **kwargs)\u001b[0m\n\u001b[1;32m    318\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    319\u001b[0m     \u001b[0;31m# Finally, create the File instance, and return it\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 320\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mFile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtitle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mroot_uep\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilters\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    321\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    322\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib64/python3.9/site-packages/tables/file.py\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, filename, mode, title, root_uep, filters, **kwargs)\u001b[0m\n\u001b[1;32m    782\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    783\u001b[0m         \u001b[0;31m# Now, it is time to initialize the File extension\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 784\u001b[0;31m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_g_new\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mparams\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    785\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    786\u001b[0m         \u001b[0;31m# Check filters and set PyTables format version for new files.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32mtables/hdf5extension.pyx\u001b[0m in \u001b[0;36mtables.hdf5extension.File._g_new\u001b[0;34m()\u001b[0m\n",
      "\u001b[0;32m/usr/lib64/python3.9/site-packages/tables/utils.py\u001b[0m in \u001b[0;36mcheck_file_access\u001b[0;34m(filename, mode)\u001b[0m\n\u001b[1;32m    155\u001b[0m         \u001b[0;31m# The file should be readable.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    156\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0maccess\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mF_OK\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 157\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mIOError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"``%s`` does not exist\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    158\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mpath\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0misfile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    159\u001b[0m             \u001b[0;32mraise\u001b[0m \u001b[0mIOError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"``%s`` is not a regular file\"\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mOSError\u001b[0m: ``data.pandas`` does not exist"
     ]
    }
   ],
   "source": [
    "#load data\n",
    "with pd.HDFStore('data.pandas', 'r',complib='zlib',complevel=9) as st:\n",
    "    data_asy = st['data_asy']\n",
    "    data_asy_d = st['data_asy_d']\n",
    "    data_sym = st['data_sym']\n",
    "    data_sym_d = st['data_sym_d']\n",
    "    pdfs = st['pdfs']\n",
    "    pdfm = st['pdfm']\n",
    "\n",
    "    pdfs_gamma = st['pdfs_gamma']  \n",
    "    pdfm_gamma = st['pdfm_gamma'] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " # Figure 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "\n",
    "fig = plt.figure(2,figsize=(sizeB*2, sizeH))\n",
    "fig.clear()\n",
    "ax2 = plt.subplot2grid((1,2),(0,1))\n",
    "ax = plt.subplot2grid((1,2),(0,0))\n",
    "\n",
    "ax.set_ylim((0,100))\n",
    "ax.set_xlim((0,2))\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('\\#labeled cells $k$')\n",
    "\n",
    "\n",
    "\n",
    "#calc mean and cdf for asym\n",
    "mean_asym = np.array([np.mean(i[0]) for i in data_asy])\n",
    "mean_asym_d = np.array([np.mean(i[0]) for i in data_asy_d])\n",
    "\n",
    "error_asy =[]\n",
    "for i in data_asy:\n",
    "    x_asy,y_asy = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_asy.append([x_asy[p95a],x_asy[p95e]])\n",
    "error_asy = np.abs(np.array(error_asy).T - mean_asym)\n",
    "\n",
    "error_asy_d =[]\n",
    "for i in data_asy_d:\n",
    "    x_asy_d,y_asy_d = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_asy_d.append([x_asy_d[p95a],x_asy_d[p95e]])\n",
    "error_asy_d = np.abs(np.array(error_asy_d).T - mean_asym_d)\n",
    "\n",
    "\n",
    "#fit Nowakowski\n",
    "lh = Nowakowski_LH(np.array([np.mean(i[0])/dargs['nCells'] for i in data_asy])[:-5:4],time_points[:-5:4])\n",
    "mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=1,s=0.1,\\\n",
    "                      error_r=0.1,error_GF=0.1,error_Tc=0.1,error_s=0.1,\\\n",
    "           limit_Tc=(0.00001,2), limit_r=(0.00001,1),limit_GF=(0,1),limit_s=(0.0001,1),\\\n",
    "           errordef=0.5,print_level=2)\n",
    "mi_old.migrad()\n",
    "res = mi_old.values\n",
    "args = {\n",
    "        'GF': res['GF'],\n",
    "        'G1': res['Tc']*((1-res['r'])*0.5),\n",
    "        'S': res['Tc']*res['r'],\n",
    "        'G2M': res['Tc']*((1-res['r'])*0.5)\n",
    "}\n",
    "\n",
    "\n",
    "ttt = np.linspace(0,2,1000)\n",
    "#plot\n",
    "\n",
    "\n",
    "plt.plot(time_points,np.array(pdfm),label='analytic',color=color_pr,zorder=4)\n",
    "\n",
    "plt.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,**args),':',color=color_no,label=\"fit\",zorder=6)\n",
    "\n",
    "plt.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=1,**dargs),label=\"Nowakowski\",color=color_no,zorder=5)\n",
    "\n",
    "plt.plot(time_points[:-5:4],[np.mean(i[0]) for i in data_asy][:-5:4],'^',color='black',label=\"fit points\",zorder=7 )\n",
    "\n",
    "plt.errorbar(time_points,mean_asym,error_asy,fmt='.',color='black',label=\"simulation\",zorder=2 )\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(list(reversed(handles)), list(reversed(labels)), loc='lower right')\n",
    "\n",
    "\n",
    "\n",
    "axi = fig.add_axes([0,0,1,1])\n",
    "# Manually set the position and relative size of the inset axes within ax1\n",
    "ip = InsetPosition(ax, [0.48,0.45,0.5,0.4])\n",
    "axi.set_axes_locator(ip)\n",
    "axi.set_xticks([]);axi.set_yticks([])\n",
    "\n",
    "axi.errorbar(time_points,mean_asym_d,error_asy_d,fmt='.',color='black',label=\"simulation\",zorder=2 )\n",
    "axi.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=1,**dargs),label=\"Nowakowski\",color=color_no,zorder=5)\n",
    "\n",
    "ax2.set_ylim((0,100))\n",
    "ax2.set_xlim((0,2))\n",
    "ax2.set_xlabel('time')\n",
    "ax2.set_ylabel('\\#labeled cells $k$')\n",
    "\n",
    "#calc mean and cdf for asym\n",
    "mean_sym = np.array([np.mean(i[0]) for i in data_sym])\n",
    "mean_sym_d = np.array([np.mean(i[0]) for i in data_sym_d])\n",
    "\n",
    "error_sym =[]\n",
    "for i in data_sym:\n",
    "    x_sym,y_sym = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_sym.append([x_sym[p95a],x_sym[p95e]])\n",
    "error_sym = np.abs(np.array(error_sym).T - mean_sym)\n",
    "\n",
    "error_sym_d =[]\n",
    "for i in data_sym_d:\n",
    "    x_sym_d,y_sym_d = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_sym_d.append([x_sym_d[p95a],x_sym_d[p95e]])\n",
    "error_sym_d = np.abs(np.array(error_sym_d).T - mean_sym_d)\n",
    "\n",
    "\n",
    "ax2.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=2,**dargs),label=\"Lefevre\",color=color_no,zorder=1)\n",
    "\n",
    "ax2.errorbar(time_points,mean_sym,error_sym,fmt='.-',color='black',label=\"simulation\",zorder=2 )\n",
    "\n",
    "\n",
    "handles, labels = ax2.get_legend_handles_labels()\n",
    "ax2.legend(list(reversed(handles)), list(reversed(labels)), loc='lower right')\n",
    "\n",
    "\n",
    "axi2 = fig.add_axes([0,0,1,1],label=\"ad\")\n",
    "# Manually set the position and relative size of the inset axes within ax1\n",
    "ip2 = InsetPosition(ax2, [0.48,0.22,0.5,0.4])\n",
    "axi2.set_xticks([]);axi2.set_yticks([])\n",
    "axi2.set_axes_locator(ip2)\n",
    "\n",
    "axi2.errorbar(time_points,mean_sym_d,error_sym_d,fmt='.',color='black',label=\"simulation\",zorder=2 )\n",
    "axi2.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=2,**dargs),label=\"Nowakowski\",color=color_no,zorder=5)\n",
    "\n",
    "\n",
    "ml = MultipleLocator(10)\n",
    "ax.yaxis.set_minor_locator(ml)\n",
    "ax2.yaxis.set_minor_locator(ml)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "ax.set_title(\"asymmetric divisoin\")\n",
    "ax2.set_title(\"symmetric divisoin\")\n",
    "bb = ax.get_position()\n",
    "plt.figtext(bb.x0-.06,bb.y1+0.02,'A',size='x-large')\n",
    "\n",
    "bb = ax2.get_position()\n",
    "plt.figtext(bb.x0-.06,bb.y1+0.02,'B',size='x-large')\n",
    "fig.subplots_adjust(left=0.08,bottom=0.16,top=0.9)\n",
    "fig.savefig(\"2.pdf\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Figure 2 gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "\n",
    "fig = plt.figure(2,figsize=(sizeB*2, sizeH))\n",
    "fig.clear()\n",
    "ax2 = plt.subplot2grid((1,2),(0,1))\n",
    "ax = plt.subplot2grid((1,2),(0,0))\n",
    "\n",
    "ax.set_ylim((0,100))\n",
    "ax.set_xlim((0,2))\n",
    "ax.set_xlabel('time')\n",
    "ax.set_ylabel('\\#labeled cells $k$')\n",
    "\n",
    "\n",
    "\n",
    "#calc mean and cdf for asym\n",
    "mean_asym = np.array([np.mean(i[0]) for i in data_asy])\n",
    "mean_asym_d = np.array([np.mean(i[0]) for i in data_asy_d])\n",
    "\n",
    "error_asy =[]\n",
    "for i in data_asy:\n",
    "    x_asy,y_asy = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_asy.append([x_asy[p95a],x_asy[p95e]])\n",
    "error_asy = np.abs(np.array(error_asy).T - mean_asym)\n",
    "\n",
    "error_asy_d =[]\n",
    "for i in data_asy_d:\n",
    "    x_asy_d,y_asy_d = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_asy_d.append([x_asy_d[p95a],x_asy_d[p95e]])\n",
    "error_asy_d = np.abs(np.array(error_asy_d).T - mean_asym_d)\n",
    "\n",
    "\n",
    "#fit Nowakowski\n",
    "lh = Nowakowski_LH(np.array([np.mean(i[0])/dargs['nCells'] for i in data_asy])[:-5:4],time_points[:-5:4])\n",
    "mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=1,s=0.1,\\\n",
    "                      error_r=0.1,error_GF=0.1,error_Tc=0.1,error_s=0.1,\\\n",
    "           limit_Tc=(0.00001,2), limit_r=(0.00001,1),limit_GF=(0,1),limit_s=(0.0001,1),\\\n",
    "           errordef=0.5,print_level=2)\n",
    "mi_old.migrad()\n",
    "res = mi_old.values\n",
    "args = {\n",
    "        'GF': res['GF'],\n",
    "        'G1': res['Tc']*((1-res['r'])*0.5),\n",
    "        'S': res['Tc']*res['r'],\n",
    "        'G2M': res['Tc']*((1-res['r'])*0.5)\n",
    "}\n",
    "\n",
    "\n",
    "ttt = np.linspace(0,2,1000)\n",
    "#plot\n",
    "\n",
    "\n",
    "plt.plot(time_points,np.array(pdfm_gamma),label='analytic',color=color_pr,zorder=4)\n",
    "\n",
    "plt.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,**args),':',color=color_no,label=\"fit\",zorder=6)\n",
    "\n",
    "plt.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=1,**dargs),label=\"Nowakowski\",color=color_no,zorder=5)\n",
    "\n",
    "plt.plot(time_points[:-5:4],[np.mean(i[0]) for i in data_asy][:-5:4],'^',color='black',label=\"fit points\",zorder=7 )\n",
    "\n",
    "plt.errorbar(time_points,mean_asym,error_asy,fmt='.',color='black',label=\"simulation\",zorder=2 )\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "handles, labels = ax.get_legend_handles_labels()\n",
    "ax.legend(list(reversed(handles)), list(reversed(labels)), loc='lower right')\n",
    "\n",
    "\n",
    "\n",
    "axi = fig.add_axes([0,0,1,1])\n",
    "# Manually set the position and relative size of the inset axes within ax1\n",
    "ip = InsetPosition(ax, [0.48,0.45,0.5,0.4])\n",
    "axi.set_axes_locator(ip)\n",
    "axi.set_xticks([]);axi.set_yticks([])\n",
    "\n",
    "axi.errorbar(time_points,mean_asym_d,error_asy_d,fmt='.',color='black',label=\"simulation\",zorder=2 )\n",
    "axi.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=1,**dargs),label=\"Nowakowski\",color=color_no,zorder=5)\n",
    "\n",
    "ax2.set_ylim((0,100))\n",
    "ax2.set_xlim((0,2))\n",
    "ax2.set_xlabel('time')\n",
    "ax2.set_ylabel('\\#labeled cells $k$')\n",
    "\n",
    "#calc mean and cdf for asym\n",
    "mean_sym = np.array([np.mean(i[0]) for i in data_sym])\n",
    "mean_sym_d = np.array([np.mean(i[0]) for i in data_sym_d])\n",
    "\n",
    "error_sym =[]\n",
    "for i in data_sym:\n",
    "    x_sym,y_sym = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_sym.append([x_sym[p95a],x_sym[p95e]])\n",
    "error_sym = np.abs(np.array(error_sym).T - mean_sym)\n",
    "\n",
    "error_sym_d =[]\n",
    "for i in data_sym_d:\n",
    "    x_sym_d,y_sym_d = ecdf(i[0])\n",
    "    p95a = int(np.round(dargs['samples']*0.025))\n",
    "    p95e = int(np.round(dargs['samples']*0.975))\n",
    "    error_sym_d.append([x_sym_d[p95a],x_sym_d[p95e]])\n",
    "error_sym_d = np.abs(np.array(error_sym_d).T - mean_sym_d)\n",
    "\n",
    "\n",
    "ax2.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=2,**dargs),label=\"Lefevre\",color=color_no,zorder=1)\n",
    "\n",
    "ax2.errorbar(time_points,mean_sym,error_sym,fmt='.-',color='black',label=\"simulation\",zorder=2 )\n",
    "\n",
    "\n",
    "handles, labels = ax2.get_legend_handles_labels()\n",
    "ax2.legend(list(reversed(handles)), list(reversed(labels)), loc='lower right')\n",
    "\n",
    "\n",
    "axi2 = fig.add_axes([0,0,1,1],label=\"ad\")\n",
    "# Manually set the position and relative size of the inset axes within ax1\n",
    "ip2 = InsetPosition(ax2, [0.48,0.22,0.5,0.4])\n",
    "axi2.set_xticks([]);axi2.set_yticks([])\n",
    "axi2.set_axes_locator(ip2)\n",
    "\n",
    "axi2.errorbar(time_points,mean_sym_d,error_sym_d,fmt='.',color='black',label=\"simulation\",zorder=2 )\n",
    "axi2.plot(ttt,dargs['nCells']*clapy.cla_det_model(ttt,mode=2,**dargs),label=\"Nowakowski\",color=color_no,zorder=5)\n",
    "\n",
    "\n",
    "ml = MultipleLocator(10)\n",
    "ax.yaxis.set_minor_locator(ml)\n",
    "ax2.yaxis.set_minor_locator(ml)\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "ax.set_title(\"asymmetric divisoin\")\n",
    "ax2.set_title(\"symmetric divisoin\")\n",
    "bb = ax.get_position()\n",
    "plt.figtext(bb.x0-.06,bb.y1+0.02,'A',size='x-large')\n",
    "\n",
    "bb = ax2.get_position()\n",
    "plt.figtext(bb.x0-.06,bb.y1+0.02,'B',size='x-large')\n",
    "fig.subplots_adjust(left=0.08,bottom=0.16,top=0.9)\n",
    "fig.savefig(\"2_gamma.pdf\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "fig = plt.figure(3,figsize=(sizeB, sizeH*1.28))\n",
    "fig.clear()\n",
    "\n",
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "\n",
    "ax1 = plt.subplot2grid((3,1),(0,0),rowspan=2)\n",
    "ax2 = plt.subplot2grid((3,1),(2,0))\n",
    "ax1.plot(-3,-3,\"k-\",label=\"solution\")\n",
    "ax1.plot(-3,-3,\"k.\",label=\"simulation\")\n",
    "\n",
    "ee=10\n",
    "ss=4\n",
    "\n",
    "def cmap(x):\n",
    "    return sb.color_palette(palette='colorblind',n_colors=3)[int((x*1000)%3)]\n",
    "\n",
    "norm = mpl.colors.Normalize(vmin=time_points[:-ee:ss][0], vmax=time_points[:-ee:ss][-1])\n",
    "\n",
    "\n",
    "for i_n,(i,d,p) in enumerate(zip(time_points[:-ee:ss],data_asy[:-ee:ss],pdfs[:-ee:ss])):\n",
    "    h = ax1.hist(d[0],density=True,bins=np.arange(X[0]-0.5,X[-1]+0.5),color=cmap(norm(i)),alpha=0.5,lw=0 )\n",
    "    ax1.plot(X,p,color=cmap(norm(i)),drawstyle='steps-mid')\n",
    "    ax1.bar(X,p,1,linewidth=0.5,color=(0,0,0,0),edgecolor=cmap(norm(i)))\n",
    "    if i_n==2:\n",
    "        ax1.annotate(\"t={:1.2f}\".format(i),xy=(np.mean(d[0]),h[0].max()) ,xytext=(np.mean(d[0])-17,h[0].max()-0.01 ),\\\n",
    "                color=cmap(norm(i)))\n",
    "    else:\n",
    "        ax1.annotate(\"t={:1.2f}\".format(i),xy=(np.mean(d[0]),h[0].max()) ,xytext=(np.mean(d[0])-10,h[0].max()+0.005 ),\\\n",
    "                color=cmap(norm(i)))\n",
    "#ax1.legend()\n",
    "ax1.set_xlim(0,100)\n",
    "ax1.set_ylim(0,0.114)\n",
    "  \n",
    "    \n",
    "ax2.plot(-3,-3,\"k-\",label=\"solution\")\n",
    "ax2.plot(-3,-3,\"k.\",label=\"simulation\")\n",
    "\n",
    "\n",
    "for i,d,p in zip(time_points[:-ee:ss],data_asy[:-ee:ss],pdfs[:-ee:ss]):\n",
    "\n",
    "    ax2.plot(X,np.cumsum(p),\"-\",color=cmap(norm(i)),drawstyle='steps-mid',lw=1 )\n",
    "    xx,y = ecdf(d[0])\n",
    "    ux,cx = np.unique(xx,return_counts=True)\n",
    "    ax2.plot(ux,[y[k:j].max() for k,j in zip(np.cumsum(cx)-cx[0],np.cumsum(cx))],'.',color=cmap(norm(i)) )\n",
    "    #ax2.annotate(\"t={}\".format(i),xy=(np.mean(xx),np.mean(y)),xytext=(np.mean(xx)-10,(2-i)*0.3+0.2 ),\\\n",
    "    #            color=cmap(norm(i)),bbox={'facecolor':'white', 'alpha':1, 'pad':2})\n",
    "\n",
    "#ax2.legend()\n",
    "ax2.set_xlim(0,100)\n",
    "ax2.set_ylim(-0.05,1.11)\n",
    "\n",
    "\n",
    "\n",
    "ax1.set_xticks([])\n",
    "ax1.set_xlabel('')\n",
    "\n",
    "ax2.set_xlabel('\\#labeled cells $k$')\n",
    "\n",
    "ax1.set_ylabel('$P(K = k)$')\n",
    "ax2.set_ylabel('$P(K\\leq k)$')\n",
    "\n",
    "divider = make_axes_locatable(ax1)\n",
    "lax = divider.append_axes(\"top\", size=\"5%\", pad=0.035)\n",
    "lax.set_xticks([])\n",
    "lax.set_yticks([])\n",
    "lax.set_frame_on(False)\n",
    "lax.set_xlim(0,100)\n",
    "a = 20\n",
    "a3 = 37\n",
    "s = 1.2\n",
    "a2 = 2.5*s\n",
    "at= 6\n",
    "\n",
    "lax.text(a+at,0,'analytic')\n",
    "lax.bar(a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.5,color=(0,0,0,0),edgecolor='black')\n",
    "lax.plot(a+s*np.array([-0.50,-0.5,0.5,0.5,1.5,1.5]),[0,2.8,2.8,5,5,0],color='black')\n",
    "lax.plot(a+a2+s*np.array([-0.,1.]),[2.5,2.5],lw=1,color='k')\n",
    "\n",
    "lax.text(a+a3+at,0,'simulation')\n",
    "lax.bar(a+a3+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.0,color=(0,0,0,0),edgecolor=None,alpha=0.5)\n",
    "lax.plot(a+a3+a2+s*np.array([0.5]),[2.5],'k.')\n",
    "\n",
    "\n",
    "bb = ax1.get_position()\n",
    "atxt = plt.figtext(bb.x0-0.09,bb.y1+0.04 ,'A',size='x-large',zorder=7)\n",
    "\n",
    "bb = ax2.get_position()\n",
    "atxt = plt.figtext(bb.x0-0.09,bb.y1+0.06 ,'B',size='x-large',zorder=7)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace=0)\n",
    "fig.savefig(\"3.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 3 gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "fig = plt.figure(3,figsize=(sizeB, sizeH*1.28))\n",
    "fig.clear()\n",
    "\n",
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "\n",
    "ax1 = plt.subplot2grid((3,1),(0,0),rowspan=2)\n",
    "ax2 = plt.subplot2grid((3,1),(2,0))\n",
    "ax1.plot(-3,-3,\"k-\",label=\"solution\")\n",
    "ax1.plot(-3,-3,\"k.\",label=\"simulation\")\n",
    "\n",
    "ee=10\n",
    "ss=4\n",
    "\n",
    "def cmap(x):\n",
    "    return sb.color_palette(palette='colorblind',n_colors=3)[int((x*1000)%3)]\n",
    "\n",
    "norm = mpl.colors.Normalize(vmin=time_points[:-ee:ss][0], vmax=time_points[:-ee:ss][-1])\n",
    "\n",
    "\n",
    "for i_n,(i,d,p) in enumerate(zip(time_points[:-ee:ss],data_asy[:-ee:ss],pdfs_gamma[:-ee:ss])):\n",
    "    h = ax1.hist(d[0],density=True,bins=np.arange(X[0]-0.5,X[-1]+0.5),color=cmap(norm(i)),alpha=0.5,lw=0 )\n",
    "    ax1.plot(X,p,color=cmap(norm(i)),drawstyle='steps-mid')\n",
    "    ax1.bar(X,p,1,linewidth=0.5,color=(0,0,0,0),edgecolor=cmap(norm(i)))\n",
    "    if i_n==2:\n",
    "        ax1.annotate(\"t={:1.2f}\".format(i),xy=(np.mean(d[0]),h[0].max()) ,xytext=(np.mean(d[0])-17,h[0].max()-0.01 ),\\\n",
    "                color=cmap(norm(i)))\n",
    "    else:\n",
    "        ax1.annotate(\"t={:1.2f}\".format(i),xy=(np.mean(d[0]),h[0].max()) ,xytext=(np.mean(d[0])-10,h[0].max()+0.005 ),\\\n",
    "                color=cmap(norm(i)))\n",
    "#ax1.legend()\n",
    "ax1.set_xlim(0,100)\n",
    "ax1.set_ylim(0,0.114)\n",
    "  \n",
    "    \n",
    "ax2.plot(-3,-3,\"k-\",label=\"solution\")\n",
    "ax2.plot(-3,-3,\"k.\",label=\"simulation\")\n",
    "\n",
    "\n",
    "for i,d,p in zip(time_points[:-ee:ss],data_asy[:-ee:ss],pdfs_gamma[:-ee:ss]):\n",
    "\n",
    "    ax2.plot(X,np.cumsum(p),\"-\",color=cmap(norm(i)),drawstyle='steps-mid',lw=1 )\n",
    "    xx,y = ecdf(d[0])\n",
    "    ux,cx = np.unique(xx,return_counts=True)\n",
    "    ax2.plot(ux,[y[k:j].max() for k,j in zip(np.cumsum(cx)-cx[0],np.cumsum(cx))],'.',color=cmap(norm(i)) )\n",
    "    #ax2.annotate(\"t={}\".format(i),xy=(np.mean(xx),np.mean(y)),xytext=(np.mean(xx)-10,(2-i)*0.3+0.2 ),\\\n",
    "    #            color=cmap(norm(i)),bbox={'facecolor':'white', 'alpha':1, 'pad':2})\n",
    "\n",
    "#ax2.legend()\n",
    "ax2.set_xlim(0,100)\n",
    "ax2.set_ylim(-0.05,1.11)\n",
    "\n",
    "\n",
    "\n",
    "ax1.set_xticks([])\n",
    "ax1.set_xlabel('')\n",
    "\n",
    "ax2.set_xlabel('\\#labeled cells $k$')\n",
    "\n",
    "ax1.set_ylabel('$P(K = k)$')\n",
    "ax2.set_ylabel('$P(K\\leq k)$')\n",
    "\n",
    "divider = make_axes_locatable(ax1)\n",
    "lax = divider.append_axes(\"top\", size=\"5%\", pad=0.035)\n",
    "lax.set_xticks([])\n",
    "lax.set_yticks([])\n",
    "lax.set_frame_on(False)\n",
    "lax.set_xlim(0,100)\n",
    "a = 20\n",
    "a3 = 37\n",
    "s = 1.2\n",
    "a2 = 2.5*s\n",
    "at= 6\n",
    "\n",
    "lax.text(a+at,0,'analytic')\n",
    "lax.bar(a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.5,color=(0,0,0,0),edgecolor='black')\n",
    "lax.plot(a+s*np.array([-0.50,-0.5,0.5,0.5,1.5,1.5]),[0,2.8,2.8,5,5,0],color='black')\n",
    "lax.plot(a+a2+s*np.array([-0.,1.]),[2.5,2.5],lw=1,color='k')\n",
    "\n",
    "lax.text(a+a3+at,0,'simulation')\n",
    "lax.bar(a+a3+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.0,color=(0,0,0,0),edgecolor=None,alpha=0.5)\n",
    "lax.plot(a+a3+a2+s*np.array([0.5]),[2.5],'k.')\n",
    "\n",
    "\n",
    "bb = ax1.get_position()\n",
    "atxt = plt.figtext(bb.x0-0.09,bb.y1+0.04 ,'A',size='x-large',zorder=7)\n",
    "\n",
    "bb = ax2.get_position()\n",
    "atxt = plt.figtext(bb.x0-0.09,bb.y1+0.06 ,'B',size='x-large',zorder=7)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace=0)\n",
    "fig.savefig(\"3_gamma.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 4 (careful takesa round 3 min for 10e6 samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T13:44:05.941404Z",
     "start_time": "2021-07-31T13:44:01.784221Z"
    }
   },
   "outputs": [],
   "source": [
    "m,s=clapy.log_params(1,0.3)\n",
    "mc,sc=clapy.log_params(1,clapy.calc_sigma_true(0.3,dargs['G1'],dargs['G2M']))\n",
    "SIZE = 10000000\n",
    "tchc  = np.random.lognormal(mean=mc,sigma=sc,size=SIZE)\n",
    "tc  = np.random.lognormal(mean=m,sigma=s,size=SIZE)\n",
    "tcd = np.random.lognormal(mean=m,sigma=s,size=SIZE)\n",
    "age = np.random.uniform(size=SIZE)\n",
    "w =np.where(age < dargs['G1']+dargs['S'],dargs['G1']*tc-age*tc,dargs['G1']*tcd+tc-age*tc )\n",
    "whc=(dargs['G1']+dargs['G2M']-age)*tchc\n",
    "x,y = ecdf(w)\n",
    "xhc,yhc = ecdf(whc)\n",
    "cdfwt = pd.Series(y,index=pd.to_datetime(x+10, unit='s'))\n",
    "cdfwtrs = cdfwt.resample('10ms').mean()\n",
    "cdfwrs =pd.Series(cdfwtrs.values ,index=np.array(cdfwtrs.index.microsecond/1000000) + np.array(cdfwtrs.index.second) - 10)\n",
    "\n",
    "\n",
    "cdfwhct = pd.Series(yhc,index=pd.to_datetime(xhc+10, unit='s'))\n",
    "cdfwhctrs = cdfwhct.resample('10ms').mean()\n",
    "cdfwhcrs =pd.Series(cdfwhctrs.values ,index=np.array(cdfwhctrs.index.microsecond/1000000) + np.array(cdfwhctrs.index.second) - 10)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-07-31T14:13:31.037718Z",
     "start_time": "2021-07-31T14:13:29.168996Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAD5CAYAAADspDPqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAArWUlEQVR4nO3deXxb5Z3v8c85kmzJi7zGcRwri1NwYjlkISXEgbK1jaGlC+1EDLdl2AydmQbuzCSXmemtuc2k08W57ZRwp21caKcFihjSaTtMEB2gbJEhJEDiyImbjUS2Y8erJO+Szrl/GAs7sWXZkX20PO/Xyy9b58jn+R3bX5/tOc+RVFVVEQQhachaFyAIwtwSoReEJCNCLwhJRoReEJKMCL0gJBkRekFIMnqtCwintLRU6xIEIS40NjZG/N6YDj1Mb2UEIRlNd+M4rd37Bx98cFoLFwQh9kQcerfbTV1dHW63ezbrEQRhlkUcepfLxf33388zzzwzm/UIgjDLprV7X1lZyYsvvjhbtQhREFRUjp7u4ujpLs62eunxDREMKpO+t9s3SGCS+UJiiuhEnsvlwmq1YrFYxr0WtDc4HMB1qpPm9l7auwd481ALPb5BAALBj+6lMqboyMk0snShGb1O5ri7h7bOPhQVTKk6rCX55GUZ0etkjCk6MtNSKFmYxccs2WSYDEiSpNUqClEWUeidTicWiwWXy4XFYmHv3r0i9DHA/lIje/67gbKUc1iywGRK5e41C1m1uhxTRgZ+RaZ3CHy9Awy3naa38zxNPh99UjrrrlvMwnwTmcEe2s+dp7npBCnnW8gYbMOry6JJKeCplw00D6UTNKSxcH4WN6wp4vrV80lPkZGNaUiS6OYRj6RIbq11OBxUVlYCIyf07r77bl566aVZL660tFRcspvE6RYP//mvj3FL5lF0OhlDdgFq0E+wtwdlsPei90spJvTmPIJ9HpQB35gZMro0M7IxjZR5i0gpWIy/u5WhlhP4u9tACQCgSDpkNRj6Nn12AVlXfRbjIiu69Cz0GTmzvs7CxKabkym39BfuylssFsxms9jF15CiKLz79I/ZlHaEwlu/TvrlVyHpDR/NHx5E9Q+hBv2ogWGQZPTZBaEts+IfIuDtQJJ16LPmIcm6CdtRlSABXydKfy9qYBhVNvAfb37Aq++f4y/LgP3P0/mHJwAwlawi76a/IKVg8ayvv3Bpwu6fOZ1OHnroIZxOZ2ja6CW7nTt34nK5Zrc6YULvP7+HFYPvk33b35NRtnFc4AHkFOPI1tecjyG3CENO4bhdcdmQSkrewpHpkwQeQJJ1GLIKSF1QgtGyHNPCZdxhu4nbPncN2+vMtF7/CEse/jXFD/wI2ZhB08+20v3Gs6hKcNJlCtqLaPdeK2L3fmL7dz5Ed+YyNj2gXWep379xkl8838D2+zdQviwfgP5T73P+dz8idf5SCjf/w0X/jITZMd2ciDMxcSbg6yZ3sAnziqs1reNz1y7jC9ct43u/OkCXd+RqQVrJaorv+7/4O5vpfvM5TesTJidCH2ea39tHt5JB6aqVWpfC/6hcwaL5mXz/VwdCfQH0mbnkf+Yv6an7D4ZaT2lcoTCRqIfe4XBQU1Mz4Ty73Y7T6Rx3jkCYHk/D25zWl5CbZdK6FHSyxLavrONcRx+7f1vP6JFiWslqMldeT/vz/yqO72NQ1ENfWVmJz+e7aLrdbqe8vJyKigocDke0m00KyvAAps5GgsWrtC4lJDszlep71/PKATe/e/1kaHruTXcS8HbgO/SKhtUJE5mz3XuXy4XZbAYQN+3MUP/peoZVPUVla7QuZZxlxdk8fOfH+cXzDRw81gaAzpRBzidsdL/2a5Shfo0rFMaKmWP6Xbt2UVpaOu5DGK/rdCMfBHIpv6xA61Iusm7FfO7YtJwfPP0unZ4BAMxrP41syqTH+RuNqxPGmrPQW61WvF4vAJmZmRfN37JlC42NjeM+hPE8TafxGeaRFwPH8xP58o2XUbIwi5onDxJUVCRZR+4NX8HzzgsEBy7uJShoY1ZO5Lnd7lDHHbvdDoDNZuPIkSM4nU5uv/32aDebHDwtKFlFWlcxKVmW+Ns71tJ8vpc9rxwHIO2yK9Gbc/G99weNqxNGRX24rMrKylA/fRgJ+0RfC9OjBgOYBjtJXVisdSlh5WQa2WJbzXd+sZ+1pQV8zJJN1lW30v3Gv5O1/lYkneiwo7WYOaYXwvN3nUNGwVy8VOtSpnRVWSGfvGoxO586SN+An4yV16EqAXpd+7QuTUCEPm4MtrvpUdIoKpqndSkRufdzVlJTdOx86iCqLgXzmk/jfVcMwBILROjjhMd9inPBLBbOy9C6lIgYU/T877vXc6KphydfOErGyusYav4T/u5WrUtLeiL0caLv3Bk6ySUrI0XrUiI2L8fE39/5cX7z6glO+YykLlhGb4PYxdeaCH2cULqbGc6YH3fDVllL8rilYgmPPvsephUb6T3yOjF8Y2dSEKGPA2owQEp/O7rc2D5zP5k7byljcDjIyx0L8Hc0M3z+jNYlJTUR+jgweuY+c8ESrUuZEVOqnr+87Qp+/WYbctEK+hre1LqkpCZCHwf83a30qiYKi/K1LmXGPl5WyOrL5/GWbwF9jfu1LiepidDHgYHOVjqC6XFz5n4yVZ9fyYvNufg7mxnubNG6nKQlQh8HvG3n6AqmsyA/XetSLsmC/HSu2VhOu1xA/5/E1l4rIvRxYLCrjQFDFsaUmH/I8JRuvbaEA31FdNSLgVS0IkIfBxRfB0FjrtZlREVBThq6JWuQ2k8S6O3RupykJEIfB/QDXUjm+D2Jd6HrP7mBrmAG5w+/pXUpSUmEPsYpw4OkBPtJzYm9gTNmavmSPM6lLqHp/be1LiUpidDHuICnHYCM/EKNK4mu7MvXYOr+E4oYOHPOidDHuICnnX41lbx5ifWsuJXXXotJHeRP7x/RupSkE/XTwXa7PfRI64qKinHzxg59feE8YWJD3W10BtMpyo7NIbJmKic/jxMphbQfqGP52tgZ3TcZRHVLH26Ya5fLhdfrpaKiQjwDbxp621vpUjKYl52mdSlRl7poJXJrA/6AonUpSSWqoQ83zLXVamX37t08+OCD44bTEsIb6GylR82Iq1tqI7Vk3dUskVs56GrSupSkMmfH9G63m/vvv5+VK1dO+AQcMQT2xIKedvzGnLi7pTYSGUvK0MkSx/a/o3UpSSWqoQ83zLXD4aCyspKqqqrQMf9YYgjsicn9XajpeVqXMStkfQrKvMvwNx1hcCigdTlJI6qhn2iY69EhsCsrK3E4HLhcLjZu3BjNZhOWEhjG4PdhyI6PcfFmYl7ZlVyub+VtlxhGa65E/ez9hcNcj762WCwTbuGFyQW9nQCk5yXWNfqx0ktWUfzqUzz37kmuWxufg4TEG3GdPoYFfJ340ZOdlxj97ieSUrgEUtLoPVVPb/+w1uUkBRH6GBbwdeFR0piXm3iX60ZJkkx6yRVYTW0cPHZe63KSggh9DBvoOk930MS8BOuYc6G0pVdQltrGOw1tWpeSFEToY1hfZzs9Sjr5CR5609IryPB3cvzYKYJB0VFntonQx7Chng765QxMqfE/eEY4+uz5yJl5LFRbOHamW+tyEp4IfQxTersImrK1LmPWSZJE2qIyPp7r4Z0GcelutonQxzB5oBs5PbHurpuM0bKCxXIr+8Vx/awToY9RqhLE4O8lJStxRswJx7hoBcb+NrrOt3O+u1/rchKaCH2MCvb2IKGSlpu4vfHGMuQXI5syWZvj4fDxdq3LSWgi9DEq4OtEQSJrXnKEXpJkjJblrM32cOh4h9blJDQR+hgV8HbiVUzk58T3WPfTYbSUUay2cOh4u3jI5SwSoY9R/V3n6VbSEv4a/VjGRWWkeJvo9/XibvNpXU7CiqsLwIqi0NLSgt/v17qUWTdoKiT7xj9n0Hee072zd4xrMBgoKipClrX//59auBTJkMLV83s5dLyDRYVmrUtKSHEV+paWFsxmc2h0nkQ2kGmkf1ghr2jhrLbj9XppaWmhuFj7O9wkWYexeDlrvT04j7dz67UlWpeUkLT/9z4Nfr8/KQIPgBJAlXSz3ozZbI6pPSfjIisLlWaOnOwgqIjj+tkQV6FPKkoQdHG1IxYVpsVlGHrOEBga5IMWj9blJKSoh95ut+N0OscNd33hvNHRdOKB1+u9aGTfSLhcLqqrq2fUhqqqSGoQSU6+0KcuWIYk69hQ2E/9yU6ty0lIczoEdlZWFhUVFdx8883RbHZWmc3mGY3ea7VaJ53n9Xqpra2dvA1FQUJF1idf6CWdAWNxKauzujlyUlyvnw1zNgT23r17cbvdk+4FTFfvgJ9Oz8CMP3oHJj6OdTqduFyuUCjHbrFdLhd33303TqeTmpqa0F7Lyy+/HHrPbbfddtEyvV4vdrud2tpa3G53aBzB0fH/x7Yx8t5neOvAe7z62h8vajMZGBdZKQ624DrViSKO66Mu7Kbk8ccf595776WpqSl0dre3t5eMjIwZNWa1WqmoqODuu+++aOu5a9cuHnvssYiWEwwq3LvjD/QPznwE1TSjnqe334xON/7/ntvtxuFwhMb2G7vFHv16dE+mvLycrKwsjhz56NFM5eXlF7VlNpux2Ww4nU4cDgdVVVXs27cvtLyxbezcuZNH/vFhhrra2PbPP+SLn981rs1kYFxchmHfHoYHP86ZVi9Li7K0LimhhA19ZmYm1dXVNDc3U1xcjKqquFwu9uzZM+H7ww2BvXLlyrCFbNmyhS1btoybNtnY9zqdzOP/+9MMDc889Kkp+osCDyPhqqiooKamhh07dlx0tWDs4J4TXT7s6em5aJmjW/Wxj/Ly+Xy4XK4JDwN6urtIQWZ0qPtkG1DUWHQ5kiRxdUE/9Sc7ROijLGzoN2/ezObNm2loaKCsrAyAhoaGSd9vs9mw2+14PJ5xQ2DbbDYqKyupra3F6/VeNGLuTGSYDGSYDJe8nAs5HA6sViu33HILZrMZt9sdOlRxuVy43W68Xi9utzv02uVykZmZidPpxOfz4XQ6sVgs4w5xRl+PXZ7HM3J2euy0qqoqHI4/MC8vj8/ccsuEbYY7X5AIJL0Bo2UFa3u7ePtkJ5+7dpnWJSUUSZ1mJ+e6ujo2bNgwW/WMU1paOu6hF6dPn2bp0qVz0raWBjpa6R8OznrHnFGx+HPtcf6G1kN1bG+5gSe/dTOynHhP+ImWC3MylbBb+nvuuQer1UpDQ0No9765uXnOQp+sVCUASXi5bizj4nJSXnuGwEA/H5zzUrJQ7OJHS9i/rCeeeAIYv3UPt3svRIekBEA2al2GplIXLEPSp1BR2MfhE+0i9FEU0SU7l8tFU1MTTU1N1NXVzXZNSS3UMScJe+ONJck6TIvKWGPuFPfXR1lEf1k2m42amhokSeK+++6b7ZqSm5q8HXMuZFxSzoIDf8TV3EkwqEx4tUWYvoh+imMfk5yTkxwDNc6Fibr4qsGRy5A6feRXJmbaVTjWmRaVo/e4kf0DnGjq0bqchBFR6H/yk59w++23Y7PZ+OlPfzrbNSWNibr4qsEAQVVGr5/8Drspu/EmiJT5i5FTTFyzoJ/DJ8QufrREFPqNGzdSVlZGWVlZqINJb2/vrBamlQu7zDqdznHdYL1e76TTHnzwwVBX2dHljPbCG/0egAcffBCv1zthF9839+3jhz99gnf2v4Xdbg9do59eN96P2o3nbrySrMNoWcEqcyeHxXF91EQU+traWh555BGqq6upra2lurqahx56aLZrCys42EfA1zXjj+Bg34TLHe0ya7VacTgcVFRUYLFYqKioYOPGjdjt9kmnAWzbtg0Y6U5rs9moqKhg7969VFRUYLPZqK6u5tFHH8VsNk/YxXf9lWvx9vaz6oqVlJeX43Q6J6yprKxs0m68Y9sd23XY54u/IaiMi60U+ps4eqYLf0A88ioaIjpbVFVVddG1eS3P4qtKkLOPfQ11aObjo0upaSz5218gyeN3oyfqMhupC7vLer3ecd10LRYLPT09F00fO18NBFCRyMrKCnVpnm433onajVfGRVZ0rzyJPjjMcXc3ZUvztC4p7kUU+ok642jZQUeSdSz6+k9Q/UMzX4Yh9aLAj7qwy+zobr7L5aKqqgrgomkul4uGhoZQ4KqqqnjhhRewWCzccsstOBwO6uvr2bFjB3fddRc/+tGPQsuBMV18PT20nDs3rouvxWKZVjfese3GezfekXHzUrmmqJ/6Ex0i9NGgxrDLL7983OtTp05pUsc3v/nNiKZFQ9+5D9TOtvZZWfZktPq5Rqrl6X9SX//Zv6jf+PGbWpcSky7MyVTEhc8pjG4pp5oWDaqqIquBpO+YcyHT0isoHDzN0dNd+ANBrcuJeyL0U7Barfz85z+fclpUKMGRjjmG6N89GM/SLluHrseNWeqjUTzK+pKJ0McQNRhAVaVpdcxJBil5RRhyi7hxfhf14nr9JYur0KuqGlPDNUebEvATRMYQpmNOtPn9/rh4hFTaZesoN7h570/i4ZaXKq4OHgsKCmhpaUFREvN67WBPB4O9PsyDKrI0N/ePy7JMQUHBnLR1KdIuW4f5nb2c7myjt3+YjLQUrUuKW1EPvd1uD10Xnuhat8PhwGKxzOiyUUZGxozH54sHR578LadOn+dL3/yO1qXEHKNlOboUI6sy2jl0vIONq4q0LiluzdkQ2DDSaaS+vj50fVkYL+A5TzBNXIeeiCTrMC1bzYacDt5tPK91OXFtzobABjhy5MiUA2QmM7m/E505OZ5HPxNpJWuwBM7wbmNbXJyHiFVzdiJvqq6tu3btorS0dNxHMlFVFeNwD6a8Qq1LiVmmklXoB7uRvW00nU/MG77mQlRDH24IbCD0oIvRu8PG2rJlC42NjeM+kkmwrwc9QbIXiGPVyegzckgpWMI187o4eKxN63LiVlRDb7PZQrd9jh0CG0b+IYw+CGKiseGT3WBnK0FVoqBIhD4cU8kqVpraOHBUhH6mon72/sIx7ce+NpvNbN++PdpNJoTO5ia6lXTKCpLkUdwzlFaympx3XqCx6zz9g37SjKIj03TFVeecROZtO4dXypyVB3gkEqNlBbIssdzUyfuio86MiNDHiKGuVoZSxPiDU5H0BkxLr+ATeefFLv4MidDHCMXXgZohrtFHIn35BpYMH+fg0VbxVNsZEKGPEakD7ehzxEm8SKRftg6dv4+cwSZONvdoXU7cEaGPAYp/iIygl+yFi7UuJS7IxnTSSlZzY34bdfXntC4n7ojQxwDfOTeSBIUlJVqXEjfSl19NKadwHmoRvfOmSYQ+BrSdPkmXks7CBeKYPlJpl1+Fwd+Loec0Z9vib5RfLYnQxwBP8xm8ulwMevHriJTOlEHasjV8Kq8Z52Gxiz8d4q8sBvi7WvCnx/497bEmY+V1LFdP8vbh6I9XmMhE6GOAvreNlDxx5n660j52JXpZxdRxlJZ2cQNOpEToNaaqKuZAF+YiceZ+umRDKpkrKrgxu4nX3mvWupy4IUKvsZ62VlKkAIVLl2ldSlzKWHkdS5XTvH3whDiLHyEReo2dO3mCIVXP/EXFWpcSl4yLVqA357Oo9zAnm8SITJEQoddYT/MZPHIOep34VcyEJMlkf7ySGzJO8tq7Z7UuJy6IvzSNDbU3MZQmhsi6FJlX3EAmvTQd2k8wmJgjJUeTCL3GDB43KfOXal1GXNOZMkkv28ha1cU74s67KUU99Ha7PTQs1lherxen04ndbp9wpNxk1Nc3QL7SzrzL4ucpsrEq96pbKNO7efON97QuJebN2RDYTqcTi8WCzWZj9+7d0Ww2bp1yNaAnyOKycq1LiXupC5YhFZWx4NxrtHSIa/bhzNkQ2JWVlaGHYBQXizPVAO3HG/Do8jCkpWtdSkJYcNPtrE89xSuvvq91KTFtzh9rVVtby44dOy6avmvXLh577LG5LkdT/taTDGcv0rqMhGFaVIY/dymyy0H/Z9eL8fMmMadDYDscDiorKyd8wk2yDYGtKCoZ/c2kWy7XupSEYvn0HazTN4qtfRhzNgS20+lk9+7d1NTUUFNTE81m41LTuU4KpG4WrhBP/Imm9JJVDGWXMPzOHvyBoNblxCRJjeG+i6WlpQm7xX/9xVcpPPD/uOzvn0bSid3QaPKdbaTtl9+g6eq/4/pPbtC6nFk33ZyI6/Qa6Th+hF5joQj8LMhcVEpPXjnDb9vF1n4CIvQa8AeCpHc1krpY7NrPluW33ccStZnXnhd9Qi4kQq+BQw1NLJVbWbxuo9alJKyM+cX0Lr2ejPo9eLz9WpcTU0ToNfDBwbcI6lJIX7Rc61IS2hW33U267Mdp/5XWpcQUEfo5pigqUtMh/POtSLJO63ISmt6URsrGr1DS9hLHDoruuaNE6OfYCXcXyzjLgjUVWpeSFEqv30Rb7lr6HY8y2CtGzQUR+jm3/439ZMiD5C5fp3UpSWPdXzzEoJrC4ce/jaqIs/ki9HOof9APJ/YRmFeKzpShdTlJIy09jdzbtqH3NnP4Vz9M+mG1ROjn0B/fPsGV+hMs/MTntC4l6Swv+xi+ir9C7z7Aid//XOtyNCVCP0dUVaXJ+Qek1DQyLhe79lq47qYNHF32VYL1L3L6hSe1LkczIvRzpO5wC+X+w2RduUmctdfQl26/hXcW2hg++HtOP/coin9I65LmnAj9HOgf9PPH3z3PfL2PeVdt0rqcpCZJEnfc+QX2We7i/NH3OP7jv2O4o0nrsuaUCP0cePo/3+VW3T6yr/kz9BnZWpeT9HQ6mbvvvJmOa7ZxuNPEmdpteA79MWlO8InQz7LX32sio34PaTl55G38otblCB+SJInP3FBG2Z0P81/D62l9/id88LOHGfigXuvSZp0I/Sx668g53t/zS9aZzmD50kNIujkfqEiYgnVZPl/9nw/w2/wq3mwy0Pz0P+H+5SMMuo8m7JZf3E8/CwJBhWf/+xhdb+6hMv0IRbZ/IK1ktdZlCVPY72rlt//1Fqv667jScArMheReeSOZy9djyI3dB4xONydRD73dbg8NgFlRURHxvInEW+iDikpdfQsvv+jkE8OvsijFw4IvbyOtZJXWpQkRUhSVA0fbeKOuAcOZt1ib8gFFcieB9HlkLltF5tJyUhdehj67AEmKjR3l6eYkqvubo0NgW61WqqurxwU73Lx4paoqPb1DnG720Hi4ge7G91ihHueruvOkrtjI/E13o8/I0bpMYRpkWeIqayFXWQvpH7yWg8fOs/e9YwydOcziAydYdvgtsuglKBsYTitAzinCNG8h6fnzSc+bjyF7HrrMPGRDqtarMqmoht7lcoXCfOEQ2OHmzdTw0BCu118lGPB/NDHMjsuFc1RVBVUdma6OvFZUBVVRQVVQFQVFCRL0BwgG/ASHhwkODxIcGoDhfnTDPjLVPvJ1PtZLAYbMC8i1biR37U0YchdEZR0F7aQZDVy7eiHXrl6IotzImVYvJ9w9HHK3EOhwo/OeI+V0B9kn95Mt95Gr6yNVCgAQQMewbCQopxLUpaLoUgnqUlDlFNAZPvzQgawHnWHkfI+kA1kGSUaSdUg6PaXXb8KcmxvV9QoberfbTU1NDU1NTdx///2haQBVVVVRLWQmQ2B3tLQQ2P8sMuFuopAm+CrM+yRQkQAJVZJRkVCl0V+OHvQpSCYT+txsjFklmPMLyC1eRGrBYnTpWdOqX4gfsiyxtCiLpUVZsH4x8NHYe/2Dfjy9w3R7B2jt7Gagu51Afy/BwT6UwX5U/wAEBpGDfmRlGMnvB2UAWQkiqQFkdfSzgowCqoqMgoTKuUWXz23oLRYLt99+e2jo6lG33XYbVqv1ol30cENgTzU89pYtW9iyZcu4aaWlpWGLL1q6lKJvPB72PYIw29KMBtKMBhbkp0NJPnCZ1iWFNaMzEcXFxRPuoocbAnuieYIgzL1pH9M7HA6ys7Ox2WwTzr9w+tjXk32PIAhzJ6LQu91uXC5X6Mk0o5fdBEGIPxHt3lssltAxfGVlZegEnyAI8WdGx/QVFRXU1dVFuxZBEObAjK7T19fXs2HDxI8Lcjgc1NfXs23btovmTbdHHkx9Bl8QhOmZ8jr9M888Q1NTEw6HIzTN5/Oxffv2Cb+nsrISp9N50fSZ9siLp264gqCF6W4Yp7xO/+ijj15SQaNmo0eeEPuCikpXdy/tLa14O9rwezoZ7vcRHOhDGepHDgyiCw5CMABKAEkNIqlBdGoQWQ2iI4g00kUq1NtSQmW0C9VIV0o+fA+M6Vo1zUqnfr9OlkIdPrMyUwkGFXx9wwwa81mY7scb0OPz9bMgY+RN6cuvJvvqz3Huqf9D3qfuIe1jaydcrhIYpvmJh8m66jOYV39ymnVPX8zc6zmTHnlCbFEUlRN/Ok1z/bsMtZ5G5ztHTrCDbKkPowQGZPqlNIZlE4reiKI3oeqNKKlpSDoDks6ArNePPNRTp0OV9SDrkSRppGuqJI37mtDrDz8DSKP/IqQPO1lKo5MjEP5NprQU/AEFFZX8QjNBRcHX7KWgsIi89CBGvwxdfeTlyEgSGPIt6DJyyNt0H8aFl0/eqs5A/qZ7Scmfm6ticxb62eiRJ2gvEAxypO4t2t99jWzPcfJkL7lSBv0ZxRiWXE5m0SbMCy1kzy/EkJkdM3emRcvqso++TgPmTfCeqW6rliQJ0+LyaJYVVtRD73A4Qtf1rVYrdrsdm82GzWbDbrfj8XhEj7wE0NvdxaH/fJa0s07SGSA9vQT9ms+St+5qls4v+mjLK8QcMYiGMC3+Ph/vPvdzstxv4pHMBEtv4opPVpKRLW4h1oqm99MLiUtVVY6/9Fv8+/8dRU2ja81dXLVpE3q9GM473ojQC1Ma6mqj4akfkOI5w+kFm7j+9jvISI/dQSKE8ETohbA6jx7g/H/8gPPBeRTcup3Prp78LLQQH0TohQmpqsrZl55leP9zvJeynsqv/xW55jStyxKiQIReuIiqKhx/7l8JNr7JO/O/zO13fZkUgzh2TxQi9MI4qhLkT0/WMHSmnuPL7+WrX/6kuPyWYETohRBVCXLsyRoGzrg4u/qv+NKtG7UuSZgFIvQCMLJLf+zX/8LQmXra1n2dL9w88V2UQvxLrD6Rwoyoqsrx535M4NRBWtZ8jZtF4BOaCL3AyT88S7DxdY6X3cMtn71W63KEWSZCn+TOvPUK6jv/TsOiP+OLt92kdTnCHBDH9Ems9dhhBl/6CQ3zNvHFr3xJnKVPEiL0Scp3/hwde77PqfS13HrfvciyCHyyEKFPQsODAxx9YjteeT43fe1v0OvEUV4yEb/tJKMoCm/t/i5SYIjV9/4j6SZx40yyEaFPMq899QvyvY3M+9L/Ir8gT+tyBA2I0CeRd155leIze1GvuZclK5ZrXY6gERH6JHHm+ClSnbvpWnIT1us/pXU5gobEibwk0OvtpeXZ76GkL2LDHfdrXY6gMbGlT3DBoML+2u9ikAKsrfoGsixukU12Ykuf4N54+ucUDpwg5893YMq4eOhxIfmE3dI//vjjADQ1NYWm9fb2zm5FQtQcev01Fp55geDGe1mw7GNalyPEiLBb+szMTKqrq2lubqa4uBhVVXG5XOzZs2eu6hNmqOX0B8iv/5Q2y01cc4M4cSd8JGzoN2/ezObNm2loaKCsbORRHg0NDXNSmDBz/X19uH/9HYZMi7jmq+LEnTBeRCfyRgMP4PF4wr7XbrfjdDonfHLt6IMrRx9vJUSfoii8s/u7yGqAtff9gzhxJ1wk7Jb+nnvuwWq10tDQENq9b25unvTZ9OEeR+31eqmurqasrIwHHnggumshhNQ98wvm9R0n48+2k5GVpXU5QgwKG/onnngCgLq6ulDQw+3eT/U46q1bt2KxWDCbzTMuWJhcw743KDi1F+/6+ykrFePTCxOLaPfe5XLR1NREU1MTdXV1M2pobNBra2svmr9r1y5KS0vHfQiRa3efJfjHH9O04Aau/NSntS5HiGERhd5ms7F7925qa2vZtGnTpO8L9zhqp9MZ2spPdLy/ZcsWGhsbx30IkRkaGODUk/9MZ2ox1971Na3LEWJcRKF3Op2Ul5fzrW99C5/PN+n7bDYbR44cwel0hh5HbbfbASgvL8ftduN0OrHZbFEoXYCRQS33/+z7yMowq+/9R/Q6ceJOCC+iHnmSJLF582YAenp6wr73wkCPvjabzVit1hmUKISzf89T5HsaSPnCI2TnZmtdjhAHItrSnz17lrq6utBxvRAbjh/YT86x3+Fd8xWWlZdN/Q2CQIShv++++3C5XOzcuZPy8vLZrkmIQHdrK/2OXbjzN7D+M7dqXY4QRyK+4ea+++6bzTqEaQgMD3Hs33YwbMjn2nu2aF2OEGfEXXZxqO6JH5Lu93FZ1U5SUgxalyPEGXE/fZw5+PweCjreJe3mvyF//jytyxHikAh9HPmg/hAZ79tpL9vM8itXa12OEKfE7n2c8HV20PP7H9CRvYabvvglrcsR4pgIfRwIBgLUP/FtgrKZa+77W/H4KeGSiNDHgbf/bRdpQx0U3vU9jEbxcArh0ohj+hh3+A//Rf45J7qb/prC4iKtyxESgAh9DGtuPEbK/l/SsuzzrNxwtdblCAlC7N7HqH5PD23PfZ+u9DJutN2hdTlCAhGhj0GKEuT9x/8ZVUplQ9VWZFnskAnRI0Ifg/Y/+VMy+pvJuePbpGeka12OkGDEJiTGNLz+MnlnX2F44/1YSpZoXY6QgEToY0jb6ZPwei1niytZd8N1WpcjJCixex8jhvp6cf/6O3iMy7juq/doXY6QwEToY4Cqqhz42XeQVZmrqh5GrxM7YMLsEaGPAe888zhZvlOkfXk7mVlieHBhdolNisYaX3+J7JMv0nfVvZQsv0zrcoQkIEKvoaajLtTXazldfDPrP/1JrcsRkoTYvddIZ0szHb/5Hu2ZV/CpO+/WuhwhiYjQa6CzpYUzP/8mHsN8rvvaVmRZ3CorzB0R+jl29tgxOvZ8l15DHuv/ejspqeJWWWFuidDPkUBQ4e3fPkfO0d/gyVzBxgceJsVo1LosIQlFPfR2ux2LxQIw7lHVU81LVH0DQxx65WUCh/aSr3TQW/Z5rv/i7WL0G0EzUQ19uOfTh5sXT1RVRVFUVFRQFFRgeDhAb28/fZ4eBjwees6fZ+D8WZSOMxQOfkCuFKCz4EoWfuGbZBXM13oVhCQX1dCHez79VM+unw2dngG27XqDwaEgq+Tj3GJ4GwkVUBndzo68Hvk8Mm3854/mj0yb6pybCYBU9IZ8lOwi9OvuZNn6jSxPy4jqugnCTMXMMf2uXbt47LHHorrMnEwjf/Pnaxn2B5H8ZahdK0GSkJBQpZEHc6qSBNKH3RXkD2P+4XvQjbxGGvmXIMmgIIM8Ml+SRubp9TLpaSbSc3JITc9E0sXMj1UQLhLVv85wz6cPNw9Gnk+/Zcv4RzSVlpZeUj2yLLFyWf6YKYsvaXmCkAii2iMv3PPpJ5onCMLck1RVVbUuYjKXuqUXhGTR2NgY8XtjOvQzVVpaOq0fQqy3I9oSbUWzHXHDjSAkGRF6QUgyIvSCkGQSMvRf//rXE6od0ZZoK5rtJOSJPEEQJpeQW3pBECYn+ovGoLm8UzHc8txuNxaLBa/Xi9l86QN2OhwO6uvr2bZt27TqiHZb0Vwvr9fLkSNHcLvdZGVlUVlZOW5+NNdrqrYiXa+E2tI7HA5qamomnGe323E6nTidzqi0FW55ozcUjXY7nu5yy8vLqaiowOFwRDxvJsItz+v1Ul1dPenPcyYqKyvx+XzTqiPabUV7vZxOJxaLBZvNxu7du8fNi/Z6hWtrOuuVUKGfqz+q2QyLy+UK/Zee6E7FyeZFuy2ArVu38sADD0RlK38pdURbNNersrIytCUvLi4eNy/a6xWuLYh8vRIq9JNJ1LDMprG119bWalhJdM3WetXW1rJjx46oLW+6bU1nvZIi9HPpUv+oLuVOxWi2NboraTabo3ZINJM6om021svhcFBZWYnH4xk3fTbWa7K2prNecXcib/SuvVETndC40Ex/+JO1NVVYysvLQz/8qqqqiNuDkbsR7XY7Ho9n3J2KNpttwnmXIlxb5eXluN1uPB4PNpvtktuCkT9Yt9uNy+XCarXO2nqFayva6+V0Otm9e3dod/vRRx+dtfUK19Z01iuhrtM7HA7sdjtbt24d94uG2T/rPdqW1+sN/fC9Xu+U/5AEYa4lVOgFQZiaOKYXhCQjQi8ISUaEXhCSjAh9HHM6nVRXV1/09UTzJ/v+qToRjV2G1+uNWo+5aC5LmB4R+jg29irERFckprpKEclVjLHvMZvNM74a4fV6x/VbuJRlTbQ8IXJxd50+2dTW1oYuDVqtVhwOBxaLBYvFgtVqjXg5brc71GmjoqIitMyGhgacTiculwubzYbH45m0DZfLhd1uZ/v27TidTmpra6mqqmLfvn1s27Yt1MboDR9j2xkdCbmiogKr1TrpshwOx7h/NJOt84XLEyIntvQxrLa2FqvVGuqBVVNTg9VqxWKxsHfv3mkta/RGjfLycp555plx0ysqKrDZbOzcuTNsG2PDNRrM0c+j/RNgZCteXl4eCjwQej26jMmWlZmZicViobKykr17905az4XLEyIntvQxrL6+PrQLbLPZ2LdvXygg0/1jH92aR7JLHWkbF/ZGzMrKwmKxkJWVNel9B6M95MIt68LvDVfPZMsTJidCH8MeeOCBUJ/qrKwstm3bRm1t7bjdX7fbHdrKjn49GhqXyzVuvsViweVy4fP5xoXF5XLhdDrZunUrHo9nyjZGv6epqSk03el0YrVaqampwWKxsHHjxovOB4wua3T3frJluVwuvF4vTU1NfPvb3x5Xz2jNFy5PiJzokSdETW1tLTabDbPZTG1t7bhbQYXYIbb0QtRYLBaOHDlCVlZW6LUQe8SWXhCSjDh7LwhJRoReEJKMCL0gJBkRekFIMiL0gpBkROgFIcn8f+llwj5JrwKwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 257.4x253.44 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "fig = plt.figure(4,figsize=(sizeB, sizeH*1.28))\n",
    "fig.clear()\n",
    "\n",
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "\n",
    "ax1 = plt.subplot2grid((3,1),(0,0),rowspan=2)\n",
    "ax2 = plt.subplot2grid((3,1),(2,0))\n",
    "\n",
    "\n",
    "bins = np.linspace(-1,1.5,100)\n",
    "mbins = bins[:-1] + bins[1] - bins[0]\n",
    "app,_ =np.histogram(whc,bins=bins,density=True)\n",
    "true,_ =np.histogram(w,bins=bins,density=True)\n",
    "ax1.plot(mbins,true,label=r'simulation')\n",
    "ax1.plot(mbins,app,label=r'approximation')\n",
    "ax1.legend(loc='center left')\n",
    "\n",
    "\n",
    "  \n",
    "ax2.plot(cdfwrs,label=r'simulation')\n",
    "ax2.plot(cdfwhcrs,label=r'approximation')\n",
    "#ax2.legend()\n",
    "\n",
    "\n",
    "ax1.set_xticks([])\n",
    "ax1.set_xlabel('')\n",
    "\n",
    "ax2.set_xlabel('cell labeling time t')\n",
    "\n",
    "ax1.set_ylabel('pdf')\n",
    "ax2.set_ylabel('cdf')\n",
    "\n",
    "\n",
    "bb = ax1.get_position()\n",
    "atxt = plt.figtext(bb.x0-0.09,bb.y1+0.04 ,'A',size='x-large',zorder=7)\n",
    "\n",
    "bb = ax2.get_position()\n",
    "atxt = plt.figtext(bb.x0-0.09,bb.y1+0.06 ,'B',size='x-large',zorder=7)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "fig.subplots_adjust(hspace=0)\n",
    "fig.savefig(\"4.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Figure 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.transforms as plttransforms\n",
    "import matplotlib.lines as pltlines\n",
    "\n",
    "ho_abc = 0.35\n",
    "ho_abc2 = 0.45\n",
    "alpha_no = 0.4\n",
    "alpha_pr = 0.3\n",
    "\n",
    "with pd.HDFStore(\"data.pandas\", 'r') as store:\n",
    "    print(store.keys())\n",
    "    data_sv = store['minuit_full_new_s100']\n",
    "    data_sv.index.rename(['model','sPopulation','sSample','N'],inplace=True)\n",
    "    data_sv['M'] = 5\n",
    "    data_sv.set_index('M', append=True, inplace=True)\n",
    "    data_sv['GF'] = 0.95\n",
    "    data_sv.set_index('GF', append=True, inplace=True)\n",
    "    data_sv = data_sv.reorder_levels(['model','sPopulation','sSample','M','GF','N'])\n",
    "    data_sv.rename({'sigma_cell' : 'sc','sigma_cell_error':'sc_error','sigma_sample':'ss','sigma_sample_error':'ss_error'}, axis=\"columns\",inplace=True)\n",
    "\n",
    "    data_s = store['minuit_full_new_s5']\n",
    "    data_s.index.rename(['model','sPopulation','sSample','N'],inplace=True)\n",
    "    data_s['M'] = 5\n",
    "    data_s.set_index('M', append=True, inplace=True)\n",
    "    data_s['GF'] = 0.95\n",
    "    data_s.set_index('GF', append=True, inplace=True)\n",
    "    data_s = data_s.reorder_levels(['model','sPopulation','sSample','M','GF','N'])\n",
    "    data_s.rename({'sigma_cell' : 'sc','sigma_cell_error':'sc_error','sigma_sample':'ss','sigma_sample_error':'ss_error'}, axis=\"columns\",inplace=True)\n",
    "\n",
    "#set dataframe\n",
    "f = plt.figure(4,figsize=(sizeB*2, sizeH*2.5))\n",
    "f.clear()\n",
    "\n",
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "fh = f.get_figheight() \n",
    "fb = f.get_figwidth() \n",
    "\n",
    "f.subplots_adjust(left=0.06,right=0.99,hspace=0.28,top=0.94)\n",
    "f.clear()\n",
    "p = {'Tc' : 1,'r' : 0.33 , 'ss':0.2,'sc' :0.3,'GFf':0.95,'m' : 5 }\n",
    "\n",
    "A = 'lognorm' \n",
    "B = 'old' \n",
    "#plot swarm wit N_data points\n",
    "N_data = 100\n",
    "#plot and calc\n",
    "data,sdata = calc_statistic(data_sv,p)\n",
    "axs = (3,10)\n",
    "\n",
    "ax_1 = plt.subplot2grid(axs,(0,0),colspan=2)\n",
    "ax_2 = plt.subplot2grid(axs,(0,2),colspan=2)\n",
    "ax_3 = plt.subplot2grid(axs,(0,4),colspan=2)\n",
    "ax_4 = plt.subplot2grid(axs,(0,6),colspan=2)\n",
    "ax_5 = plt.subplot2grid(axs,(0,8),colspan=2)\n",
    "\n",
    "ax_6 = plt.subplot2grid(axs,(1,0),colspan=2)\n",
    "ax_7 = plt.subplot2grid(axs,(1,2),colspan=2)\n",
    "ax_8 = plt.subplot2grid(axs,(1,4),colspan=2)\n",
    "ax_9 = plt.subplot2grid(axs,(1,6),colspan=2)\n",
    "ax_10 = plt.subplot2grid(axs,(1,8),colspan=2)\n",
    "        \n",
    "ax_11 = plt.subplot2grid(axs,(2,1),colspan=4)\n",
    "ax_12 = plt.subplot2grid(axs,(2,6),colspan=4)\n",
    "\n",
    "ax_10.set_xticks([0,0.3,1])\n",
    "\n",
    "types = ['Tc','r','GFf','ss','sc']\n",
    "axes_0 = [ax_1,ax_2,ax_3,ax_4,ax_5]\n",
    "axes_1 = [ax_6,ax_7,ax_8,ax_9,ax_10]\n",
    "names = [r'$\\tau$',r'$f$',r'$g$',r'$\\sigma$',r'$\\hat{\\sigma}_c$']\n",
    "\n",
    "\n",
    "hist = dict()\n",
    "mhist = dict()\n",
    "scalef = []\n",
    "for typ,ax_s in zip(types,axes_0):\n",
    "        try:\n",
    "            hist[typ+\"o\"] = np.histogram(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "        except:\n",
    "            pass\n",
    "        hist[typ+\"f\"] = np.histogram(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "    \n",
    "for typ,ax_s in zip(types,axes_0):\n",
    "        try:\n",
    "            bins = np.linspace(hist[typ+'o'][1][0],hist[typ+'o'][1][-1],len(hist[typ+'o'][1])*2)\n",
    "            hist[typ+\"o\"] = ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                      bins=bins,alpha=alpha_no ,density=True,color=color_no,ec=color_no,linewidth=0.5)\n",
    "            ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                      bins=bins,density=True,histtype='step',color=color_no,linewidth=1.0)\n",
    "            mhist[typ+\"o\"] = hist[typ+\"o\"][0].max()\n",
    "        except:\n",
    "            mhist[typ+\"o\"] = 0\n",
    "            pass\n",
    "        bins = np.linspace(hist[typ+'f'][1][0],hist[typ+'f'][1][-1],len(hist[typ+'f'][1])*2)\n",
    "        ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                  bins=bins,density=True,histtype='step',color=color_pr,linewidth=1.0)\n",
    "        hist[typ+\"f\"] = ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                  bins=bins,alpha=alpha_pr ,density=True,color=color_pr,edgecolor=color_pr,linewidth=0.5)\n",
    "\n",
    "        mhist[typ+\"f\"] = hist[typ+\"f\"][0].max()\n",
    "        \n",
    "        scalef.append(mhist[typ+\"f\"] < mhist[typ+\"o\"])\n",
    "        \n",
    "scalef = np.max( scalef ) \n",
    "for typ,ax_s,name in zip(types,axes_0,names):\n",
    "    if scalef:\n",
    "        a = 'o'\n",
    "        b = 'f'\n",
    "    else:\n",
    "        a = 'f'\n",
    "        b = 'o'\n",
    "    if mhist[typ+a] == 0:\n",
    "        ax_s.set_ylim(0,mhist[typ+b]*1/0.9)\n",
    "        ymax = mhist[typ+b]*1/0.9\n",
    "    else:\n",
    "        ax_s.set_ylim(0,mhist[typ+a]*1/0.9)\n",
    "        ymax = mhist[typ+a]*1/0.9\n",
    "    try:\n",
    "        xo = np.mean(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "        ind = np.sum(hist[typ+'o'][1]<xo)-1\n",
    "        ax_s.plot([xo,xo],[0, ymax*0.249 ],\"-\",color=color_no,lw=1.2,zorder=1.1)      \n",
    "        ax_s.plot([xo,xo],[0, ymax*0.25],\"-\",lw=2,color='.15',zorder=1)\n",
    "    except:\n",
    "        pass\n",
    "    aco = sdata[typ+'_p'].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    acn = sdata[typ+'_p'].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    ax_s.set_title(r'{}'.format(name)  )\n",
    "    xf = np.mean(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "    ind = np.sum(hist[typ+'f'][1]<xf)-1\n",
    "    ax_s.plot([xf,xf],[0, ymax*0.249 ],\"-\",color=color_pr,lw=1.2,zorder=1.1)      \n",
    "    ax_s.plot([xf,xf],[0, ymax*0.25 ],\"-\",color='.15',lw=2,zorder=1)        \n",
    "    #ax_s.plot([p[typ],p[typ]],[0,ymax*0.2],\"-\",color='gray',lw=4,zorder=8)\n",
    "    ax_s.set_yticks(np.linspace(0,1,5)*ymax)\n",
    "    ax_s.set_yticklabels([])\n",
    "    \n",
    "    \n",
    "\n",
    " \n",
    "\n",
    "data,sdata = calc_statistic(data_s,p)\n",
    "\n",
    "\n",
    "\n",
    "hist = dict()\n",
    "mhist = dict()\n",
    "scalef = []\n",
    "for typ,ax_s in zip(types,axes_1):\n",
    "        try:\n",
    "            hist[typ+\"o\"] = np.histogram(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "        except:\n",
    "            pass\n",
    "        hist[typ+\"f\"] = np.histogram(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "    \n",
    "\n",
    "    \n",
    "for typ,ax_s in zip(types,axes_1):\n",
    "        try:\n",
    "            bins = np.linspace(hist[typ+'o'][1][0],hist[typ+'o'][1][-1],len(hist[typ+'o'][1])*2)\n",
    "            hist[typ+\"o\"] = ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                      bins=bins,alpha=alpha_no ,density=True,color=color_no,ec=color_no,linewidth=0.5)\n",
    "            ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                      bins=bins,density=True,histtype='step',color=color_no,linewidth=1.0)\n",
    "            mhist[typ+\"o\"] = hist[typ+\"o\"][0].max()\n",
    "        except:\n",
    "            mhist[typ+\"o\"] = 0\n",
    "            pass\n",
    "        bins = np.linspace(hist[typ+'f'][1][0],hist[typ+'f'][1][-1],len(hist[typ+'f'][1])*2)\n",
    "        ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                  bins=bins,density=True,histtype='step',color=color_pr,linewidth=1.0)\n",
    "        hist[typ+\"f\"] = ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                  bins=bins,alpha=alpha_pr ,density=True,color=color_pr,ec=color_pr,lw=0.5)\n",
    "\n",
    "        mhist[typ+\"f\"] = hist[typ+\"f\"][0].max()\n",
    "        \n",
    "        scalef.append([mhist[typ+\"f\"] < mhist[typ+\"o\"]])\n",
    "        \n",
    "scalef = np.max( scalef )\n",
    "for typ,ax_s,name in zip(types,axes_1,names):\n",
    "    if scalef:\n",
    "        a = 'o'\n",
    "        b = 'f'\n",
    "    else:\n",
    "        a = 'f'\n",
    "        b = 'o'\n",
    "    if mhist[typ+a] == 0:\n",
    "        ax_s.set_ylim(0,mhist[typ+b]*1/0.9)\n",
    "        ymax = mhist[typ+b]*1/0.9\n",
    "    else:\n",
    "        ax_s.set_ylim(0,mhist[typ+a]*1/0.9)\n",
    "        ymax = mhist[typ+a]*1/0.9\n",
    "    try:\n",
    "        xo = np.mean(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "        ind = np.sum(hist[typ+'o'][1]<xo)-1\n",
    "        ax_s.plot([xo,xo],[0, ymax*0.249 ],\"-\",color=color_no,lw=1.2,zorder=1.1)      \n",
    "        ax_s.plot([xo,xo],[0, ymax*0.25],\"-\",lw=2,color='.15',zorder=1)\n",
    "    except:\n",
    "        pass\n",
    "    aco = sdata[typ+'_p'].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    acn = sdata[typ+'_p'].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    #ax_s.set_title(r'{}'.format(name)  )\n",
    "    xf = np.mean(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "    ind = np.sum(hist[typ+'f'][1]<xf)-1\n",
    "    ax_s.plot([xf,xf],[0, ymax*0.249 ],\"-\",color=color_pr,lw=1.2,zorder=1.1)      \n",
    "    ax_s.plot([xf,xf],[0, ymax*0.25 ],\"-\",color='.15',lw=2,zorder=1)      \n",
    "    #ax_s.plot([p[typ],p[typ]],[0,ymax*0.2],\"-\",color='gray',lw=4,zorder=8)\n",
    "    ax_s.set_yticks(np.linspace(0,1,5)*ymax)\n",
    "    ax_s.set_yticklabels([])\n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "with pd.HDFStore('data.pandas','r',complevel=9) as store:\n",
    "     allimnew = store['s100_n10'].loc[idx[:,:0.2:,:],:].copy()\n",
    "\n",
    "allimnew.reset_index(inplace=True)\n",
    "allimnew.rename({\"level_0\":\"model\"}, axis=\"columns\",inplace=True)\n",
    "allimnew.set_index(['model','start','leng','N'],inplace=True)\n",
    "allimnew.sort_index(inplace=True)\n",
    "\n",
    "        \n",
    "im = plot2d_pandas(allimnew.loc[idx[B]],'Tc','start','leng',aspect='auto',vmin=0.75,vmax=1.25,\\\n",
    "        cmap=plt.cm.bwr,axis=ax_12)\n",
    "\n",
    "cbar = plt.colorbar(im,ax=ax_12)\n",
    "tl = cbar.ax.get_yticklabels()\n",
    "plt.draw()\n",
    "tl = [float(i.get_text().replace('$','')) for i in tl]\n",
    "\n",
    "cbar.ax.set_yticklabels(['{:2.0f}'.format( (i-1)*100 )+\"\\%\" for i in tl])\n",
    "\n",
    "im = plot2d_pandas(allimnew.loc[idx[A]],'Tc','start','leng',aspect='auto',vmin=0.75,vmax=1.25,\\\n",
    "        cmap=plt.cm.bwr,axis=ax_11)\n",
    "#plt.colorbar(im,ax=ax_11,extend='neither')\n",
    "\n",
    "\n",
    "ax_11.set_ylabel('measurement duration')\n",
    "ax_12.set_ylabel('measurement duration')\n",
    "ax_11.set_xlabel('time of first measurement')\n",
    "ax_12.set_xlabel('time of first measurement')\n",
    "\n",
    "#add lines\n",
    "for azoom,aall in zip(axes_0,axes_1):\n",
    "    azoom.set_zorder(10)\n",
    "    aall.set_zorder(6)\n",
    "    zoom_effect02(azoom,aall,color='black',ls=':',zorder=0)\n",
    "\n",
    "    \n",
    "#add text\n",
    "ax_1.set_ylabel('probability density')\n",
    "bb = ax_1.get_position()\n",
    "atxt = plt.figtext(bb.x0 - ho_abc/fb,bb.y1 + 0.05/fh,'A - 100 samples',size='x-large',zorder=7)\n",
    "\n",
    "ax_6.set_ylabel('probability density')\n",
    "bb = ax_6.get_position()\n",
    "btxt = plt.figtext(bb.x0 - ho_abc/fb,bb.y1 + 0.05/fh,'B - 5 samples',size='x-large',zorder=7)\n",
    "btxt.set_bbox(dict(facecolor='white', alpha=1, edgecolor=None,pad=0))\n",
    "\n",
    "\n",
    "\n",
    "#ax_tmp = plt.subplot2grid(axs,(2,0))\n",
    "#ax_tmp.set_axis_off()\n",
    "#bb = ax_tmp.get_position()\n",
    "#plt.figtext(bb.x0 - ho_abc/fb,bb.y1 + 0.05/fh,'C',size='x-large')\n",
    "\n",
    "bb = ax_11.get_position()\n",
    "plt.figtext(bb.x0 - 0.7/fb,bb.y1 + 0.05/fh,'C',size='x-large')\n",
    "bb = ax_12.get_position()\n",
    "plt.figtext(bb.x0 - 0.7/fb,bb.y1 + 0.05/fh,'D',size='x-large')\n",
    "\n",
    "plt.draw()\n",
    "        \n",
    "tmp = ax_4.get_xticklabels()\n",
    "tmp[-2].set_horizontalalignment('right')\n",
    "ax_4.set_xticklabels(tmp)\n",
    "\n",
    "tmp = ax_5.get_xticklabels()\n",
    "tmp[1].set_horizontalalignment('left')\n",
    "ax_5.set_xticklabels(tmp)\n",
    "\n",
    "tmp = ax_8.get_xticklabels()\n",
    "tmp[-2].set_horizontalalignment('right')\n",
    "ax_8.set_xticklabels(tmp)\n",
    "\n",
    "for ax_s in axes_0:\n",
    "    for t in ax_s.get_xticklabels():\n",
    "        #try:\n",
    "        #    bbox = t.get_window_extent().transformed(ax_s.transData.inverted())\n",
    "        #except:\n",
    "        #    continue\n",
    "        #rect = mpl.patches.Rectangle((bbox.x0,bbox.y0), bbox.x1 - bbox.x0 ,bbox.y1 - bbox.y0 , ec=\"none\", color='white',zorder=0)\n",
    "        #rect.set_clip_on(False)\n",
    "        #ax_s.add_patch(rect)\n",
    "        t.set_bbox(dict(facecolor='white', alpha=1, edgecolor=None,pad=0))\n",
    "\n",
    "\n",
    "for typ,ax_s in zip(types,axes_0):\n",
    "    labels = []\n",
    "    trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "    for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "        la = l.get_text()\n",
    "        if np.abs(t - p[typ])<0.005 and typ!='sc':\n",
    "            lold = la.replace('$','')\n",
    "            #lfont = l.get_font_properties()\n",
    "            #lfont.set_weight('bold')\n",
    "            #l.set_fontproperties(font)\n",
    "            la = '$\\\\mathbf{'+lold+'}$'\n",
    "            line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "            ax_s.add_line(line)\n",
    "        labels.append(la)\n",
    "    ax_s.set_xticklabels(labels)  \n",
    "    \n",
    "for typ,ax_s in zip(types,axes_1):\n",
    "    labels = []\n",
    "    trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "    for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "        la = l.get_text()\n",
    "        if np.abs(t - p[typ])<0.005 and typ!='sc':\n",
    "            lold = la.replace('$','')\n",
    "            #lfont = l.get_font_properties()\n",
    "            #lfont.set_weight('bold')\n",
    "            #l.set_fontproperties(font)\n",
    "            la = '$\\\\mathbf{'+lold+'}$'\n",
    "            line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "            ax_s.add_line(line)\n",
    "        labels.append(la)\n",
    "    ax_s.set_xticklabels(labels)\n",
    "\n",
    "ax_s = axes_0[-1]\n",
    "trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "found=False\n",
    "for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "    la = l.get_text()\n",
    "    if np.abs(t - clapy.calc_sigma_true(p[typ],0.5,0.2))<0.005:\n",
    "        lold = la.replace('$','')\n",
    "        #lfont = l.get_font_properties()\n",
    "        #lfont.set_weight('bold')\n",
    "        #l.set_fontproperties(font)\n",
    "        #a = '$\\\\mathbf{'+lold+'}$'\n",
    "        a = '$\\\\math{'+lold+'}$'\n",
    "        line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "        ax_s.add_line(line)\n",
    "        found=True\n",
    "    labels.append(la)\n",
    "nt = clapy.calc_sigma_true(p[typ],0.5,0.2)\n",
    "if not found:\n",
    "    line = pltlines.Line2D([nt,nt],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "    ax_s.add_line(line)\n",
    "\n",
    "ax_s = axes_1[-1]\n",
    "trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "found=False\n",
    "for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "    la = l.get_text()\n",
    "    if np.abs(t - clapy.calc_sigma_true(p[typ],0.5,0.2))<0.005:\n",
    "        lold = la.replace('$','')\n",
    "        #lfont = l.get_font_properties()\n",
    "        #lfont.set_weight('bold')\n",
    "        #l.set_fontproperties(font)\n",
    "        #a = '$\\\\mathbf{'+lold+'}$'\n",
    "        a = '$\\\\math{'+lold+'}$'\n",
    "        line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "        ax_s.add_line(line)\n",
    "        found=True\n",
    "    labels.append(la)\n",
    "nt = clapy.calc_sigma_true(p[typ],0.5,0.2)\n",
    "if not found:\n",
    "    line = pltlines.Line2D([nt,nt],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "    ax_s.add_line(line)    \n",
    "    \n",
    "    \n",
    "bb = ax_2.get_position()\n",
    "lax = plt.axes([bb.x0 - 0.2/fb,bb.y1 + 0.2/fh ,0.4,0.014])\n",
    "lax.set_xticks([])\n",
    "lax.set_yticks([])\n",
    "lax.set_frame_on(False)\n",
    "lax.set_xlim(0,100)\n",
    "a = 20\n",
    "a3 = 60\n",
    "s = 1.2\n",
    "at= 4\n",
    "\n",
    "lax.text(a+at,0,'stochastic model',size='large')\n",
    "lax.bar(a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.5,color='w',edgecolor=color_pr,alpha=alpha_pr)\n",
    "lax.plot(a+s*np.array([-0.50,-0.5,0.5,0.5,1.5,1.5]),[0,2.8,2.8,5,5,0],color=color_pr)\n",
    "lax.bar(a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.0,color=color_pr,edgecolor=None,alpha=alpha_pr)\n",
    "\n",
    "lax.text(a3+a+at,0,'Nowakowski model',size='large')\n",
    "lax.bar(a3+a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.5,color='w',edgecolor=color_no,alpha=alpha_no)\n",
    "lax.plot(a3+a+s*np.array([-0.50,-0.5,0.5,0.5,1.5,1.5]),[0,2.8,2.8,5,5,0],color=color_no)\n",
    "lax.bar(a3+a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.0,color=color_no,edgecolor=None,alpha=alpha_no)\n",
    "\n",
    "bb = ax_12.get_position()\n",
    "ax_11.set_position([ax_11.get_position().x0,bb.y0,bb.width,bb.height])\n",
    "\n",
    "ax_11.set_title('Bias stochastic model')\n",
    "ax_12.set_title('Bias Nowakowski model')\n",
    "\n",
    "\n",
    "plt.draw()\n",
    "f.savefig(\"5.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# figure5 gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.transforms as plttransforms\n",
    "import matplotlib.lines as pltlines\n",
    "\n",
    "ho_abc = 0.35\n",
    "ho_abc2 = 0.45\n",
    "alpha_no = 0.4\n",
    "alpha_pr = 0.3\n",
    "\n",
    "with pd.HDFStore(\"data.pandas\", 'r') as store:\n",
    "    print(store.keys())\n",
    "    data_sv = store['minuit_full_new_s100']\n",
    "    data_sv.index.rename(['model','sPopulation','sSample','N'],inplace=True)\n",
    "    data_sv['M'] = 5\n",
    "    data_sv.set_index('M', append=True, inplace=True)\n",
    "    data_sv['GF'] = 0.95\n",
    "    data_sv.set_index('GF', append=True, inplace=True)\n",
    "    data_sv = data_sv.reorder_levels(['model','sPopulation','sSample','M','GF','N'])\n",
    "    data_sv.rename({'sigma_cell' : 'sc','sigma_cell_error':'sc_error','sigma_sample':'ss','sigma_sample_error':'ss_error'}, axis=\"columns\",inplace=True)\n",
    "\n",
    "    data_s = store['minuit_full_new_s5']\n",
    "    data_s.index.rename(['model','sPopulation','sSample','N'],inplace=True)\n",
    "    data_s['M'] = 5\n",
    "    data_s.set_index('M', append=True, inplace=True)\n",
    "    data_s['GF'] = 0.95\n",
    "    data_s.set_index('GF', append=True, inplace=True)\n",
    "    data_s = data_s.reorder_levels(['model','sPopulation','sSample','M','GF','N'])\n",
    "    data_s.rename({'sigma_cell' : 'sc','sigma_cell_error':'sc_error','sigma_sample':'ss','sigma_sample_error':'ss_error'}, axis=\"columns\",inplace=True)\n",
    "\n",
    "#set dataframe\n",
    "f = plt.figure(4,figsize=(sizeB*2, sizeH*2.5))\n",
    "f.clear()\n",
    "\n",
    "sb.set(context=\"paper\",style='ticks',font=font,font_scale=fonts)\n",
    "fh = f.get_figheight() \n",
    "fb = f.get_figwidth() \n",
    "\n",
    "f.subplots_adjust(left=0.06,right=0.99,hspace=0.28,top=0.94)\n",
    "f.clear()\n",
    "p = {'Tc' : 1,'r' : 0.33 , 'ss':0.2,'sc' :0.3,'GFf':0.95,'m' : 5 }\n",
    "\n",
    "A = 'gamma' \n",
    "B = 'old' \n",
    "#plot swarm wit N_data points\n",
    "N_data = 100\n",
    "#plot and calc\n",
    "data,sdata = calc_statistic(data_sv,p)\n",
    "axs = (3,10)\n",
    "\n",
    "ax_1 = plt.subplot2grid(axs,(0,0),colspan=2)\n",
    "ax_2 = plt.subplot2grid(axs,(0,2),colspan=2)\n",
    "ax_3 = plt.subplot2grid(axs,(0,4),colspan=2)\n",
    "ax_4 = plt.subplot2grid(axs,(0,6),colspan=2)\n",
    "ax_5 = plt.subplot2grid(axs,(0,8),colspan=2)\n",
    "\n",
    "ax_6 = plt.subplot2grid(axs,(1,0),colspan=2)\n",
    "ax_7 = plt.subplot2grid(axs,(1,2),colspan=2)\n",
    "ax_8 = plt.subplot2grid(axs,(1,4),colspan=2)\n",
    "ax_9 = plt.subplot2grid(axs,(1,6),colspan=2)\n",
    "ax_10 = plt.subplot2grid(axs,(1,8),colspan=2)\n",
    "        \n",
    "ax_11 = plt.subplot2grid(axs,(2,1),colspan=4)\n",
    "ax_12 = plt.subplot2grid(axs,(2,6),colspan=4)\n",
    "\n",
    "ax_10.set_xticks([0,0.3,1])\n",
    "\n",
    "types = ['Tc','r','GFf','ss','sc']\n",
    "axes_0 = [ax_1,ax_2,ax_3,ax_4,ax_5]\n",
    "axes_1 = [ax_6,ax_7,ax_8,ax_9,ax_10]\n",
    "names = [r'$\\tau$',r'$f$',r'$g$',r'$\\sigma$',r'$\\hat{\\sigma}_c$']\n",
    "\n",
    "\n",
    "hist = dict()\n",
    "mhist = dict()\n",
    "scalef = []\n",
    "for typ,ax_s in zip(types,axes_0):\n",
    "        try:\n",
    "            hist[typ+\"o\"] = np.histogram(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "        except:\n",
    "            pass\n",
    "        hist[typ+\"f\"] = np.histogram(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "    \n",
    "for typ,ax_s in zip(types,axes_0):\n",
    "        try:\n",
    "            bins = np.linspace(hist[typ+'o'][1][0],hist[typ+'o'][1][-1],len(hist[typ+'o'][1])*2)\n",
    "            hist[typ+\"o\"] = ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                      bins=bins,alpha=alpha_no ,density=True,color=color_no,ec=color_no,linewidth=0.5)\n",
    "            ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                      bins=bins,density=True,histtype='step',color=color_no,linewidth=1.0)\n",
    "            mhist[typ+\"o\"] = hist[typ+\"o\"][0].max()\n",
    "        except:\n",
    "            mhist[typ+\"o\"] = 0\n",
    "            pass\n",
    "        bins = np.linspace(hist[typ+'f'][1][0],hist[typ+'f'][1][-1],len(hist[typ+'f'][1])*2)\n",
    "        ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                  bins=bins,density=True,histtype='step',color=color_pr,linewidth=1.0)\n",
    "        hist[typ+\"f\"] = ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                  bins=bins,alpha=alpha_pr ,density=True,color=color_pr,edgecolor=color_pr,linewidth=0.5)\n",
    "\n",
    "        mhist[typ+\"f\"] = hist[typ+\"f\"][0].max()\n",
    "        \n",
    "        scalef.append(mhist[typ+\"f\"] < mhist[typ+\"o\"])\n",
    "        \n",
    "scalef = np.max( scalef ) \n",
    "for typ,ax_s,name in zip(types,axes_0,names):\n",
    "    if scalef:\n",
    "        a = 'o'\n",
    "        b = 'f'\n",
    "    else:\n",
    "        a = 'f'\n",
    "        b = 'o'\n",
    "    if mhist[typ+a] == 0:\n",
    "        ax_s.set_ylim(0,mhist[typ+b]*1/0.9)\n",
    "        ymax = mhist[typ+b]*1/0.9\n",
    "    else:\n",
    "        ax_s.set_ylim(0,mhist[typ+a]*1/0.9)\n",
    "        ymax = mhist[typ+a]*1/0.9\n",
    "    try:\n",
    "        xo = np.mean(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "        ind = np.sum(hist[typ+'o'][1]<xo)-1\n",
    "        ax_s.plot([xo,xo],[0, ymax*0.249 ],\"-\",color=color_no,lw=1.2,zorder=1.1)      \n",
    "        ax_s.plot([xo,xo],[0, ymax*0.25],\"-\",lw=2,color='.15',zorder=1)\n",
    "    except:\n",
    "        pass\n",
    "    aco = sdata[typ+'_p'].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    acn = sdata[typ+'_p'].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    ax_s.set_title(r'{}'.format(name)  )\n",
    "    xf = np.mean(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "    ind = np.sum(hist[typ+'f'][1]<xf)-1\n",
    "    ax_s.plot([xf,xf],[0, ymax*0.249 ],\"-\",color=color_pr,lw=1.2,zorder=1.1)      \n",
    "    ax_s.plot([xf,xf],[0, ymax*0.25 ],\"-\",color='.15',lw=2,zorder=1)        \n",
    "    #ax_s.plot([p[typ],p[typ]],[0,ymax*0.2],\"-\",color='gray',lw=4,zorder=8)\n",
    "    ax_s.set_yticks(np.linspace(0,1,5)*ymax)\n",
    "    ax_s.set_yticklabels([])\n",
    "    \n",
    "    \n",
    "\n",
    " \n",
    "\n",
    "data,sdata = calc_statistic(data_s,p)\n",
    "\n",
    "\n",
    "\n",
    "hist = dict()\n",
    "mhist = dict()\n",
    "scalef = []\n",
    "for typ,ax_s in zip(types,axes_1):\n",
    "        try:\n",
    "            hist[typ+\"o\"] = np.histogram(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "        except:\n",
    "            pass\n",
    "        hist[typ+\"f\"] = np.histogram(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],density=True)\n",
    "    \n",
    "\n",
    "    \n",
    "for typ,ax_s in zip(types,axes_1):\n",
    "        try:\n",
    "            bins = np.linspace(hist[typ+'o'][1][0],hist[typ+'o'][1][-1],len(hist[typ+'o'][1])*2)\n",
    "            hist[typ+\"o\"] = ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                      bins=bins,alpha=alpha_no ,density=True,color=color_no,ec=color_no,linewidth=0.5)\n",
    "            ax_s.hist(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                      bins=bins,density=True,histtype='step',color=color_no,linewidth=1.0)\n",
    "            mhist[typ+\"o\"] = hist[typ+\"o\"][0].max()\n",
    "        except:\n",
    "            mhist[typ+\"o\"] = 0\n",
    "            pass\n",
    "        bins = np.linspace(hist[typ+'f'][1][0],hist[typ+'f'][1][-1],len(hist[typ+'f'][1])*2)\n",
    "        ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                  bins=bins,density=True,histtype='step',color=color_pr,linewidth=1.0)\n",
    "        hist[typ+\"f\"] = ax_s.hist(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]],\\\n",
    "                                  bins=bins,alpha=alpha_pr ,density=True,color=color_pr,ec=color_pr,lw=0.5)\n",
    "\n",
    "        mhist[typ+\"f\"] = hist[typ+\"f\"][0].max()\n",
    "        \n",
    "        scalef.append([mhist[typ+\"f\"] < mhist[typ+\"o\"]])\n",
    "        \n",
    "scalef = np.max( scalef )\n",
    "for typ,ax_s,name in zip(types,axes_1,names):\n",
    "    if scalef:\n",
    "        a = 'o'\n",
    "        b = 'f'\n",
    "    else:\n",
    "        a = 'f'\n",
    "        b = 'o'\n",
    "    if mhist[typ+a] == 0:\n",
    "        ax_s.set_ylim(0,mhist[typ+b]*1/0.9)\n",
    "        ymax = mhist[typ+b]*1/0.9\n",
    "    else:\n",
    "        ax_s.set_ylim(0,mhist[typ+a]*1/0.9)\n",
    "        ymax = mhist[typ+a]*1/0.9\n",
    "    try:\n",
    "        xo = np.mean(data[typ].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "        ind = np.sum(hist[typ+'o'][1]<xo)-1\n",
    "        ax_s.plot([xo,xo],[0, ymax*0.249 ],\"-\",color=color_no,lw=1.2,zorder=1.1)      \n",
    "        ax_s.plot([xo,xo],[0, ymax*0.25],\"-\",lw=2,color='.15',zorder=1)\n",
    "    except:\n",
    "        pass\n",
    "    aco = sdata[typ+'_p'].loc[idx[B,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    acn = sdata[typ+'_p'].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]]\n",
    "    #ax_s.set_title(r'{}'.format(name)  )\n",
    "    xf = np.mean(data[typ].loc[idx[A,p['ss'],p['sc'],p['m'],p['GFf']]])\n",
    "    ind = np.sum(hist[typ+'f'][1]<xf)-1\n",
    "    ax_s.plot([xf,xf],[0, ymax*0.249 ],\"-\",color=color_pr,lw=1.2,zorder=1.1)      \n",
    "    ax_s.plot([xf,xf],[0, ymax*0.25 ],\"-\",color='.15',lw=2,zorder=1)      \n",
    "    #ax_s.plot([p[typ],p[typ]],[0,ymax*0.2],\"-\",color='gray',lw=4,zorder=8)\n",
    "    ax_s.set_yticks(np.linspace(0,1,5)*ymax)\n",
    "    ax_s.set_yticklabels([])\n",
    "    \n",
    "    \n",
    "\n",
    "    \n",
    "\n",
    "with pd.HDFStore('data.pandas','r',complevel=9) as store:\n",
    "     allimnew = store['s100_n10'].loc[idx[:,:0.2:,:],:].copy()\n",
    "\n",
    "allimnew.reset_index(inplace=True)\n",
    "allimnew.rename({\"level_0\":\"model\"}, axis=\"columns\",inplace=True)\n",
    "allimnew.set_index(['model','start','leng','N'],inplace=True)\n",
    "allimnew.sort_index(inplace=True)\n",
    "\n",
    "        \n",
    "im = plot2d_pandas(allimnew.loc[idx[B]],'Tc','start','leng',aspect='auto',vmin=0.75,vmax=1.25,\\\n",
    "        cmap=plt.cm.bwr,axis=ax_12)\n",
    "\n",
    "cbar = plt.colorbar(im,ax=ax_12)\n",
    "tl = cbar.ax.get_yticklabels()\n",
    "plt.draw()\n",
    "tl = [float(i.get_text().replace('$','')) for i in tl]\n",
    "\n",
    "cbar.ax.set_yticklabels(['{:2.0f}'.format( (i-1)*100 )+\"\\%\" for i in tl])\n",
    "\n",
    "im = plot2d_pandas(allimnew.loc[idx[A]],'Tc','start','leng',aspect='auto',vmin=0.75,vmax=1.25,\\\n",
    "        cmap=plt.cm.bwr,axis=ax_11)\n",
    "#plt.colorbar(im,ax=ax_11,extend='neither')\n",
    "            \n",
    "\n",
    "\n",
    "ax_11.set_ylabel('measurement duration')\n",
    "ax_12.set_ylabel('measurement duration')\n",
    "ax_11.set_xlabel('time of first measurement')\n",
    "ax_12.set_xlabel('time of first measurement')\n",
    "\n",
    "#add lines\n",
    "for azoom,aall in zip(axes_0,axes_1):\n",
    "    azoom.set_zorder(10)\n",
    "    aall.set_zorder(6)\n",
    "    zoom_effect02(azoom,aall,color='black',ls=':',zorder=0)\n",
    "\n",
    "    \n",
    "#add text\n",
    "ax_1.set_ylabel('probability density')\n",
    "bb = ax_1.get_position()\n",
    "atxt = plt.figtext(bb.x0 - ho_abc/fb,bb.y1 + 0.05/fh,'A - 100 samples',size='x-large',zorder=7)\n",
    "\n",
    "ax_6.set_ylabel('probability density')\n",
    "bb = ax_6.get_position()\n",
    "btxt = plt.figtext(bb.x0 - ho_abc/fb,bb.y1 + 0.05/fh,'B - 5 samples',size='x-large',zorder=7)\n",
    "btxt.set_bbox(dict(facecolor='white', alpha=1, edgecolor=None,pad=0))\n",
    "\n",
    "\n",
    "\n",
    "#ax_tmp = plt.subplot2grid(axs,(2,0))\n",
    "#ax_tmp.set_axis_off()\n",
    "#bb = ax_tmp.get_position()\n",
    "#plt.figtext(bb.x0 - ho_abc/fb,bb.y1 + 0.05/fh,'C',size='x-large')\n",
    "\n",
    "bb = ax_11.get_position()\n",
    "plt.figtext(bb.x0 - 0.7/fb,bb.y1 + 0.05/fh,'C',size='x-large')\n",
    "bb = ax_12.get_position()\n",
    "plt.figtext(bb.x0 - 0.7/fb,bb.y1 + 0.05/fh,'D',size='x-large')\n",
    "\n",
    "plt.draw()\n",
    "        \n",
    "tmp = ax_4.get_xticklabels()\n",
    "tmp[-2].set_horizontalalignment('right')\n",
    "ax_4.set_xticklabels(tmp)\n",
    "\n",
    "tmp = ax_5.get_xticklabels()\n",
    "tmp[1].set_horizontalalignment('left')\n",
    "ax_5.set_xticklabels(tmp)\n",
    "\n",
    "tmp = ax_8.get_xticklabels()\n",
    "tmp[-2].set_horizontalalignment('right')\n",
    "ax_8.set_xticklabels(tmp)\n",
    "\n",
    "for ax_s in axes_0:\n",
    "    for t in ax_s.get_xticklabels():\n",
    "        #try:\n",
    "        #    bbox = t.get_window_extent().transformed(ax_s.transData.inverted())\n",
    "        #except:\n",
    "        #    continue\n",
    "        #rect = mpl.patches.Rectangle((bbox.x0,bbox.y0), bbox.x1 - bbox.x0 ,bbox.y1 - bbox.y0 , ec=\"none\", color='white',zorder=0)\n",
    "        #rect.set_clip_on(False)\n",
    "        #ax_s.add_patch(rect)\n",
    "        t.set_bbox(dict(facecolor='white', alpha=1, edgecolor=None,pad=0))\n",
    "\n",
    "\n",
    "for typ,ax_s in zip(types,axes_0):\n",
    "    labels = []\n",
    "    trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "    for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "        la = l.get_text()\n",
    "        if np.abs(t - p[typ])<0.005 and typ!='sc':\n",
    "            lold = la.replace('$','')\n",
    "            #lfont = l.get_font_properties()\n",
    "            #lfont.set_weight('bold')\n",
    "            #l.set_fontproperties(font)\n",
    "            la = '$\\\\mathbf{'+lold+'}$'\n",
    "            line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "            ax_s.add_line(line)\n",
    "        labels.append(la)\n",
    "    ax_s.set_xticklabels(labels)  \n",
    "    \n",
    "for typ,ax_s in zip(types,axes_1):\n",
    "    labels = []\n",
    "    trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "    for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "        la = l.get_text()\n",
    "        if np.abs(t - p[typ])<0.005 and typ!='sc':\n",
    "            lold = la.replace('$','')\n",
    "            #lfont = l.get_font_properties()\n",
    "            #lfont.set_weight('bold')\n",
    "            #l.set_fontproperties(font)\n",
    "            la = '$\\\\mathbf{'+lold+'}$'\n",
    "            line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "            ax_s.add_line(line)\n",
    "        labels.append(la)\n",
    "    ax_s.set_xticklabels(labels)\n",
    "\n",
    "ax_s = axes_0[-1]\n",
    "trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "found=False\n",
    "for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "    la = l.get_text()\n",
    "    if np.abs(t - clapy.calc_sigma_true(p[typ],0.5,0.2))<0.005:\n",
    "        lold = la.replace('$','')\n",
    "        #lfont = l.get_font_properties()\n",
    "        #lfont.set_weight('bold')\n",
    "        #l.set_fontproperties(font)\n",
    "        #a = '$\\\\mathbf{'+lold+'}$'\n",
    "        a = '$\\\\math{'+lold+'}$'\n",
    "        line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "        ax_s.add_line(line)\n",
    "        found=True\n",
    "    labels.append(la)\n",
    "nt = clapy.calc_sigma_true(p[typ],0.5,0.2)\n",
    "if not found:\n",
    "    line = pltlines.Line2D([nt,nt],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "    ax_s.add_line(line)\n",
    "\n",
    "ax_s = axes_1[-1]\n",
    "trans = plttransforms.blended_transform_factory(ax_s.transData, ax_s.transAxes)\n",
    "found=False\n",
    "for l,t in zip(ax_s.get_xticklabels(),ax_s.get_xticks()):\n",
    "    la = l.get_text()\n",
    "    if np.abs(t - clapy.calc_sigma_true(p[typ],0.5,0.2))<0.005:\n",
    "        lold = la.replace('$','')\n",
    "        #lfont = l.get_font_properties()\n",
    "        #lfont.set_weight('bold')\n",
    "        #l.set_fontproperties(font)\n",
    "        #a = '$\\\\mathbf{'+lold+'}$'\n",
    "        a = '$\\\\math{'+lold+'}$'\n",
    "        line = pltlines.Line2D([t,t],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "        ax_s.add_line(line)\n",
    "        found=True\n",
    "    labels.append(la)\n",
    "nt = clapy.calc_sigma_true(p[typ],0.5,0.2)\n",
    "if not found:\n",
    "    line = pltlines.Line2D([nt,nt],[0,-0.045],color='0.15',transform=trans,clip_on=False,lw=2,solid_capstyle='butt' )\n",
    "    ax_s.add_line(line)    \n",
    "    \n",
    "    \n",
    "bb = ax_2.get_position()\n",
    "lax = plt.axes([bb.x0 - 0.2/fb,bb.y1 + 0.2/fh ,0.4,0.014])\n",
    "lax.set_xticks([])\n",
    "lax.set_yticks([])\n",
    "lax.set_frame_on(False)\n",
    "lax.set_xlim(0,100)\n",
    "a = 20\n",
    "a3 = 60\n",
    "s = 1.2\n",
    "at= 4\n",
    "\n",
    "lax.text(a+at,0,'stochastic model',size='large')\n",
    "lax.bar(a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.5,color='w',edgecolor=color_pr,alpha=alpha_pr)\n",
    "lax.plot(a+s*np.array([-0.50,-0.5,0.5,0.5,1.5,1.5]),[0,2.8,2.8,5,5,0],color=color_pr)\n",
    "lax.bar(a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.0,color=color_pr,edgecolor=None,alpha=alpha_pr)\n",
    "\n",
    "lax.text(a3+a+at,0,'Nowakowski model',size='large')\n",
    "lax.bar(a3+a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.5,color='w',edgecolor=color_no,alpha=alpha_no)\n",
    "lax.plot(a3+a+s*np.array([-0.50,-0.5,0.5,0.5,1.5,1.5]),[0,2.8,2.8,5,5,0],color=color_no)\n",
    "lax.bar(a3+a+s*np.array([0,1]),[2.8,5],s*1,linewidth=0.0,color=color_no,edgecolor=None,alpha=alpha_no)\n",
    "\n",
    "bb = ax_12.get_position()\n",
    "ax_11.set_position([ax_11.get_position().x0,bb.y0,bb.width,bb.height])\n",
    "\n",
    "ax_11.set_title('Bias stochastic model')\n",
    "ax_12.set_title('Bias Nowakowski model')\n",
    "\n",
    "\n",
    "plt.draw()\n",
    "f.savefig(\"5_gamma.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "12px",
    "width": "252px"
   },
   "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": "block",
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}