paper_simulations.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import pandas as pd\n",
    "import warnings\n",
    "\n",
    "import cython\n",
    "%load_ext Cython\n",
    "from iminuit import Minuit\n",
    "idx = pd.IndexSlice\n",
    "\n",
    "import clapy\n",
    "import clasim\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Data for time series and histogramms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "asym_dist = clapy.dist()\n",
    "\n",
    "dargs = {\n",
    "        'samples': 10000,\n",
    "        'nCells': 100,\n",
    "        'mCells': 100,\n",
    "        'GF': 0.95,\n",
    "        'G1': 0.46,\n",
    "        'S': 0.33,\n",
    "        'G2M': 0.21,\n",
    "        'sCells' : 0.3,\n",
    "        'sSamples' : 0.2\n",
    "}\n",
    "dargs_d = dargs.copy()\n",
    "dargs_d['sSamples'] = 1e-6\n",
    "dargs_d['sCells'] = 1e-6\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",
    "\n",
    "\n",
    "pdfs = list()\n",
    "pdfm = list()\n",
    "data_asy = list()\n",
    "data_sym = list()\n",
    "\n",
    "data_asy_d = list()\n",
    "data_sym_d = list()\n",
    "\n",
    "for t in time_points:\n",
    "    #simulations asymetric\n",
    "    data_asy.append(np.array(clasim.run(seed=int(np.random.rand()*1000),mode=1,times=[t],**dargs)) )\n",
    "    data_asy_d.append(np.array(clasim.run(seed=int(np.random.rand()*1000),mode=1,times=[t],**dargs_d)) )\n",
    "    #simulations symetric\n",
    "    data_sym.append(np.array(clasim.run(seed=int(np.random.rand()*1000),mode=3,times=[t],**dargs)) )\n",
    "    data_sym_d.append(np.array(clasim.run(seed=int(np.random.rand()*1000),mode=3,times=[t],**dargs_d)) )\n",
    "    #pdfs\n",
    "    pdfs.append( [asym_dist.pmf_f(dargs['nCells'],dTC,dFS,dargs['GF'],dargs['sCells'],dargs['sSamples'],t,i) for i in X] )\n",
    "    #means\n",
    "    pdfm.append( asym_dist.pmf_mean(dargs['nCells'],dTC,dFS,dargs['GF'],dargs['sCells'],dargs['sSamples'],t) )\n",
    "\n",
    "    \n",
    "with pd.HDFStore('paper_data.pandas',complevel=9) as st:\n",
    "    st['data_asy'] = pd.Series(data_asy)\n",
    "    st['data_sym'] = pd.Series(data_sym)\n",
    "    \n",
    "    st['data_asy_d'] = pd.Series(data_asy_d)\n",
    "    st['data_sym_d'] = pd.Series(data_sym_d)\n",
    "    \n",
    "    st['pdfs']     = pd.Series(pdfs)\n",
    "    st['pdfm']     = pd.Series(pdfm)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gernerate data for parameter recovery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython --force\n",
    "# distutils: language = c++\n",
    "\n",
    "#use --annotate if you wonder what kind of code it generates \n",
    "cimport cython\n",
    "import numpy as np\n",
    "cimport numpy as np #overwritten those from python with cython\n",
    "from libc.math cimport exp, M_PI, sqrt, log\n",
    "from iminuit.util import describe, make_func_code\n",
    "from libcpp.map cimport map     \n",
    "import scipy as sp\n",
    "\n",
    "\n",
    "\n",
    "@cython.embedsignature(True)#dump the signatre so describe works\n",
    "cpdef fitfunc(double t,double Tc,double r,double GF):\n",
    "    cdef double res = 0\n",
    "    if t<(Tc-Tc*r):\n",
    "        res = GF/Tc*(t+Tc*r)\n",
    "    else:\n",
    "        res =  GF\n",
    "    return res \n",
    "'''\n",
    "@cython.embedsignature(True)#dump the signatre so describe works\n",
    "cpdef double mypdf(double x, double mu, double std):\n",
    "    #cpdef means generate both c function and python function\n",
    "    cdef double norm = 1./(sqrt(2*M_PI*std))\n",
    "    cdef double ret = exp(-1*(x-mu)*(x-mu)/(2.*std))*norm\n",
    "    return ret\n",
    "'''\n",
    "@cython.embedsignature(True)#dump the signatre so describe works\n",
    "cpdef double mypdfln(double x, double mu, double std):\n",
    "    #cpdef means generate both c function and python function\n",
    "    cdef double norm = (sqrt(2*M_PI*std*std))\n",
    "    cdef double ret = (-1*(x-mu)*(x-mu)/(2.*std*std))-log(norm)\n",
    "    return ret\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "cdef class Nowakowski_LH:\n",
    "    cdef np.ndarray data\n",
    "    cdef np.ndarray err\n",
    "    cdef np.ndarray t\n",
    "    cdef int ndata\n",
    "    \n",
    "    def __init__(self,data,t):\n",
    "        self.data = data\n",
    "        self.t = t\n",
    "        self.ndata = len(data)\n",
    "    \n",
    "    @cython.embedsignature(True)#you need this to dump function signature in docstring\n",
    "    def compute(self, double Tc,double r,double GF,double s):\n",
    "        #this line is a cast not a copy. Let cython knows mydata will spit out double\n",
    "        cdef np.ndarray[np.double_t, ndim=1] mydata = self.data\n",
    "        cdef np.ndarray[np.double_t, ndim=1] myt = self.t\n",
    "        cdef double loglh = 0.\n",
    "        cdef double lmu = 0.\n",
    "        cdef double 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"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## with 5 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_args = {\n",
    "        'samples': 5,\n",
    "        'nCells': 100,\n",
    "        'mCells': 100,\n",
    "        'GF': 0.95,\n",
    "        'G1': 0.46,\n",
    "        'S': 0.33,\n",
    "        'G2M': 0.21,\n",
    "        'sCells' : 0.3,\n",
    "        'sSamples' : 0.2 \n",
    "}\n",
    "warnings.filterwarnings('ignore')\n",
    "calc_error = False\n",
    "times=np.linspace(0.01,1.5,5)\n",
    "default_args['times'] = times\n",
    "#default_args['samples'] = int(np.round(120/(mm)))\n",
    "\n",
    "\n",
    "index = pd.MultiIndex.from_product([[],[]],names=['std','index'])\n",
    "# Priors for unknown model parameters\n",
    "\n",
    "tmplist = []\n",
    "tmplist2= []\n",
    "for num in range(1000):\n",
    "    Y,X = np.array(clasim.run(seed=2*num,mode=1,**default_args))\n",
    "    \n",
    "\n",
    "    lh = Nowakowski_LH(Y*1.0/default_args['nCells'],X)\n",
    "    #lh.compute(0.4,5,3)\n",
    "\n",
    "    mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=0.95,s=0.2,\\\n",
    "               error_s=0.1,error_r=0.1,error_GF=0.1,error_Tc=0.1,\\\n",
    "               limit_Tc=(0,2), limit_r=(0,1),limit_GF=(0,1),limit_s=(0,1),\\\n",
    "               errordef=0.5,print_level=0)\n",
    "    mi_old.migrad(ncall=999999999);\n",
    "    s = dict(mi_old.values)\n",
    "    for ii in mi_old.errors:\n",
    "        s.update({ii+\"_error\" : mi_old.errors[ii]})\n",
    "    s.update({'model' : 'old'})\n",
    "\n",
    "    if calc_error:\n",
    "        try:\n",
    "            tmp = mi_old.minos(sigma=2)                    \n",
    "            for i in tmp:\n",
    "                s.update({i+'_hpd46_l' :  tmp[i]['lower'] + mi_old.values[i] }) \n",
    "                s.update({i+'_hpd46_u' :  tmp[i]['upper'] + mi_old.values[i] }) \n",
    "        except:\n",
    "            print('error in errorfind',s['model'])\n",
    "\n",
    "    tmplist2.append(num)\n",
    "    tmplist.append(s)\n",
    "\n",
    "\n",
    "\n",
    "    Y,X = np.array(clasim.run(seed=2*num+1,mode=1,**default_args))\n",
    "\n",
    "    lh = clapy.asym_lh(Y,X,100)\n",
    "    mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=0.95,sigma_cell=0.3,sigma_sample=0.2,\\\n",
    "       error_sigma_cell=0.1,error_r=0.1,error_GF=0.1,error_Tc=0.1,error_sigma_sample=0.1,\\\n",
    "       limit_Tc=(0.00001,2), limit_r=(0.00001,1),limit_GF=(0,1),limit_sigma_cell=(0.00001,1),limit_sigma_sample=(0.00001,1),\\\n",
    "       errordef=0.5,print_level=0)\n",
    "    mi_old.migrad();\n",
    "    s = dict(mi_old.values)\n",
    "    for ii in mi_old.errors:\n",
    "        s.update({ii+\"_error\" : mi_old.errors[ii]})\n",
    "    s.update({'model' : 'full'})\n",
    "\n",
    "    if calc_error:\n",
    "        try:\n",
    "            tmp = mi_old.minos(sigma=2)                    \n",
    "            for i in tmp:\n",
    "                s.update({i+'_hpd46_l' :  tmp[i]['lower'] + mi_old.values[i] }) \n",
    "                s.update({i+'_hpd46_u' :  tmp[i]['upper'] + mi_old.values[i] }) \n",
    "        except:\n",
    "            print('error in errorfind',s['model'])\n",
    "\n",
    "\n",
    "    tmplist2.append(num)\n",
    "    tmplist.append(s)   \n",
    "    print('finiscged',num,s['model'],\"t\")\n",
    "\n",
    "stat5 = pd.DataFrame(tmplist,index=tmplist2)\n",
    "stat5.rename_axis({\"GF\":\"GFf\"}, axis=\"columns\",inplace=True)\n",
    "stat5.index.rename('N',inplace=True)\n",
    "stat5['sSamples'] = 0.2\n",
    "stat5['sCells'] = 0.3\n",
    "stat5.set_index('sSamples', append=True, inplace=True)\n",
    "stat5.set_index('sCells', append=True, inplace=True)\n",
    "stat5.set_index('model', append=True, inplace=True)    \n",
    "\n",
    "stat5 =  stat5.reorder_levels(['model','sSamples','sCells','N']) \n",
    "\n",
    "warnings.filterwarnings('default')\n",
    "\n",
    "stat5 = stat5.sort_index()\n",
    "\n",
    "with pd.HDFStore('paper_data.pandas',complib='zlib',complevel=9) as st:\n",
    "    st['minuit_full_new_s5_para'] = default_args\n",
    "    st['minuit_full_new_s5'] = stat5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## with 100 samples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_args = {\n",
    "        'samples': 100,\n",
    "        'nCells': 100,\n",
    "        'mCells': 100,\n",
    "        'GF': 0.95,\n",
    "        'G1': 0.46,\n",
    "        'S': 0.33,\n",
    "        'G2M': 0.21,\n",
    "        'sCells' : 0.3,\n",
    "        'sSamples' : 0.2 \n",
    "}\n",
    "warnings.filterwarnings('ignore')\n",
    "calc_error = False\n",
    "times=np.linspace(0.01,1.5,5)\n",
    "default_args['times'] = times\n",
    "\n",
    "#default_args['samples'] = int(np.round(120/(mm)))\n",
    "\n",
    "\n",
    "index = pd.MultiIndex.from_product([[],[]],names=['std','index'])\n",
    "# Priors for unknown model parameters\n",
    "SM = 456426\n",
    "tmplist = []\n",
    "tmplist2= []\n",
    "for num in range(1000):\n",
    "    Y,X = np.array(clasim.run(seed=num*2+SM,mode=1,**default_args))\n",
    "\n",
    "    lh = Nowakowski_LH(Y*1.0/default_args['nCells'],X)\n",
    "    #lh.compute(0.4,5,3)\n",
    "\n",
    "    mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=0.95,s=0.2,\\\n",
    "               error_s=0.1,error_r=0.1,error_GF=0.1,error_Tc=0.1,\\\n",
    "               limit_Tc=(0,2), limit_r=(0,1),limit_GF=(0,1),limit_s=(0,1),\\\n",
    "               errordef=0.5,print_level=0)\n",
    "    mi_old.migrad(ncall=999999999);\n",
    "    s = dict(mi_old.values)\n",
    "    for ii in mi_old.errors:\n",
    "        s.update({ii+\"_error\" : mi_old.errors[ii]})\n",
    "    s.update({'model' : 'old'})\n",
    "\n",
    "    if calc_error:\n",
    "        try:\n",
    "            tmp = mi_old.minos(sigma=2)                    \n",
    "            for i in tmp:\n",
    "                s.update({i+'_hpd46_l' :  tmp[i]['lower'] + mi_old.values[i] }) \n",
    "                s.update({i+'_hpd46_u' :  tmp[i]['upper'] + mi_old.values[i] }) \n",
    "        except:\n",
    "            print('error in errorfind',s['model'])\n",
    "\n",
    "    tmplist2.append(num)\n",
    "    tmplist.append(s)\n",
    "\n",
    "\n",
    "\n",
    "    Y,X = np.array(clasim.run(seed=num*2+SM+1,mode=1,**default_args))\n",
    "\n",
    "\n",
    "    lh = clapy.asym_lh(Y,X,100)\n",
    "    mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=0.95,sigma_cell=0.3,sigma_sample=0.2,\\\n",
    "       error_sigma_cell=0.1,error_r=0.1,error_GF=0.1,error_Tc=0.1,error_sigma_sample=0.1,\\\n",
    "       limit_Tc=(0.00001,2), limit_r=(0.00001,1),limit_GF=(0,1),limit_sigma_cell=(0.00001,1),limit_sigma_sample=(0.00001,1),\\\n",
    "       errordef=0.5,print_level=0)\n",
    "    mi_old.migrad();\n",
    "    s = dict(mi_old.values)\n",
    "    for ii in mi_old.errors:\n",
    "        s.update({ii+\"_error\" : mi_old.errors[ii]})\n",
    "    s.update({'model' : 'full'})\n",
    "\n",
    "    if calc_error:\n",
    "        try:\n",
    "            tmp = mi_old.minos(sigma=2)                    \n",
    "            for i in tmp:\n",
    "                s.update({i+'_hpd46_l' :  tmp[i]['lower'] + mi_old.values[i] }) \n",
    "                s.update({i+'_hpd46_u' :  tmp[i]['upper'] + mi_old.values[i] }) \n",
    "        except:\n",
    "            print('error in errorfind',s['model'])\n",
    "\n",
    "\n",
    "    tmplist2.append(num)\n",
    "    tmplist.append(s)   \n",
    "    print('finiscged',num,s['model'],\"t\")\n",
    "\n",
    "stat100 = pd.DataFrame(tmplist,index=tmplist2)\n",
    "stat100.rename_axis({\"GF\":\"GFf\"}, axis=\"columns\",inplace=True)\n",
    "stat100.index.rename('N',inplace=True)\n",
    "stat100['sSamples'] = 0.2\n",
    "stat100['sCells'] = 0.3\n",
    "stat100.set_index('sSamples', append=True, inplace=True)\n",
    "stat100.set_index('sCells', append=True, inplace=True)\n",
    "stat100.set_index('model', append=True, inplace=True)    \n",
    "\n",
    "stat100 =  stat100.reorder_levels(['model','sSamples','sCells','N']) \n",
    "\n",
    "warnings.filterwarnings('default')\n",
    "\n",
    "stat100 = stat100.sort_index()\n",
    "\n",
    "with pd.HDFStore('paper_data.pandas',complib='zlib',complevel=9) as st:\n",
    "    st['minuit_full_new_s100_para'] = default_args\n",
    "    st['minuit_full_new_s100'] = stat100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## for different initial conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "default_args = {\n",
    "        'samples': 100,\n",
    "        'nCells': 100,\n",
    "        'mCells': 100,\n",
    "        'GF': 0.95,\n",
    "        'G1': 0.46,\n",
    "        'S': 0.33,\n",
    "        'G2M': 0.21,\n",
    "        'sCells' : 0.3,\n",
    "        'sSamples' : 0.2 \n",
    "}\n",
    "warnings.filterwarnings('ignore')\n",
    "calc_error = False\n",
    "allim = []\n",
    "sss = 0\n",
    "SM = 928639\n",
    "for start in np.linspace(0.01,0.2,19*1+1): #range(5,31,5):\n",
    "    for leng in np.linspace(0.5,2,15*1+1):\n",
    "        for n in range(10):\n",
    "            times=np.linspace(start,start+leng,5)\n",
    "            Y,X = np.array(clasim.run(seed=sss+SM,mode=1,times=times,**default_args))\n",
    "            sss = sss+1\n",
    "            lh = Nowakowski_LH(Y*1.0/default_args['nCells'],X)\n",
    "            \n",
    "\n",
    "            mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=0.95,s=0.2,\\\n",
    "                       error_s=0.1,error_r=0.1,error_GF=0.1,error_Tc=0.1,\\\n",
    "                       limit_Tc=(0,2), limit_r=(0,1),limit_GF=(0,1),limit_s=(0,1),\\\n",
    "                       errordef=0.5,print_level=0)\n",
    "            mi_old.migrad(ncall=999999999);\n",
    "            s = dict(mi_old.values)\n",
    "            for ii in mi_old.errors:\n",
    "                s.update({ii+\"_error\" : mi_old.errors[ii]})\n",
    "            s.update({'N' : n})\n",
    "            if calc_error:\n",
    "                try:\n",
    "                    tmp = mi_old.minos(sigma=2)                    \n",
    "                    for i in tmp:\n",
    "                        s.update({i+'_hpd46_l' :  tmp[i]['lower'] + mi_old.values[i] }) \n",
    "                        s.update({i+'_hpd46_u' :  tmp[i]['upper'] + mi_old.values[i] }) \n",
    "                except:\n",
    "                    print('error in errorfind',s['model'])\n",
    "            s.update({'leng' : leng})\n",
    "            s.update({'start' : start})\n",
    "            nowak = s.copy()\n",
    "\n",
    "\n",
    "\n",
    "            Y,X = np.array(clasim.run(seed=sss+SM,mode=1,times=times,**default_args))\n",
    "            sss = sss+1\n",
    "            \n",
    "            lh = clapy.asym_lh(Y,X,100)\n",
    "            mi_old = Minuit(lh.compute, Tc=1.0, r=0.3,GF=0.95,sigma_cell=0.3,sigma_sample=0.2,\\\n",
    "               error_sigma_cell=0.1,error_r=0.1,error_GF=0.1,error_Tc=0.1,error_sigma_sample=0.1,\\\n",
    "               limit_Tc=(0.00001,2), limit_r=(0.00001,1),limit_GF=(0,1),limit_sigma_cell=(0.00001,1),limit_sigma_sample=(0.00001,1),\\\n",
    "               errordef=0.5,print_level=0)\n",
    "            mi_old.migrad();\n",
    "            s = dict(mi_old.values)\n",
    "            for ii in mi_old.errors:\n",
    "                s.update({ii+\"_error\" : mi_old.errors[ii]})\n",
    "            s.update({'N' : n})\n",
    "            s.update({'leng' : leng})\n",
    "            s.update({'start' : start})\n",
    "\n",
    "            if calc_error:\n",
    "                try:\n",
    "                    tmp = mi_old.minos(sigma=2)                    \n",
    "                    for i in tmp:\n",
    "                        s.update({i+'_hpd46_l' :  tmp[i]['lower'] + mi_old.values[i] }) \n",
    "                        s.update({i+'_hpd46_u' :  tmp[i]['upper'] + mi_old.values[i] }) \n",
    "                except:\n",
    "                    print('error in errorfind',s['model'])\n",
    "\n",
    "\n",
    "            print('finiscged')\n",
    "            print(sss+SM,nowak,s)\n",
    "            stat = pd.DataFrame([nowak,s],index=['old','full'])\n",
    "            stat.set_index('start', append=True, inplace=True)\n",
    "            stat.set_index('leng', append=True, inplace=True)\n",
    "            stat.set_index('N', append=True, inplace=True)    \n",
    "\n",
    "            allim.append( stat )\n",
    "            #reorder_levels(['start','leng','model']) )\n",
    "\n",
    "warnings.filterwarnings('default')\n",
    "\n",
    "allimnew = pd.concat(allim).sort_index()\n",
    "\n",
    "with pd.HDFStore('paper_data.pandas',complib='zlib',complevel=9) as st:\n",
    "    st['s100_n10'] = allimnew\n",
    "    st['s100_n10_para'] = default_args"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allimnew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {
    "height": "49.1333px",
    "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
}