siRNA_knockdown_calibrated_nonfixed.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## siRNA knockdown calibrated ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook fits analytical functions to the gene silencing traces without fixing the protein maturation and degradation rates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Motivation\n",
    "siRNA (small interfering RNA) triggers specific degradation of mRNA. The RISC (RNA-induced silencing complex), which consists of siRNA and some proteins, cuts mRNA containing a strand sequence complementary to the sequence of the siRNA. Reducing protein expression by adding siRNA is called gene knockdown.\n",
    "\n",
    "Gene knockdown is a promising approach for the treatment of some diseases, e.g. cancer. The aim of this project is to study the influence of siRNA on gene expression and mRNA degradation on the single cell level to promote development of siRNA-based medical treatments.\n",
    "\n",
    "This is done by fitting the solutions of the differential equations describing the expression network to measured fluorescence traces of cells transfected with a GFP mRNA and a RFP mRNA, where siRNA specific for GFP mRNA is added and RFP is used as a reference.\n",
    "\n",
    "Among the fit parameters, there is the GFP mRNA degradation rate $\\delta_\\text{g}$ and the RFP mRNA degradation rate $\\delta_\\text{r}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solutions look like this:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f_\\text{red}(t) =\n",
    "m_\\text{r}\\,k_\\text{tl} \\left(\n",
    "\\frac{1}{\\beta_\\text{r}-\\delta_\\text{r}+k_\\text{m,r}} \\mathrm{e}^{-(\\beta_\\text{r}+k_\\text{m,r})(t-t_0)}\n",
    "-\\frac{1}{\\beta_\\text{r} - \\delta_\\text{r}} \\mathrm{e}^{-\\beta_\\text{r} (t-t_0)}\n",
    "+\\frac{k_\\text{m,r}}{(\\beta_\\text{r}-\\delta_\\text{r}) (\\beta_\\text{r}-\\delta_\\text{r}+k_\\text{m,r})} \\mathrm{e}^{-\\delta_\\text{r} (t-t_0)}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$f_\\text{green}(t) =\n",
    "m_\\text{g}\\,k_\\text{tl} \\left(\n",
    "\\frac{1}{\\beta_\\text{g}-\\delta_\\text{g}+k_\\text{m,g}} \\mathrm{e}^{-(\\beta_\\text{g}+k_\\text{m,g})(t-t_0)}\n",
    "-\\frac{1}{\\beta_\\text{g} - \\delta_\\text{g}} \\mathrm{e}^{-\\beta_\\text{g} (t-t_0)}\n",
    "+\\frac{k_\\text{m,g}}{(\\beta_\\text{g}-\\delta_\\text{g}) (\\beta_\\text{g}-\\delta_\\text{g}+k_\\text{m,g})} \\mathrm{e}^{-\\delta_\\text{g} (t-t_0)}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Notebook structure\n",
    "The notebook has the following structure:\n",
    "\n",
    "At first, the model functions are defined and the data is loaded. The next section contains code for fitting the two models separately. The next section contains code for fitting the two traces in one run with parameters shared among the models.\n",
    "\n",
    "Fitting requires that the result list `R` is defined, which can be done by running the corresponding cell. When `R` has been populated by fitting, the results can be plotted. There are cells for plotting the results of the separate fit, the results of the combined fit, and the pure parameter distributions of all fits.\n",
    "\n",
    "Additionally, there are cells for saving and loading paramaters by python’s `pickle` module."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules needed\n",
    "\n",
    "# Standard library\n",
    "from collections import OrderedDict\n",
    "from copy import deepcopy\n",
    "import inspect\n",
    "import os\n",
    "import pickle\n",
    "import re\n",
    "import sys\n",
    "\n",
    "# Scientific stack\n",
    "import numpy as np\n",
    "np.seterr(divide='print')\n",
    "import pandas as pd\n",
    "import scipy as sc\n",
    "import scipy.optimize as so\n",
    "import scipy.stats as ss\n",
    "from sklearn.cluster import DBSCAN\n",
    "from sklearn.neighbors import KernelDensity\n",
    "\n",
    "# Matplotlib\n",
    "%matplotlib inline\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "import matplotlib.cm\n",
    "from matplotlib.gridspec import GridSpec\n",
    "import matplotlib.lines as mlin\n",
    "import matplotlib.patches as mptch\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Notebook utilities\n",
    "import IPython\n",
    "import ipywidgets as wdg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define utility functions\n",
    "def getTimeStamp():\n",
    "    \"\"\"Returns a human-readable string representation of the current time\"\"\"\n",
    "    import time\n",
    "    from datetime import datetime\n",
    "    return datetime.now().strftime(\"%Y-%m-%d–%H%M%S\")\n",
    "\n",
    "\n",
    "def getOutpath(filename='', timestamp=None):\n",
    "    \"\"\"Returns (and creates, if necessary) the path to a directory\n",
    "    called “out” inside the current directory.\n",
    "    If `filename` is given, the filename is appended to the output directory.\n",
    "    A timestamp will be added to the filename if `timestamp != ''`.\n",
    "    If timestamp is `None`, the current timestamp is used.\n",
    "    \"\"\"\n",
    "    # Create output directory\n",
    "    outpath = os.path.join(os.getcwd(), 'out')\n",
    "    if not os.path.isdir(outpath) and not os.path.lexists(outpath):\n",
    "        os.mkdir(outpath)\n",
    "\n",
    "    # If requested, build filename\n",
    "    if len(filename) > 0:\n",
    "        if timestamp == None:\n",
    "            timestamp = getTimeStamp()\n",
    "        outpath = os.path.join(outpath, ((timestamp + '_') if len(timestamp) > 0 else '') + filename)\n",
    "    return outpath"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(t, t0, m, k, beta, delta, offset):\n",
    "    \"\"\"General expression model function\"\"\"\n",
    "\n",
    "    f = np.zeros(np.shape(t))\n",
    "    idx_after = (t > t0)\n",
    "    dt = t[idx_after] - t0\n",
    "\n",
    "    f1 = np.exp(- (beta + k) * dt) / (beta - delta + k)\n",
    "    f2 = - np.exp(- beta * dt) / (beta - delta)\n",
    "    f3 = k * np.exp(- delta * dt) / (beta - delta) / (beta - delta + k)\n",
    "\n",
    "    f[idx_after] = (f1 + f2 + f3) * m\n",
    "\n",
    "    return f + offset\n",
    "\n",
    "def red(t, tr, m_ktl, kmr, betr, deltr, offr):\n",
    "    \"\"\"Model function for RFP data\"\"\"\n",
    "    return model(t=t, t0=tr, m=m_ktl, k=kmr, beta=betr, delta=deltr, offset=offr)\n",
    "\n",
    "def green(t, tg, m_ktl, kmg, betg, deltg, offg):\n",
    "    \"\"\"Model function for GFP data\"\"\"\n",
    "    return model(t=t, t0=tg, m=m_ktl, k=kmg, beta=betg, delta=deltg, offset=offg)\n",
    "\n",
    "def combined(t, tr, tg, m_ktl, kmr, kmg, betr, betg, deltr, deltg, offr, offg):\n",
    "    \"\"\"Model function for a combined fit of red and green data\"\"\"\n",
    "    return np.stack(\n",
    "        (model(t=t, t0=tr, m=m_ktl, k=kmr, beta=betr, delta=deltr, offset=offr),\n",
    "         model(t=t, t0=tg, m=m_ktl, k=kmg, beta=betg, delta=deltg, offset=offg)),\n",
    "        axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Set default parameter values\n",
    "m_ktl_0 = 2e4\n",
    "\n",
    "tr_0 = 4.5\n",
    "kmr_0 = 0.30\n",
    "betr_0 = 4.21e-2\n",
    "deltr_0 = 0.03\n",
    "offr_0 = 0\n",
    "\n",
    "tg_0 = 2\n",
    "kmg_0 = 1.28\n",
    "betg_0 = 5.22e-3\n",
    "deltg_0 = 2\n",
    "offg_0 = 0\n",
    "\n",
    "MAX_m_ktl = 5e5\n",
    "MAX_tr = 30\n",
    "MAX_tg = 30\n",
    "MAX_kmr = 30\n",
    "MAX_kmg = 30\n",
    "MAX_betr = 10\n",
    "MAX_betg = 10\n",
    "MAX_deltr = .2\n",
    "MAX_deltg = 10\n",
    "\n",
    "MIN_m_ktl = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class FitParameters:\n",
    "    \"\"\"FitParameters facilitates managing values and bounds of fit parameters\"\"\"\n",
    "    def __init__(self, fun, independent=[], fixed=[]):\n",
    "        # Store function\n",
    "        self.fun = fun\n",
    "\n",
    "        # Get parameters of fun\n",
    "        params = inspect.signature(self.fun).parameters\n",
    "\n",
    "        # Build data frame of parameters\n",
    "        self.df = pd.DataFrame(columns=['value', 'min', 'max'],\n",
    "                               index=[p for p in params.keys()],\n",
    "                               dtype=np.float64)\n",
    "\n",
    "        # Set “independent” and “fixed” flag\n",
    "        self.df.add(pd.DataFrame(columns=['independent', 'fixed'], dtype=np.bool))\n",
    "        for p in self.df.index.values:\n",
    "            self.df.loc[p, 'independent'] = p in independent\n",
    "            self.df.loc[p, 'fixed'] = p in fixed\n",
    "\n",
    "        # Set default parameters\n",
    "        for p in self.df.index.values:\n",
    "            if params[p].default == inspect.Parameter.empty:\n",
    "                if self.df.loc[p, 'independent']:\n",
    "                    self.df.loc[p, 'value'] = np.NaN\n",
    "                else:\n",
    "                    self.df.loc[p, 'value'] = 0\n",
    "            else:\n",
    "                self.df.loc[p, 'value'] = params[p].default\n",
    "\n",
    "    def set(self, p, **props):\n",
    "        \"\"\"Allows user to change parameter properties\"\"\"\n",
    "        if p not in self.df.index.values:\n",
    "            raise KeyError(\"Unknown parameter name: {}\".format(p))\n",
    "\n",
    "        for prop, val in props.items():\n",
    "            if prop == 'value':\n",
    "                self.df.loc[p, 'value'] = val\n",
    "            elif prop == 'min':\n",
    "                self.df.loc[p, 'min'] = val\n",
    "            elif prop == 'max':\n",
    "                self.df.loc[p, 'max'] = val\n",
    "            elif prop == 'independent':\n",
    "                self.df.loc[p, 'independent'] = val\n",
    "            elif prop == 'fixed':\n",
    "                self.df.loc[p, 'fixed'] = val\n",
    "            else:\n",
    "                raise KeyError(\"Illegal parameter property: {}\".format(prop))\n",
    "\n",
    "    def eval_params(self, params=[], **vals):\n",
    "        \"\"\"Returns parameters for evaluating the function.\n",
    "\n",
    "        Arguments:\n",
    "        params: optional list of values of free parameters\n",
    "        vals: dictionary of parameter values\n",
    "\n",
    "        If a value for a parameter is specified in both `params` and `vals`,\n",
    "        the value from `vals` is used.\n",
    "        Values for independent parameters must be specified in `vals`.\"\"\"\n",
    "        # Add additional values from `params` to vals\n",
    "        if len(params) != 0:\n",
    "            par_names = self.names()\n",
    "            if np.size(par_names) != len(params):\n",
    "                raise ValueError(\"Wrong number of parameters given ({})\".format(len(params)))\n",
    "            for pn, pv in zip(par_names, params):\n",
    "                if pn not in vals:\n",
    "                    vals[pn] = pv\n",
    "\n",
    "        # Fill values unspecified so far from `self.df`\n",
    "        for p in self.df.index.values:\n",
    "            if p not in vals:\n",
    "                if self.df.loc[p, 'independent']:\n",
    "                    raise ValueError(\"Independent parameter `{}` not specified\".format(p))\n",
    "                else:\n",
    "                    vals[p] = self.df.loc[p, 'value']\n",
    "        return vals\n",
    "\n",
    "    def fixed_params(self):\n",
    "        \"\"\"Returns a dictionary of all fixed parameters and their values.\"\"\"\n",
    "        return self.df.loc[self.df.index[self.df.loc[:,'fixed']],'value'].to_dict()\n",
    "\n",
    "    def eval(self, params=[], **vals):\n",
    "        \"\"\"Evaluates the function.\n",
    "\n",
    "        Arguments:\n",
    "        params: optional list of values of free parameters\n",
    "        vals: dictionary of parameter values\n",
    "\n",
    "        If a value for a parameter is specified in both `params` and `vals`,\n",
    "        the value from `vals` is used.\n",
    "        Values for independent parameters must be specified in `vals`.\"\"\"\n",
    "        return self.fun(**self.eval_params(params, **vals))\n",
    "\n",
    "    def freeIdx(self):\n",
    "        \"\"\"Returns a list of names of free parameters\"\"\"\n",
    "        return [p for p in self.df.index.values\n",
    "                if not (self.df.loc[p, 'independent'] or self.df.loc[p, 'fixed'])]\n",
    "\n",
    "    def bounds(self):\n",
    "        \"\"\"Returns a list of bound tuples of free parameters\n",
    "        for use in scipy.optimize.minimize\"\"\"\n",
    "        bnds = []\n",
    "        for p in self.freeIdx():\n",
    "            # Get parameter bounds\n",
    "            min_val = self.df.loc[p, 'min']\n",
    "            max_val = self.df.loc[p, 'max']\n",
    "\n",
    "            # Replace missing values with default minimum and maximum values\n",
    "            if np.isnan(min_val):\n",
    "                min_val = None\n",
    "            if np.isnan(max_val):\n",
    "                max_val = None\n",
    "\n",
    "            # Append to bounds list\n",
    "            bnds.append((min_val, max_val))\n",
    "        return bnds\n",
    "\n",
    "    def initial(self, noise=0):\n",
    "        \"\"\"Returns a numpy.ndarray of initial values for use in scipy.optimize.minimize\n",
    "\n",
    "        Arguments:\n",
    "        :param noise: The amount of random noise to add to the initial values\n",
    "\n",
    "        If ``noise == 0``, no noise is added. If ``noise > 0``, a value in the interval\n",
    "        ``(value - noise, value + noise]`` is added.\"\"\"\n",
    "        vals =  self.df.loc[self.freeIdx(), 'value'].values.copy()\n",
    "        if noise and noise > 0:\n",
    "            vals *= 1 + (np.random.rand(*vals.shape) - 0.5) * (2.0 * noise)\n",
    "        return vals\n",
    "\n",
    "    def index(self, p):\n",
    "        \"\"\"Returns the index of a given parameter in the parameter vector\"\"\"\n",
    "        idx = np.flatnonzero(self.df.index.values == p)\n",
    "        if len(idx) == 0:\n",
    "            raise KeyError(\"Unknown parameter name: {}\".format(p))\n",
    "        return idx[0]\n",
    "\n",
    "    def names(self, onlyFree=True):\n",
    "        \"\"\"Returns an array of the parameter names.\n",
    "\n",
    "        If `onlyFree == True`, only free parameters are returned.\n",
    "        Else, all parameters (including independent and fixed parameters) are returned.\"\"\"\n",
    "        if onlyFree:\n",
    "            return np.array(self.freeIdx(), dtype=np.object_)\n",
    "        else:\n",
    "            return self.df.index.values.copy()\n",
    "\n",
    "    def free_values(self, values):\n",
    "        \"\"\"Returns a dictionary of all free parameters and their values.\n",
    "\n",
    "        Arguments:\n",
    "            values -- an iterable with the values of the free parameters\n",
    "\n",
    "        The entries in `values` must have the same order as the parameters\n",
    "        returned by `FitParameters.names`.\n",
    "        \"\"\"\n",
    "        return {p: v for p, v in zip(self.freeIdx(), values)}\n",
    "\n",
    "    def copy(self):\n",
    "        \"\"\"Returns a deep copy of this instance\"\"\"\n",
    "        return deepcopy(self)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separate models for scipy.optimize.minimize\n",
    "red_p = FitParameters(red, independent='t')\n",
    "red_p.set('tr', min=0, max=MAX_tr, value=tr_0)\n",
    "red_p.set('m_ktl', min=MIN_m_ktl, max=MAX_m_ktl, value=m_ktl_0)\n",
    "red_p.set('kmr', min=0, max=MAX_kmr, value=kmr_0)\n",
    "red_p.set('betr', min=0.001, max=MAX_betr, value=betr_0)\n",
    "red_p.set('deltr', min=0.001, max=MAX_deltr, value=deltr_0)\n",
    "red_p.set('offr', value=offr_0)\n",
    "\n",
    "green_p = FitParameters(green, independent='t')\n",
    "green_p.set('tg', min=0, max=MAX_tg, value=tg_0)\n",
    "green_p.set('m_ktl', min=MIN_m_ktl, max=MAX_m_ktl, value=m_ktl_0)\n",
    "green_p.set('kmg', min=0, max=MAX_kmg, value=kmg_0)\n",
    "green_p.set('betg', min=0.001, max=MAX_betg, value=betg_0)\n",
    "green_p.set('deltg', min=0.001, max=MAX_deltg, value=deltg_0)\n",
    "green_p.set('offg', value=offg_0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Combined model for scipy.optimize.minimize\n",
    "combined_p = FitParameters(combined, independent = 't')\n",
    "combined_p.set('tr', min=0, max=MAX_tr, value=tr_0)\n",
    "combined_p.set('tg', min=0, max=MAX_tg, value=tg_0)\n",
    "combined_p.set('m_ktl', min=MIN_m_ktl, max=MAX_m_ktl, value=m_ktl_0)\n",
    "combined_p.set('kmr', min=0, max=MAX_kmr, value=kmr_0)\n",
    "combined_p.set('kmg', min=0, max=MAX_kmg, value=kmg_0)\n",
    "combined_p.set('betr', min=0, max=MAX_betr, value=betr_0)\n",
    "combined_p.set('betg', min=0, max=MAX_betg, value=betg_0)\n",
    "combined_p.set('deltr', min=0, max=MAX_deltr, value=deltr_0)\n",
    "combined_p.set('deltg', min=0, max=MAX_deltg, value=deltg_0)\n",
    "combined_p.set('offr', value=offr_0)\n",
    "combined_p.set('offg', value=offg_0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Jacobian\n",
    "To increase the efficiency of fitting, the Jacobian matrix of the objective function is provided to the optimization routine.\n",
    "If the objective function is a typical negative log-likelihood function with normal distribution of residuals\n",
    "$$\n",
    " L(\\theta) = \\sum_{t\\in T} \\frac{1}{2\\sigma_t^2} \\big(D_t - f(t\\mid\\theta)\\big)^2 \\text{,}\n",
    "$$\n",
    "where $D_t$ is the measured data at time $t$ and $f(t\\mid\\theta)$ is the value of the model function at time $t$ with parameters $\\theta$, the Jacobian is:\n",
    "$$\\begin{align}\n",
    "\\nabla L(\\theta) &= \\nabla \\sum_{t\\in T} \\frac{1}{2\\sigma_t^2} \\big(D_t - f\\left(t\\,\\middle|\\,\\theta\\right)\\big)^2 \\\\\n",
    "&= \\sum_{t\\in T} \\nabla \\frac{1}{2\\sigma_t^2} \\big( D_t - f\\left(t\\,\\middle|\\,\\theta\\right) \\big)^2 \\\\\n",
    "&= \\sum_{t\\in T} \\frac{2}{2\\sigma_t^2} \\big( D_t - f\\left(t\\,\\middle|\\,\\theta\\right) \\big) \\nabla\\big( D_t - f\\left(t\\,\\middle|\\,\\theta\\right) \\big) \\\\\n",
    "&= \\sum_{t\\in T} \\frac{1}{\\sigma_t^2} \\big( D_t - f\\left(t\\,\\middle|\\,\\theta\\right) \\big)\\big(\\nabla D_t - \\nabla  f\\left(t\\,\\middle|\\,\\theta\\right)\\big) \\\\\n",
    "&= -\\sum_{t\\in T} \\frac{1}{\\sigma_t^2} \\big(D_t - f\\left(t\\,\\middle|\\,\\theta\\right)\\big) \\nabla f\\left(t\\,\\middle|\\,\\theta\\right) \\\\\n",
    "\\end{align}$$\n",
    "We see that for calculating the Jacobian of the objective function we need the Jacobian of the model function.\n",
    "\n",
    "We use the general expression model function\n",
    "$$\n",
    "f\\left(t \\,\\middle|\\, t_0, m, k, \\beta, \\delta, a\\right) = a + m \\left(\\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)} + \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{k + \\beta - \\delta} - \\frac{\\mathrm{e}^{- \\beta \\left(t - t_{0}\\right)}}{\\beta - \\delta}\\right) \\text{,}\n",
    "$$\n",
    "where $t_0$ is the mRNA expression onset time, $m$ is the product of initial mRNA amount and translation rate, $k$ is the maturation rate, $\\beta$ is the protein degradation rate, $\\delta$ is the mRNA degradation rate, and $a$ is a vertical offset."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Jacobian $\\nabla f\\left(t \\,\\middle|\\, t_0, m, k, \\beta, \\delta, a\\right)$ of the general expression model function is the vector of the derivatives with respect to the various parameters:\n",
    "$$\\begin{align*}\n",
    "\\frac{\\partial f}{\\partial t_0} &= m \\left(\\frac{k \\delta \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)} - \\frac{\\beta \\mathrm{e}^{- \\beta \\left(t - t_{0}\\right)}}{\\beta - \\delta} + \\frac{\\left(k + \\beta\\right) \\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{k + \\beta - \\delta}\\right)\\\\\n",
    "\\frac{\\partial f}{\\partial m} &= \\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)} + \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{k + \\beta - \\delta} - \\frac{\\mathrm{e}^{- \\beta \\left(t - t_{0}\\right)}}{\\beta - \\delta}\\\\\n",
    "\\frac{\\partial f}{\\partial k} &= m \\left(- \\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)^{2}} + \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{k + \\beta - \\delta} \\left(- t + t_{0}\\right) - \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{\\left(k + \\beta - \\delta\\right)^{2}} + \\frac{\\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)}\\right)\\\\\n",
    "\\frac{\\partial f}{\\partial \\beta} &= m \\left(- \\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)^{2}} - \\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right)^{2} \\left(k + \\beta - \\delta\\right)} + \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{k + \\beta - \\delta} \\left(- t + t_{0}\\right) - \\frac{\\mathrm{e}^{- \\beta \\left(t - t_{0}\\right)}}{\\beta - \\delta} \\left(- t + t_{0}\\right) - \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{\\left(k + \\beta - \\delta\\right)^{2}} + \\frac{\\mathrm{e}^{- \\beta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right)^{2}}\\right)\\\\\n",
    "\\frac{\\partial f}{\\partial \\delta} &= m \\left(\\frac{k \\left(- t + t_{0}\\right) \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)} + \\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right) \\left(k + \\beta - \\delta\\right)^{2}} + \\frac{k \\mathrm{e}^{- \\delta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right)^{2} \\left(k + \\beta - \\delta\\right)} + \\frac{\\mathrm{e}^{\\left(- k - \\beta\\right) \\left(t - t_{0}\\right)}}{\\left(k + \\beta - \\delta\\right)^{2}} - \\frac{\\mathrm{e}^{- \\beta \\left(t - t_{0}\\right)}}{\\left(\\beta - \\delta\\right)^{2}}\\right)\\\\\n",
    "\\frac{\\partial f}{\\partial a} &= 1\\\\\n",
    "\\end{align*}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def general_jacobian(t, t0, m, k, beta, delta, a):\n",
    "    \"\"\"Returns the Jacobi matrix of the general expression model function\n",
    "    with time along axis=0 and parameters along axis=1\"\"\"\n",
    "\n",
    "    # Initialize Jacobian\n",
    "    jac = np.zeros((np.size(t), 6))\n",
    "    #jac = np.zeros((np.size(t), 4))\n",
    "\n",
    "    # Get time after onset “kink”\n",
    "    after_t0 = (t > t0)\n",
    "\n",
    "    # Define abbreviations for frequent terms\n",
    "    dt = t[after_t0] - t0\n",
    "    bmd = beta - delta\n",
    "    kbd = k + bmd\n",
    "    kpb = k + beta\n",
    "\n",
    "    # Derive w.r.t. t0\n",
    "    jac[after_t0, 0] = m * (k * delta * np.exp(-delta * dt) / bmd / kbd \n",
    "                          - beta * np.exp(-beta * dt) / bmd\n",
    "                          + kpb * np.exp(-kpb * dt) / kbd)\n",
    "\n",
    "    # Derive w.r.t. m\n",
    "    jac[after_t0, 1] = (k * np.exp(-delta * dt) / bmd / kbd\n",
    "                      + np.exp(-kpb * dt) / kbd\n",
    "                      - np.exp(-beta * dt) / bmd)\n",
    "\n",
    "    # Derive w.r.t. k\n",
    "    jac[after_t0, 2] = m * (-k * np.exp(-delta * dt) / bmd / kbd**2\n",
    "                           - dt * np.exp(-kpb * dt) / kbd\n",
    "                           - np.exp(-kpb * dt) / kbd**2\n",
    "                           + np.exp(-delta * dt) / bmd / kbd)\n",
    "\n",
    "    # Derive w.r.t. beta\n",
    "    jac[after_t0, 3] = m * (-k * np.exp(-delta * dt) / bmd / kbd**2\n",
    "                           - k * np.exp(-delta * dt) / bmd**2 / kbd\n",
    "                           - dt * np.exp(-kpb * dt) / kbd\n",
    "                           + dt * np.exp(-beta * dt) / bmd\n",
    "                           - np.exp(-kpb * dt) / kbd**2\n",
    "                           + np.exp(-beta * dt) / bmd**2)\n",
    "\n",
    "    # Derive w.r.t. delta\n",
    "    jac[after_t0, 4] = m * (-k * dt * np.exp(-delta * dt) / bmd / kbd\n",
    "                           + k * np.exp(-delta * dt) / bmd / kbd**2\n",
    "                           + k * np.exp(-delta * dt) / bmd**2 / kbd\n",
    "                           + np.exp(-kpb * dt) / kbd**2\n",
    "                           - np.exp(-beta * dt) / bmd**2)\n",
    "\n",
    "    # Derive w.r.t. a\n",
    "    jac[:, 5] = 1\n",
    "\n",
    "    return jac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def red_jacobian(t, tr, m_ktl, kmr, betr, deltr, offr):\n",
    "    \"\"\"Wrapper function for Jacobian of red model function\"\"\"\n",
    "    return general_jacobian(t=t, t0=tr, m=m_ktl, k=kmr, beta=betr, delta=deltr, a=offr)\n",
    "\n",
    "def green_jacobian(t, tg, m_ktl, kmg, betg, deltg, offg):\n",
    "    \"\"\"Wrapper function for Jacobian of red model function\"\"\"\n",
    "    return general_jacobian(t=t, t0=tg, m=m_ktl, k=kmg, beta=betg, delta=deltg, a=offg)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in data and prepare result list"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate kernel density estimation of parameter distributions\n",
    "def parameter_KDE(par_tab):\n",
    "    \"\"\"Returns a kernel density estimation of the parameter values for plotting\"\"\"\n",
    "    dens_res = 200\n",
    "    bw_div = 15\n",
    "\n",
    "    par_dist = {}\n",
    "\n",
    "    for par_name in par_tab.columns:\n",
    "        # Get parameter values\n",
    "        par_vals = par_tab.loc[:,par_name].values\n",
    "        par_vals = par_vals.reshape((par_vals.size, 1))\n",
    "\n",
    "        # Test parameter values for validity\n",
    "        if np.any(np.logical_not(np.isfinite(par_vals))):\n",
    "            print(\"Warning: invalid values encountered for “{}”\".format(par_name))\n",
    "            par_vals = par_vals[np.isfinite(par_vals)]\n",
    "            if par_vals.size > 0:\n",
    "                # Calculate distribution of valid entries\n",
    "                par_vals = par_vals.reshape((par_vals.size, 1))\n",
    "            else:\n",
    "                # No valid entries found; cancel distribution calculation\n",
    "                par_dist[par_name] = {'val': [], 'prob': []}\n",
    "                continue\n",
    "\n",
    "        # Get parameter extrema and bandwidth\n",
    "        par_min = np.min(par_vals)\n",
    "        par_max = np.max(par_vals)\n",
    "        bw = (par_max - par_min) / bw_div\n",
    "\n",
    "        # Avoid error if only one parameter value found\n",
    "        if bw == 0.:\n",
    "            par_dist[par_name] = {'val': [], 'prob': []}\n",
    "            continue\n",
    "\n",
    "        # Get kernel density estimation of parameter values\n",
    "        kde = KernelDensity(kernel='epanechnikov', bandwidth=bw).fit(par_vals)\n",
    "        par_x = np.linspace(par_min, par_max, dens_res).reshape((dens_res, 1))\n",
    "        par_dens = np.exp(kde.score_samples(par_x))\n",
    "\n",
    "        # Adjust values for nicer plotting (KDE >= 0, edges == 0)\n",
    "        #par_dens[par_dens < 0] = 0\n",
    "        if par_dens[0] != 0:\n",
    "            par_dens = np.insert(par_dens, 0, 0)\n",
    "            par_x = np.insert(par_x, 0, par_min)\n",
    "        if par_dens[-1] != 0:\n",
    "            par_dens = np.append(par_dens, 0)\n",
    "            par_x = np.append(par_x, par_max)\n",
    "\n",
    "        # Insert KDE into dict\n",
    "        par_dist[par_name] = {'val': par_x.flatten(), 'prob': par_dens.flatten()}\n",
    "    return par_dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kde(ax, dist, label, clr_face='b', clr_edge='k', mark=None):\n",
    "    \"\"\"Plots the current parameter value in relation to the distribution\n",
    "    in the whole dataset.\"\"\"\n",
    "    ax.fill_betweenx(dist['val'], dist['prob'], color=clr_face)\n",
    "    if mark != None:\n",
    "        ax.axhline(y=mark, color=clr_edge)\n",
    "    ax.set_xticks([])\n",
    "    ax.spines['left'].set_position('zero')\n",
    "    for s in [ax.spines[pos] for pos in ['bottom', 'right', 'top']]:\n",
    "        s.set_visible(False)\n",
    "    ax.set_title(label)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def getDataLabel(d, filename=False):\n",
    "    \"\"\"Returns a nicely formatted name for the data info dict `d`.\n",
    "    `d` must have the keys \"measurement\", \"sample\" and \"condition\",\n",
    "    such as the elements of `D`.\n",
    "    Set `filename=True` for a filename-friendly output.\"\"\"\n",
    "    if filename:\n",
    "        return \"{0[measurement]}_{0[sample]}_{0[condition]}\".format(d)\n",
    "    return \"{0[sample]}: {0[condition]} [{0[measurement]}]\".format(d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def average_data():\n",
    "    \"\"\"Create average dataset\"\"\"\n",
    "    # Define helper functions\n",
    "    def _add_to_Ddict(sample, condition, meas, seq, i):\n",
    "        \"\"\"Helper function: add average to Ddict\"\"\"\n",
    "        nonlocal Ddict\n",
    "        if sample not in Ddict:\n",
    "            Ddict[sample] = {}\n",
    "        if condition not in Ddict[sample]:\n",
    "            Ddict[sample][condition] = {}\n",
    "        if meas not in Ddict[sample][condition]:\n",
    "            Ddict[sample][condition][meas] = {}\n",
    "        Ddict[sample][condition][meas][seq] = i\n",
    "\n",
    "    def _check_t_continuity(idx):\n",
    "        '''DEBUG: Ensure that time vectors do not \"intercalate\"'''\n",
    "        has_true_started = False\n",
    "        has_true_finished = False\n",
    "        for x in idx:\n",
    "            if x:\n",
    "                if not has_true_started:\n",
    "                    has_true_started = True\n",
    "                elif has_true_finished:\n",
    "                    print(\"Intercalating time vector detected\", file=sys.stderr)\n",
    "                    import pdb\n",
    "                    pdb.set_trace()\n",
    "            elif has_true_started:\n",
    "                has_true_finished = True\n",
    "\n",
    "    def _get_average_data(X):\n",
    "        \"\"\"Average over data in X\"\"\"\n",
    "        new_nTraces = 0\n",
    "        new_t = np.unique([t for x in X for t in x['t']])\n",
    "        t_weight = np.zeros((new_t.size, 1))\n",
    "        new_rfp = np.zeros((new_t.size, 1))\n",
    "        new_gfp = np.zeros((new_t.size, 1))\n",
    "        new_rfp_error = np.zeros((new_t.size, 1))\n",
    "        new_gfp_error = np.zeros((new_t.size, 1))\n",
    "        for x in X:\n",
    "            idx = np.isin(new_t, x['t'])\n",
    "            _check_t_continuity(idx) #DEBUG\n",
    "            n = x['nTraces']\n",
    "            new_nTraces += n\n",
    "            try:\n",
    "                weight_x = x['t_weight']\n",
    "            except KeyError:\n",
    "                weight_x = n\n",
    "            t_weight[idx, :] += weight_x\n",
    "            new_rfp[idx, :] += x['rfp'] * weight_x\n",
    "            new_gfp[idx, :] += x['gfp'] * weight_x\n",
    "            new_rfp_error[idx, :] += x['rfp_error'] * weight_x\n",
    "            new_gfp_error[idx, :] += x['gfp_error'] * weight_x\n",
    "        new_rfp /= t_weight\n",
    "        new_gfp /= t_weight\n",
    "        new_rfp_error /= t_weight\n",
    "        new_gfp_error /= t_weight\n",
    "\n",
    "        return {\n",
    "            'nTraces': new_nTraces,\n",
    "            't': new_t,\n",
    "            't_weight': t_weight,\n",
    "            'rfp': new_rfp,\n",
    "            'gfp': new_gfp,\n",
    "            'rfp_error': new_rfp_error,\n",
    "            'gfp_error': new_gfp_error,\n",
    "        }\n",
    "\n",
    "    # Obtain grouping information\n",
    "    global D\n",
    "    Ddict = {}\n",
    "    for i, d in enumerate(D):\n",
    "        if 'test' in d['condition']:\n",
    "            continue\n",
    "        try:\n",
    "            meas, seq = d['measurement'].split('_')\n",
    "        except TypeError:\n",
    "            continue\n",
    "        _add_to_Ddict(d['sample'], d['condition'], meas, seq, i)\n",
    "\n",
    "    # Create new list of averaged data\n",
    "    Dnew = []\n",
    "    for sample, sample_dict in Ddict.items():\n",
    "        for condition, condition_dict in sample_dict.items():\n",
    "            l_all = []\n",
    "            for meas, meas_dict in condition_dict.items():\n",
    "                l_meas = []\n",
    "                for seq, i in meas_dict.items():\n",
    "                    # Get averages per experiment/sequence\n",
    "                    d_seq = {\n",
    "                        'sample': \"{}_average\".format(sample),\n",
    "                        'condition': condition,\n",
    "                        'measurement': \"{}_{}\".format(meas, seq),\n",
    "                        'file': \"\",\n",
    "                        't': D[i]['t'],\n",
    "                        'rfp': np.average(D[i]['rfp'], axis=1).reshape((-1, 1)),\n",
    "                        'gfp': np.average(D[i]['gfp'], axis=1).reshape((-1, 1)),\n",
    "                        'rfp_error': np.average(D[i]['rfp_error'], axis=1).reshape((-1, 1)),\n",
    "                        'gfp_error': np.average(D[i]['gfp_error'], axis=1).reshape((-1, 1)),\n",
    "                        'nTraces': D[i]['nTraces'],\n",
    "                    }\n",
    "                    Dnew.append(d_seq)\n",
    "\n",
    "                    l_meas.append({\n",
    "                        'nTraces': d_seq['nTraces'],\n",
    "                        't': d_seq['t'],\n",
    "                        'rfp': d_seq['rfp'],\n",
    "                        'gfp': d_seq['gfp'],\n",
    "                        'rfp_error': d_seq['rfp_error'],\n",
    "                        'gfp_error': d_seq['gfp_error'],\n",
    "                    })\n",
    "\n",
    "                # Get averages per condition and day\n",
    "                d_meas = _get_average_data(l_meas)\n",
    "                d_meas['sample'] = \"{}_average\".format(sample)\n",
    "                d_meas['condition'] = condition\n",
    "                d_meas['measurement'] = \"{}_all\".format(meas)\n",
    "                d_meas['file'] = \"\"\n",
    "                Dnew.append(d_meas)\n",
    "\n",
    "                l_all.append({\n",
    "                    'nTraces': d_meas['nTraces'],\n",
    "                    't': d_meas['t'],\n",
    "                    't_weight': d_meas['t_weight'],\n",
    "                    'rfp': d_meas['rfp'],\n",
    "                    'gfp': d_meas['gfp'],\n",
    "                    'rfp_error': d_meas['rfp_error'],\n",
    "                    'gfp_error': d_meas['gfp_error'],\n",
    "                })\n",
    "\n",
    "            # Get averages per condition\n",
    "            d_all = _get_average_data(l_all)\n",
    "            d_all['sample'] = \"{}_average\".format(sample)\n",
    "            d_all['condition'] = condition\n",
    "            d_all['measurement'] = \"all\"\n",
    "            d_all['file'] = \"\"\n",
    "            Dnew.append(d_all)\n",
    "    return Dnew"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Prepare data loading\n",
    "\n",
    "# Define available files\n",
    "datafiles = [\n",
    "    {\n",
    "    #    \"sample\": \"A549\",\n",
    "    #    \"condition\": \"control\",\n",
    "    #    \"measurement\": \"Test\",\n",
    "    #    \"file\": \"data/A549_control_test.xlsx\"\n",
    "    #}, {\n",
    "        \"sample\": \"A549\",\n",
    "        \"condition\": \"siRNA\",\n",
    "        \"measurement\": \"2016-01-09_seq3\",\n",
    "        \"file\": \"data/2016-01-09_seq3_A549_siRNA_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"A549\",\n",
    "        \"condition\": \"control\",\n",
    "        \"measurement\": \"2016-01-09_seq5\",\n",
    "        \"file\": \"data/2016-01-09_seq5_A549_Control_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"A549\",\n",
    "        \"condition\": \"siRNA\",\n",
    "        \"measurement\": \"2016-12-20_seq3\",\n",
    "        \"file\": \"data/2016-12-20_seq3_A549_siRNA_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"A549\",\n",
    "        \"condition\": \"control\",\n",
    "        \"measurement\": \"2016-12-20_seq4\",\n",
    "        \"file\": \"data/2016-12-20_seq4_A549_control_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"siRNA\",\n",
    "        \"measurement\": \"2017-05-26_seq10\",\n",
    "        \"file\": \"data/2017-05-26_seq10_Huh7_siRNA_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"siRNA\",\n",
    "        \"measurement\": \"2017-05-26_seq11\",\n",
    "        \"file\": \"data/2017-05-26_seq11_Huh7_siRNA_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"siRNA\",\n",
    "        \"measurement\": \"2017-06-02_seq4\",\n",
    "        \"file\": \"data/2017-06-02_seq4_Huh7_siRNA_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"siRNA\",\n",
    "        \"measurement\": \"2017-06-02_seq5\",\n",
    "        \"file\": \"data/2017-06-02_seq5_Huh7_siRNA_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"control\",\n",
    "        \"measurement\": \"2017-05-26_seq6\",\n",
    "        \"file\": \"data/2017-05-26_seq6_Huh7_control_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"control\",\n",
    "        \"measurement\": \"2017-05-26_seq7\",\n",
    "        \"file\": \"data/2017-05-26_seq7_Huh7_control_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"control\",\n",
    "        \"measurement\": \"2017-06-02_seq6\",\n",
    "        \"file\": \"data/2017-06-02_seq6_Huh7_control_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"control\",\n",
    "        \"measurement\": \"2017-06-02_seq7\",\n",
    "        \"file\": \"data/2017-06-02_seq7_Huh7_control_#molecules.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"test\",\n",
    "        \"measurement\": \"2017-06-02_seq6,5\",\n",
    "        \"file\": \"data/2017-06-02_seq6,5_Huh7_test30_#molecules_exemplary13.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"test100control\",\n",
    "        \"measurement\": \"2017-06-02_seq6\",\n",
    "        \"file\": \"data/2017-06-02_seq6_Huh7_test100_#molecules_exemplary13.xlsx\"\n",
    "    }, {\n",
    "        \"sample\": \"Huh7\",\n",
    "        \"condition\": \"test100silenced\",\n",
    "        \"measurement\": \"2017-06-02_seq5\",\n",
    "        \"file\": \"data/2017-06-02_seq5_Huh7_test100_#molecules_exemplary13.xlsx\"\n",
    "    }\n",
    "]\n",
    "\n",
    "# By default, mark all files for loading\n",
    "load_idcs = range(len(datafiles))\n",
    "\n",
    "# Define function for loading data\n",
    "def load_data_from_files():\n",
    "    \"\"\"Loads data from specified files into `D`.\n",
    "    Requires `load_idcs` to hold a list of indices to `datafiles`.\"\"\"\n",
    "    global D\n",
    "    D = []\n",
    "    for i in load_idcs:\n",
    "        # Show message\n",
    "        print(\"Loading file: {}\".format(datafiles[i][\"file\"]))\n",
    "\n",
    "        # Read sheets from excel file\n",
    "        X = pd.read_excel(datafiles[i]['file'], dtype=np.float64, sheet_name=[\n",
    "            '#RFP', '#GFP_corrected', '#RFP_error', '#GFP_error'])\n",
    "\n",
    "        # Write data into easy-to-access structure\n",
    "        d = {}\n",
    "        d['sample'] = datafiles[i]['sample']\n",
    "        d['condition'] = datafiles[i]['condition']\n",
    "        d['measurement'] = datafiles[i]['measurement']\n",
    "        d['file'] = datafiles[i]['file']\n",
    "        d['t'] = X['#RFP'].values[:,0].flatten()\n",
    "        #d['rfp'] = X['RFP'].values[:,1:]\n",
    "        #d['gfp'] = X['GFP_corrected'].values[:,1:]\n",
    "        d['rfp'] = X['#RFP'].values[:,1:]\n",
    "        d['gfp'] = X['#GFP_corrected'].values[:,1:]\n",
    "        d['rfp_error'] = X['#RFP_error'].values[:,1:]\n",
    "        d['gfp_error'] = X['#GFP_error'].values[:,1:]\n",
    "        d['nTraces'] = np.shape(d['rfp'])[1]\n",
    "        d['iFile'] = i\n",
    "        D.append(d)\n",
    "    print(\"Done.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read in data from excel sheets\n",
    "\n",
    "# Prompt user for files to load\n",
    "lbl = wdg.Label('Select the files to load:')\n",
    "lbl.layout.width = 'initial'\n",
    "entries = []\n",
    "for f in datafiles:\n",
    "    entries.append(\"{} {}: {}\".format(\n",
    "        f['sample'], f['condition'], f['file']))\n",
    "sel_entry = wdg.SelectMultiple(options=entries, rows=len(entries))\n",
    "sel_entry.layout.width = 'initial'\n",
    "bload = wdg.Button(description='Load')\n",
    "bselall = wdg.Button(description='Select all')\n",
    "bselnone = wdg.Button(description='Select none')\n",
    "\n",
    "# Define callbacks\n",
    "def sel_all_files(_):\n",
    "    sel_entry.value = entries\n",
    "def sel_no_files(_):\n",
    "    sel_entry.value = ()\n",
    "def load_button_clicked(_):\n",
    "    global load_idcs\n",
    "    load_idcs = [entries.index(r) for r in sel_entry.value]\n",
    "    vb.close()\n",
    "    load_data_from_files()\n",
    "bselall.on_click(sel_all_files)\n",
    "bselnone.on_click(sel_no_files)\n",
    "bload.on_click(load_button_clicked)\n",
    "\n",
    "# Finally, show the widgets\n",
    "vb = wdg.VBox((lbl, sel_entry, wdg.HBox((bload,bselall,bselnone))))\n",
    "IPython.display.display(vb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# (optional) Calculate average of data and discard original data\n",
    "D = average_data()"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Insert here the timestamp you want to enforce for the parameters.\n",
    "\n",
    "# `require_timestamp` is the enforced timestamp as string or\n",
    "# `None` for loading the latest file.\n",
    "# Example:\n",
    "# to enforce timestamp \"2012-03-04–123456\",\n",
    "# execute: require_timestamp = \"2012-03-04–123456\"\n",
    "require_timestamp = None"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Assess data files to load\n",
    "outpath = getOutpath()\n",
    "file_re = re.compile(r\"(?P<date>\\d{4}-\\d{2}-\\d{2}–\\d{6})_fixed_distribution_moments\\.xlsx$\")\n",
    "\n",
    "# Search file for fixed parameters\n",
    "# After this cell, `fixparfile` will be a string of the path to the file\n",
    "# holding the fixed parameter values, or `None` if no file was found.\n",
    "fixparfile = None\n",
    "ts = None\n",
    "files_found = {}\n",
    "\n",
    "# Iteratively search parameter files and add them to dictionary\n",
    "for f in os.listdir(outpath):\n",
    "    m = file_re.match(f)\n",
    "    if not m:\n",
    "        continue\n",
    "    files_found[m.group('date')] = os.path.join(outpath, m.string)\n",
    "\n",
    "# Select newest (or preferred) timestamp or set `fixparfile` to None,\n",
    "# if no file was found\n",
    "if len(files_found) > 0:\n",
    "    if require_timestamp:\n",
    "        if require_timestamp in files_found.keys():\n",
    "            # Make `files_found` point to the requested file\n",
    "            fixparfile = files_found[require_timestamp]\n",
    "            ts = require_timestamp\n",
    "        else:\n",
    "            print(\"Requested timestamp “{}” not found.\".format(require_timestamp))\n",
    "    else:\n",
    "        ts = sorted(files_found.keys(), reverse=True)[0]\n",
    "        fixparfile = files_found[ts]\n",
    "\n",
    "if fixparfile is not None:\n",
    "    print(\"Found parameter data with timestamp “{}”:\\n{}\".format(ts, fixparfile))\n",
    "else:\n",
    "    print(\"No parameter data found.\")\n",
    "\n",
    "# DEBUG\n",
    "#for d, c in datafiles.items():\n",
    "#    print(d + \": \" + str(c))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Load fixed parameters\n",
    "F = pd.read_excel(fixparfile, dtype=np.float64, sheet_name=['red', 'green'])\n",
    "\n",
    "red_p.set('kmr', fixed=True, value=F['red'].loc['kmr', 'mean'])\n",
    "red_p.set('betr', fixed=True, value=F['red'].loc['betr', 'mean'])\n",
    "\n",
    "green_p.set('kmg', fixed=True, value=F['green'].loc['kmg', 'mean'])\n",
    "green_p.set('betg', fixed=True, value=F['green'].loc['betg', 'mean'])\n",
    "\n",
    "combined_p.set('kmr', fixed=True, value=F['red'].loc['kmr', 'mean'])\n",
    "combined_p.set('betr', fixed=True, value=F['red'].loc['betr', 'mean'])\n",
    "combined_p.set('kmg', fixed=True, value=F['green'].loc['kmg', 'mean'])\n",
    "combined_p.set('betg', fixed=True, value=F['green'].loc['betg', 'mean'])"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Set no fixed parameters\n",
    "red_p.set('kmr', fixed=False)\n",
    "red_p.set('betr', fixed=False)\n",
    "\n",
    "green_p.set('kmg', fixed=False)\n",
    "green_p.set('betg', fixed=False)\n",
    "\n",
    "combined_p.set('kmr', fixed=False)\n",
    "combined_p.set('betr', fixed=False)\n",
    "combined_p.set('kmg', fixed=False)\n",
    "combined_p.set('betg', fixed=False)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Optional: enforce shorter times\n",
    "t_max = 15\n",
    "for d in D:\n",
    "    idx = d['t'] <= t_max\n",
    "    d['t'] = d['t'][idx]\n",
    "    d['rfp'] = d['rfp'][idx,:]\n",
    "    d['gfp'] = d['gfp'][idx,:]\n",
    "    d['rfp_error'] = d['rfp_error'][idx,:]\n",
    "    d['gfp_error'] = d['gfp_error'][idx,:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Provide output tables\n",
    "\n",
    "# Initialize result dictionary\n",
    "R = []\n",
    "\n",
    "# Get a list of fit parameters\n",
    "par_names = green_p.names().tolist()\n",
    "par_names.extend(p for p in red_p.names() if p not in par_names)\n",
    "par_names.sort()\n",
    "\n",
    "# Iteratively populate the result dictionary\n",
    "for k in range(len(D)):\n",
    "    R.insert(k, {})\n",
    "    nTraces = np.shape(D[k]['gfp'])[1]\n",
    "    nTimes = np.shape(D[k]['gfp'])[0]\n",
    "    tpl_traces = np.empty((nTimes, nTraces))\n",
    "    tpl_traces.fill(np.NaN)\n",
    "\n",
    "    R[k]['green'] = {}\n",
    "    R[k]['green']['params'] = pd.DataFrame(index=np.arange(nTraces), columns=green_p.names(), dtype='float64')\n",
    "    R[k]['green']['fit'] = np.copy(tpl_traces)\n",
    "    R[k]['green']['success'] = np.zeros(nTraces, dtype=np.bool_)\n",
    "    R[k]['green']['chisq'] = np.full(nTraces, np.NaN)\n",
    "\n",
    "    R[k]['red'] = {}\n",
    "    R[k]['red']['params'] = pd.DataFrame(index=np.arange(nTraces), columns=red_p.names(), dtype='float64')\n",
    "    R[k]['red']['fit'] = np.copy(tpl_traces)\n",
    "    R[k]['red']['success'] = np.zeros(nTraces, dtype=np.bool_)\n",
    "    R[k]['red']['chisq'] = np.full(nTraces, np.NaN)\n",
    "\n",
    "    R[k]['combined'] = {}\n",
    "    R[k]['combined']['params'] = pd.DataFrame(index=np.arange(nTraces), columns=combined_p.names(), dtype='float64')\n",
    "    #R[k]['combined']['fit'] = {}\n",
    "    R[k]['combined']['success'] = np.zeros(nTraces, dtype=np.bool_)\n",
    "    R[k]['combined']['chisq'] = np.full((nTraces, 2), np.NaN)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pickle or load fitting results\n",
    "Pickling is only reasonable if the result list `R` has already been populated by fitting (see below)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pickle fit results for future sessions\n",
    "outfile = getOutpath('fit_results.pickled')\n",
    "with open(outfile, 'wb') as f:\n",
    "    pickle.dump(R, f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Load pickled results (requires file suffix “.pickled”)\n",
    "pickfiles = [f for f in os.listdir(getOutpath()) if f.lower().endswith('.pickled')]\n",
    "pickfiles.sort(reverse=True)\n",
    "\n",
    "lbl = wdg.Label('Select the file to load:')\n",
    "lbl.layout.width = 'initial'\n",
    "rad = wdg.RadioButtons(options=pickfiles)\n",
    "but = wdg.Button(description='Load')\n",
    "vb = wdg.VBox([lbl, rad, but])\n",
    "IPython.display.display(vb)\n",
    "\n",
    "def clicked_on_but(b):\n",
    "    global R\n",
    "    with open(getOutpath(rad.value, ''), 'rb') as f:\n",
    "        R = pickle.load(f)\n",
    "    print('Loaded: ' + rad.value)\n",
    "    vb.close()\n",
    "but.on_click(clicked_on_but)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write results to XLSX\n",
    "if len(R) != len(D):\n",
    "    raise ValueError(\"R and D must have the same length!\")\n",
    "\n",
    "for k in range(len(D)):\n",
    "    # Collect information about this file\n",
    "    measurement = D[k]['measurement']\n",
    "    sample = D[k]['sample']\n",
    "    condition = D[k]['condition']\n",
    "    time = D[k]['t']\n",
    "    rfp_raw = D[k]['rfp']\n",
    "    gfp_raw = D[k]['gfp']\n",
    "    rfp_error = D[k]['rfp_error']\n",
    "    gfp_error = D[k]['gfp_error']\n",
    "\n",
    "    rfp_fit = R[k]['red']['fit']\n",
    "    gfp_fit = R[k]['green']['fit']\n",
    "    rfp_params = R[k]['red']['params']\n",
    "    gfp_params = R[k]['green']['params']\n",
    "    rfp_chsq = R[k]['red']['chisq']\n",
    "    gfp_chsq = R[k]['green']['chisq']\n",
    "\n",
    "    try:\n",
    "        hasCombined = True\n",
    "        cmb_fit_rfp = R[k]['combined']['fit']['red']\n",
    "        cmb_fit_gfp = R[k]['combined']['fit']['green']\n",
    "        cmb_params = R[k]['combined']['params']\n",
    "        cmb_chsq = R[k]['combined']['chisq']\n",
    "    except:\n",
    "        hasCombined = False\n",
    "\n",
    "    # Write data to file\n",
    "    file = getOutpath(\"{}_{}_{}.xlsx\".format(sample, measurement, condition))\n",
    "    xlsx_writer = pd.ExcelWriter(file, engine='xlsxwriter')\n",
    "\n",
    "    pd.DataFrame(time).to_excel(xlsx_writer, sheet_name=\"t\")\n",
    "    pd.DataFrame(rfp_raw).to_excel(xlsx_writer, sheet_name=\"RFP_raw\")\n",
    "    pd.DataFrame(gfp_raw).to_excel(xlsx_writer, sheet_name=\"GFP_raw\")\n",
    "    pd.DataFrame(rfp_error).to_excel(xlsx_writer, sheet_name=\"RFP_error\")\n",
    "    pd.DataFrame(gfp_error).to_excel(xlsx_writer, sheet_name=\"GFP_error\")\n",
    "\n",
    "    pd.DataFrame(rfp_fit).to_excel(xlsx_writer, sheet_name=\"RFP_fit\")\n",
    "    pd.DataFrame(gfp_fit).to_excel(xlsx_writer, sheet_name=\"GFP_fit\")\n",
    "    pd.DataFrame(rfp_chsq).to_excel(xlsx_writer, sheet_name=\"RFP_chisq\")\n",
    "    pd.DataFrame(gfp_chsq).to_excel(xlsx_writer, sheet_name=\"GFP_chisq\")\n",
    "    rfp_params.to_excel(xlsx_writer, sheet_name=\"RFP_params\")\n",
    "    gfp_params.to_excel(xlsx_writer, sheet_name=\"GFP_params\")\n",
    "\n",
    "    if hasCombined:\n",
    "        pd.DataFrame(cmb_fit_rfp).to_excel(xlsx_writer, sheet_name=\"RFP_fit_cmb\")\n",
    "        pd.DataFrame(cmb_fit_gfp).to_excel(xlsx_writer, sheet_name=\"GFP_fit_cmb\")\n",
    "        cmb_params.to_excel(xlsx_writer, sheet_name=\"params_cmb\")\n",
    "        pd.DataFrame(cmb_chsq).to_excel(xlsx_writer, sheet_name=\"chisq_cmb\")\n",
    "    \n",
    "    xlsx_writer.save()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit and plot separate models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotSeparate(ds, tr, pdf=None, par_kde=None):\n",
    "    \"\"\"Fits and plots the data, treating RFP and GFP separately.\n",
    "\n",
    "    Keyword arguments:\n",
    "    ds -- the dictionary key of the dataset\n",
    "    tr -- the index of the trace in the dataset to be processed\n",
    "    pdf -- a PdfPages object to which the figure is written if it is not None\n",
    "    par_kde -- if containing dict of values of parameter distributions, plot distributions\n",
    "    \"\"\"\n",
    "    global i_plot\n",
    "\n",
    "    # Plot fit results\n",
    "    fig = plt.figure()\n",
    "\n",
    "    if par_kde != None:\n",
    "        fig.set_figwidth(1.6 * fig.get_figwidth())\n",
    "\n",
    "        pn_red = [p for p in ['m_ktl', 'tr', 'kmr', 'betr', 'deltr', 'offr']\n",
    "                 if not red_p.df.loc[p,'fixed']]\n",
    "        pn_green = [p for p in ['m_ktl', 'tg', 'kmg', 'betg', 'deltg', 'offg']\n",
    "                 if not green_p.df.loc[p,'fixed']]\n",
    "\n",
    "        grid = (2, max(len(pn_red), len(pn_green)))\n",
    "        gs = GridSpec(grid[0], grid[1])\n",
    "\n",
    "        # Plot green parameters\n",
    "        for pi, label in enumerate(pn_green):\n",
    "            ax = plt.subplot(gs.new_subplotspec((0, pi)))\n",
    "            data = par_kde['green'][label]\n",
    "            clr_face = '#00ff0055'\n",
    "            clr_edge = '#009900ff'\n",
    "            curr_val = R[ds]['green']['params'].loc[tr,label]\n",
    "            plot_kde(ax, data, label, clr_face, clr_edge, curr_val)\n",
    "\n",
    "        # Plot red parameters\n",
    "        for pi, label in enumerate(pn_red):\n",
    "            ax = plt.subplot(gs.new_subplotspec((1, pi)))\n",
    "            data = par_kde['red'][label]\n",
    "            clr_face = '#ff000055'\n",
    "            clr_edge = '#990000ff'\n",
    "            curr_val = R[ds]['red']['params'].loc[tr,label]\n",
    "            plot_kde(ax, data, label, clr_face, clr_edge, curr_val)\n",
    "\n",
    "        # Adjust subplot layout\n",
    "        gs.tight_layout(fig, pad=0, rect=(0.5, 0, 1, 1))\n",
    "\n",
    "        # Create axes for fit\n",
    "        gs_fit = GridSpec(1, 1)\n",
    "        ax = fig.add_subplot(gs_fit[0])\n",
    "        gs_fit.tight_layout(fig, pad=0, rect=(0, 0, 0.5, 1))\n",
    "\n",
    "    else:\n",
    "        ax = fig.gca()\n",
    "\n",
    "    p_tr = ax.axvline(R[ds]['red']['params']['tr'][tr], label='RFP onset',\n",
    "                       color='#ff0000', linewidth=.5, linestyle='--')\n",
    "    p_tg = ax.axvline(R[ds]['green']['params']['tg'][tr], label='GFP onset',\n",
    "                      color='#00ff00', linewidth=.5, linestyle='--')\n",
    "    p_fr, = ax.plot(D[ds]['t'], R[ds]['red']['fit'][:,tr], '-', label='RFP (fit)', color='#ff0000', linewidth=1)\n",
    "    p_fg, = ax.plot(D[ds]['t'], R[ds]['green']['fit'][:,tr], '-', label='GFP (fit)', color='#00ff00', linewidth=1)\n",
    "    p_dr, = ax.plot(D[ds]['t'], D[ds]['rfp'][:,tr], '-', label='RFP (measured)', color='#990000', linewidth=.5)\n",
    "    p_dg, = ax.plot(D[ds]['t'], D[ds]['gfp'][:,tr], '-', label='GFP (measured)', color='#009900', linewidth=.5)\n",
    "\n",
    "    # Format plot\n",
    "    ax.set_xlabel('Time [h]')\n",
    "    ax.set_ylabel('Number of molecules [10³]')\n",
    "    ax.set_title('{} #{:03d}\\n(separate fit)'.format(getDataLabel(D[ds]), tr))\n",
    "    ax.legend(handles=[p_dg, p_fg, p_tg, p_dr, p_fr, p_tr])\n",
    "\n",
    "    # Write figure to pdf\n",
    "    if pdf != None:\n",
    "        pdf.savefig(fig, bbox_inches='tight')\n",
    "\n",
    "    # Show and close figure\n",
    "    if i_plot < 200:\n",
    "        i_plot += 1\n",
    "        plt.show(fig)\n",
    "    fig.clear()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Fit traces separately\n",
    "for ds in range(len(D)):\n",
    "    nTraces = np.shape(D[ds]['rfp'])[1]\n",
    "\n",
    "    for tr in range(nTraces):\n",
    "        print('Fitting „{}“ #{:03d}/{:03d} …'.format(getDataLabel(D[ds]), tr, nTraces))\n",
    "        \n",
    "        # Prepare data\n",
    "        time = D[ds]['t']\n",
    "        data_red = D[ds]['rfp'][:,tr].flatten()\n",
    "        data_green = D[ds]['gfp'][:,tr].flatten()\n",
    "\n",
    "        wght_red = D[ds]['rfp_error'][:,tr]**2\n",
    "        wght_green = D[ds]['gfp_error'][:,tr]**2\n",
    "        \n",
    "        # Adjust parameter properties for onset time and offset\n",
    "        red_p.set('tr', max=time[data_red.argmax()])\n",
    "        red_p.set('offr',\n",
    "                       min=data_red[:10].min(),\n",
    "                       max=data_red[:10].max(),\n",
    "                       value=np.median(data_red[:10]))\n",
    "        green_p.set('tg', max=time[data_green.argmax()])\n",
    "        green_p.set('offg',\n",
    "                       min=data_green[:10].min(),\n",
    "                       max=data_green[:10].max(),\n",
    "                       value=np.median(data_green[:10]))\n",
    "\n",
    "        # Assess fixed parameters\n",
    "        red_fixed = red_p.df.loc[red_p.df.index[red_p.df.loc[:,'fixed']],'value'].to_dict()\n",
    "        green_fixed = green_p.df.loc[green_p.df.index[green_p.df.loc[:,'fixed']],'value'].to_dict()\n",
    "\n",
    "        # Objective function (closure)\n",
    "        def objective_red(params):\n",
    "            \"\"\"Objective function for RFP in separate model\"\"\"\n",
    "            cur_val = red_p.eval(params, t=time)\n",
    "            lik = np.sum(.5 * (data_red - cur_val)**2 / wght_red)\n",
    "            return lik\n",
    "\n",
    "        # Jacobian/gradient (closure)\n",
    "        def gradient_red(params):\n",
    "            \"\"\"Gradient for RFP in separate model\"\"\"\n",
    "            J = red_jacobian(**red_p.eval_params(params, t=time))\n",
    "            residuals = (data_red - red_p.eval(params, t=time)).reshape((np.size(time),1))\n",
    "            vrnc = wght_red.reshape(np.shape(residuals))\n",
    "            return -np.sum(J * residuals / vrnc, axis=0).flatten()\n",
    "\n",
    "        # Fit the data\n",
    "        result = sc.optimize.minimize(objective_red,\n",
    "                                      red_p.initial(),\n",
    "                                      method='TNC',# one of: 'SLSQP' 'TNC' 'L-BFGS-B'\n",
    "                                      bounds=red_p.bounds(),\n",
    "                                      jac=gradient_red,\n",
    "                                      options={'disp':True,\n",
    "                                               'maxiter': 10000}\n",
    "                                     )\n",
    "\n",
    "        # Print result\n",
    "        print(\"\\tRed success {}: {}\".format(result.success, result.message))\n",
    "\n",
    "        # Save results to R\n",
    "        res_red = red_p.free_values(result.x)\n",
    "        R[ds]['red']['params'].iloc[tr] = res_red#result.x\n",
    "        best_fit = red(time, **res_red, **red_fixed)\n",
    "        R[ds]['red']['fit'][:,tr] = best_fit\n",
    "        R[ds]['red']['success'][tr] = result.success\n",
    "        R[ds]['red']['chisq'][tr] = np.sum((best_fit - data_red)**2)\n",
    "\n",
    "        # Fit green data\n",
    "        def objective_green(params):\n",
    "            \"\"\"Objective function for green model\"\"\"\n",
    "            cur_val = green_p.eval(params, t=time)\n",
    "            lik = np.sum(.5 * (data_green - cur_val)**2 / wght_green)\n",
    "            return lik\n",
    "\n",
    "        def gradient_green(params):\n",
    "            \"\"\"Gradient for GFP in separate model\"\"\"\n",
    "            J = green_jacobian(**green_p.eval_params(params, t=time))\n",
    "            residuals = (data_green - green_p.eval(params, t=time)).reshape((np.size(time),1))\n",
    "            vrnc = wght_green.reshape(np.shape(residuals))\n",
    "            return -np.sum(J * residuals / vrnc, axis=0).flatten()\n",
    "\n",
    "        result = sc.optimize.minimize(objective_green,\n",
    "                                      green_p.initial(),\n",
    "                                      method='TNC',\n",
    "                                      bounds=green_p.bounds(),\n",
    "                                      jac=gradient_green,\n",
    "                                      options={'disp': True,\n",
    "                                               'maxiter': 10000})\n",
    "        print(\"\\tGreen success {}: {}\".format(result.success, result.message))\n",
    "\n",
    "        res_green = green_p.free_values(result.x)\n",
    "        R[ds]['green']['params'].iloc[tr] = res_green\n",
    "        best_fit = green(time, **res_green, **green_fixed)\n",
    "        R[ds]['green']['fit'][:,tr] = best_fit\n",
    "        R[ds]['green']['success'][tr] = result.success\n",
    "        R[ds]['green']['chisq'][tr] = np.sum((best_fit - data_green)**2)\n",
    "\n",
    "        # DEBUG\n",
    "        #if tr >= 2:\n",
    "        #    print(\"Breaking loop for debugging purposes\")\n",
    "        #    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot results of separate fit\n",
    "ts = getTimeStamp()\n",
    "i_plot = 0\n",
    "for ds in range(len(D)):\n",
    "    par_kde = {}\n",
    "    for t in ('red', 'green'):\n",
    "        par_kde[t] = parameter_KDE(R[ds][t]['params'])\n",
    "    pdffile = os.path.join(getOutpath(), '{}_separate_{}.pdf'.format(ts, getDataLabel(D[ds], True)))\n",
    "    with PdfPages(pdffile) as pdf:\n",
    "        for tr in range(np.shape(D[ds]['rfp'])[1]):\n",
    "            plotSeparate(ds, tr, pdf, par_kde)\n",
    "\n",
    "            # DEBUG\n",
    "            #if tr >= 2:\n",
    "            #    print(\"Break loop\")\n",
    "            #    break"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit and plot combined model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotCombined(ds, tr, pdf=None, par_kde=None):\n",
    "    \"\"\"Fits and plots the data, treating RFP and GFP together.\n",
    "    \n",
    "    Keyword arguments:\n",
    "    ds -- the dictionary key of the dataset\n",
    "    tr -- the index of the trace in the dataset to be processed\n",
    "    pdf -- a PdfPages object to which the figure is written if it is not None\n",
    "    par_kde -- if containing dict of values of parameter distributions, plot distributions\n",
    "    \"\"\"\n",
    "\n",
    "    # Plot fit results\n",
    "    fig = plt.figure()\n",
    "    #fig.set_tight_layout(False)\n",
    "\n",
    "    if par_kde != None:\n",
    "        fig.set_figwidth(1.6 * fig.get_figwidth())\n",
    "\n",
    "        pn_both = ['m_ktl']\n",
    "        pn_red = [p for p in ['tr', 'kmr', 'betr', 'deltr', 'offr']\n",
    "                 if not red_p.df.loc[p, 'fixed']]\n",
    "        pn_green = [p for p in ['tg', 'kmg', 'betg', 'deltg', 'offg']\n",
    "                   if not green_p.df.loc[p, 'fixed']]\n",
    "\n",
    "        # Plot combined parameters\n",
    "        grid = (2, 1+max(len(pn_red), len(pn_green)))\n",
    "        gs = GridSpec(grid[0], grid[1])\n",
    "\n",
    "        #for pi, label in enumerate(pn_both):\n",
    "        pi = 0\n",
    "        label = pn_both[pi]\n",
    "        ax = plt.subplot(gs.new_subplotspec((pi, 0), rowspan=2))\n",
    "        data = par_kde['combined'][label]\n",
    "        clr_face = '#0000ff55'\n",
    "        clr_edge = '#000099ff'\n",
    "        curr_val = R[ds]['combined']['params'].loc[tr,label]\n",
    "        plot_kde(ax, data, label, clr_face, clr_edge, curr_val)\n",
    "\n",
    "        # Plot green parameters\n",
    "        for pi, label in enumerate(pn_green):\n",
    "            ax = plt.subplot(gs.new_subplotspec((0, pi+1)))\n",
    "            data = par_kde['combined'][label]\n",
    "            clr_face = '#00ff0055'\n",
    "            clr_edge = '#009900ff'\n",
    "            curr_val = R[ds]['combined']['params'].loc[tr,label]\n",
    "            plot_kde(ax, data, label, clr_face, clr_edge, curr_val)\n",
    "\n",
    "        # Plot red parameters\n",
    "        for pi, label in enumerate(pn_red):\n",
    "            ax = plt.subplot(gs.new_subplotspec((1, pi+1)))\n",
    "            data = par_kde['combined'][label]\n",
    "            clr_face = '#ff000055'\n",
    "            clr_edge = '#990000ff'\n",
    "            curr_val = R[ds]['combined']['params'].loc[tr,label]\n",
    "            plot_kde(ax, data, label, clr_face, clr_edge, curr_val)\n",
    "\n",
    "        # Adjust subplot layout\n",
    "        gs.tight_layout(fig, pad=0, rect=(0.5, 0, 1, 1))\n",
    "\n",
    "        # Create axes for fit\n",
    "        gs_fit = GridSpec(1, 1)\n",
    "        ax = fig.add_subplot(gs_fit[0])\n",
    "        gs_fit.tight_layout(fig, pad=0, rect=(0, 0, 0.5, 1))\n",
    "\n",
    "    else:\n",
    "        ax = fig.gca()\n",
    "\n",
    "    #wr = np.sqrt(D[ds]['rfp'][:,tr])\n",
    "    #ax.fill_between(D[ds]['t'], D[ds]['rfp'][:,tr]-wr, D[ds]['rfp'][:,tr]+wr, color='#ff000033')\n",
    "    #wg = np.sqrt(D[ds]['gfp'][:,tr])\n",
    "    #ax.fill_between(D[ds]['t'], D[ds]['gfp'][:,tr]-wg, D[ds]['gfp'][:,tr]+wg, color='#00ff0033')\n",
    "\n",
    "    p_tr = ax.axvline(R[ds]['combined']['params']['tr'][tr], label='RFP onset',\n",
    "                       color='#ff0000', linewidth=.5, linestyle='--')\n",
    "    p_tg = ax.axvline(R[ds]['combined']['params']['tg'][tr], label='GFP onset',\n",
    "                      color='#00ff00', linewidth=.5, linestyle='--')\n",
    "    p_fr, = ax.plot(D[ds]['t'], R[ds]['combined']['fit']['red'][tr], '-', label='RFP (fit)', color='#ff0000', linewidth=1)\n",
    "    p_fg, = ax.plot(D[ds]['t'], R[ds]['combined']['fit']['green'][tr], '-', label='GFP (fit)', color='#00ff00', linewidth=1)\n",
    "    p_dr, = ax.plot(D[ds]['t'], D[ds]['rfp'][:,tr], '-', label='RFP (measured)', color='#990000', linewidth=.5)\n",
    "    p_dg, = ax.plot(D[ds]['t'], D[ds]['gfp'][:,tr], '-', label='GFP (measured)', color='#009900', linewidth=.5)\n",
    "\n",
    "    # Format plot\n",
    "    ax.set_xlabel('Time [h]')\n",
    "    ax.set_ylabel('Fluorescence intensity [a.u.]')\n",
    "    ax.set_title('{} {} [{}] #{:03d}\\n(combined fit)'.format(\n",
    "        D[ds]['sample'], D[ds]['condition'], D[ds]['measurement'], tr))\n",
    "    ax.legend(handles=[p_dg, p_fg, p_tg, p_dr, p_fr, p_tr])\n",
    "\n",
    "    # Write figure to pdf\n",
    "    if pdf != None:\n",
    "        pdf.savefig(fig, bbox_inches='tight')\n",
    "\n",
    "    # Show and close figure\n",
    "    plt.show(fig)\n",
    "    fig.clear()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Fit combined model\n",
    "combined_fixed = combined_p.fixed_params()\n",
    "for ds in range(len(D)):\n",
    "    R[ds]['combined']['fit'] = {'red': [], 'green': []}\n",
    "    nTraces = np.shape(D[ds]['rfp'])[1]\n",
    "\n",
    "    for tr in range(nTraces):\n",
    "        print('Fitting „{}“ #{:03d}/{:03d}. '.format(getDataLabel(D[ds]), tr, nTraces))#, end='')\n",
    "\n",
    "        # Get the data for fitting\n",
    "        time = D[ds]['t']\n",
    "        data = np.stack([D[ds]['rfp'][:,tr], D[ds]['gfp'][:,tr]], axis=1)\n",
    "        wght = np.stack([D[ds]['rfp_error'][:,tr], D[ds]['gfp_error'][:,tr]], axis=1)**2\n",
    "        \n",
    "        # Adjust parameter properties for onset time and offset\n",
    "        combined_p.set('tr', max=time[data[:,0].argmax()])\n",
    "        combined_p.set('offr',\n",
    "                       min=data[:10,0].min(),\n",
    "                       max=data[:10,0].max(),\n",
    "                       value=np.median(data[:10,0]))\n",
    "        combined_p.set('tg', max=time[data[:,1].argmax()])\n",
    "        combined_p.set('offg',\n",
    "                       min=data[:10,1].min(),\n",
    "                       max=data[:10,1].max(),\n",
    "                       value=np.median(data[:10,1]))\n",
    "\n",
    "        # Get amplitude correction\n",
    "        amp_red = data[:,0].max() - data[:,0].min()\n",
    "        amp_green = data[:,1].max() - data[:,1].min()\n",
    "        amp_correct = amp_red - amp_green\n",
    "\n",
    "        # Fit the data\n",
    "        def objective_fcn(params):\n",
    "            \"\"\"Objective function for combined model\"\"\"\n",
    "            cur_val = combined_p.eval(params, t=time)\n",
    "            chisq = np.sum(.5 * (data - cur_val)**2 / wght)\n",
    "            return chisq\n",
    "\n",
    "        def gradient_combined(params):\n",
    "            \"\"\"Gradient for combined model\"\"\"\n",
    "            J_red = red_jacobian(**{p: v for p, v in combined_p.eval_params(params, t=time).items()\n",
    "                                    if p in {'t', 'm_ktl', 'tr', 'kmr', 'betr', 'deltr', 'offr'}})\n",
    "            J_green = green_jacobian(**{p: v for p, v in combined_p.eval_params(params, t=time).items()\n",
    "                                    if p in {'t', 'm_ktl', 'tg', 'kmg', 'betg', 'deltg', 'offg'}})\n",
    "            residuals = (data - combined_p.eval(params, t=time)) / wght\n",
    "            Jr = -np.sum(residuals[:,np.newaxis,0] * J_red, axis=0)\n",
    "            Jg = -np.sum(residuals[:,np.newaxis,1] * J_green, axis=0)\n",
    "\n",
    "            # Assemble gradient according to parameter order:\n",
    "            # (tr, tg, m_ktl, kmr, kmg, betr, betg, deltr, deltg, offr, offg)\n",
    "            #return np.array([Jr[0], Jg[0], Jr[1]+Jg[1], Jr[2], Jg[2], Jr[3], Jg[3],\n",
    "            #                 Jr[4], Jg[4], Jr[5], Jg[5]])\n",
    "            return np.array([Jr[0], Jg[0], Jr[1]+Jg[1], Jr[2], Jg[2], Jr[3], Jg[3]])\n",
    "\n",
    "        result = sc.optimize.minimize(objective_fcn,\n",
    "                                      combined_p.initial(),\n",
    "                                      method='TNC',#'L-BFGS-B','TNC'\n",
    "                                      bounds=combined_p.bounds(),\n",
    "                                      jac=gradient_combined,\n",
    "                                      options={'disp': True,\n",
    "                                               'maxiter': 10000})\n",
    "\n",
    "        # Save results to R\n",
    "        res_combined = combined_p.free_values(result.x)\n",
    "        best_fit = combined(time, **res_combined, **combined_fixed)\n",
    "        tr_idx = R[ds]['combined']['params'].index[tr]\n",
    "\n",
    "        R[ds]['combined']['params'].loc[tr_idx, list(res_combined.keys())] = res_combined\n",
    "        #R[ds]['combined']['params'].loc[tr_idx, combined_fixed.keys()] = combined_fixed\n",
    "\n",
    "        R[ds]['combined']['fit']['red'].insert(tr, best_fit[:,0])\n",
    "        R[ds]['combined']['fit']['green'].insert(tr, best_fit[:,1])\n",
    "\n",
    "        R[ds]['combined']['success'][tr] = result.success\n",
    "        R[ds]['combined']['chisq'][tr,:] = np.sum((best_fit - data)**2, axis=0)\n",
    "\n",
    "        # Print result\n",
    "        print(\"\\tSuccess {} after {} iterations: {}\".format(\n",
    "            result.success, result.nit, result.message))\n",
    "\n",
    "        # DEBUG\n",
    "        #if tr >= 20:\n",
    "        #    print(\"Breaking loop for debugging purposes\")\n",
    "        #    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot results of combined fit\n",
    "ts = getTimeStamp()\n",
    "for ds in range(len(D)):\n",
    "    par_kde = {}\n",
    "    for t in ('combined',):\n",
    "        par_kde[t] = parameter_KDE(R[ds][t]['params'])\n",
    "    pdffile = os.path.join(getOutpath(), '{}_combined_{}.pdf'.format(ts, getDataLabel(D[ds], True)))\n",
    "    with PdfPages(pdffile) as pdf:\n",
    "        for tr in range(np.shape(D[ds]['rfp'])[1]):\n",
    "            plotCombined(ds, tr, pdf, par_kde=par_kde)\n",
    "            \n",
    "            # DEBUG\n",
    "            #if tr >= 40:\n",
    "            #    print(\"Forcing break\")\n",
    "            #    break"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Compare chisquare\n",
    "print(\"     {:10s} {:10s}\".format('red', 'green'))\n",
    "for i in range(nTraces):\n",
    "    print(\"{:03d}: {:10.0f} {:10.0f}\".format(\n",
    "        i,\n",
    "        R[0]['red']['chisq'][i] - R[0]['combined']['chisq'][i,0],\n",
    "        R[0]['green']['chisq'][i] - R[0]['combined']['chisq'][i,1]\n",
    "    ))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter correlations ($t_0$, $mk_\\text{tl}$, $\\delta$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set this to True to plot all figures in this section log-log, else set to False\n",
    "want_loglog = True\n",
    "want_ellipse = True\n",
    "want_trim = False\n",
    "want_clean = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def do_PCA(p1, p2, n_sigma=1, draw=True, isLog=False):\n",
    "    \"\"\"\n",
    "    Calculates properties of a confidence ellipse by PCA.\n",
    "\n",
    "    Input:\n",
    "        p1 -- 1dim array of observations of first (horizontal) dimension\n",
    "        p2 -- 1dim array of observations of second (vertical) dimension\n",
    "        n_sigma -- number of confidence intervals to use as ellipse radius\n",
    "        draw -- boolean flag whether to create Ellipse artist\n",
    "        isLog -- if True, all returned values refer to the log10 of p1 and p2\n",
    "\n",
    "    Returns dictionary with fields:\n",
    "        center -- tuple (x,y) of ellipse center coordinates\n",
    "        width -- total (horizontal) width of ellipse\n",
    "        height -- total (vertical) height of ellipse\n",
    "        angle -- rotation angle of ellipse in degrees\n",
    "        artist -- matplotlib.patches.Ellipse artist, or None if draw==False\n",
    "    \"\"\"\n",
    "    # Get logarithm of values\n",
    "    if isLog:\n",
    "        p1 = np.log10(p1)\n",
    "        p2 = np.log10(p2)\n",
    "\n",
    "    # Get eigenvalues and eigenvectors of covariance matrix\n",
    "    eigvals, eigvecs = np.linalg.eig(np.cov(p1, p2))\n",
    "\n",
    "    # Sort eigenvalues and eigenvectors\n",
    "    if eigvals[1] > eigvals[0]:\n",
    "        eigvals = eigvals[::-1]\n",
    "        eigvecs = eigvecs[:,::-1]\n",
    "\n",
    "    # Calculate ellipse properties\n",
    "    widths = n_sigma * np.sqrt(eigvals)\n",
    "    theta = np.degrees(np.arctan2(*eigvecs[::-1,0]))\n",
    "    mean1 = np.mean(p1)\n",
    "    mean2 = np.mean(p2)\n",
    "\n",
    "    # Build ellipse\n",
    "    if draw:\n",
    "        ellipse = mptch.Ellipse(xy=(mean1, mean2), width=widths[0],\n",
    "                                height=widths[1], angle=theta)\n",
    "    else:\n",
    "        ellipse = None\n",
    "\n",
    "    return {\"center\": (mean1, mean2),\n",
    "            \"width\": widths[0],\n",
    "            \"height\": widths[1],\n",
    "            \"angle\": theta,\n",
    "            \"artist\": ellipse}\n",
    "\n",
    "def construct_ellipse(width, height, center=(0,0), theta=0, isLog=False):\n",
    "    \"\"\"\n",
    "    Calculates points on the ellipse line.\n",
    "\n",
    "    Input:\n",
    "        width -- length of horizontal semi-axis\n",
    "        height -- length of vertical semi-axis\n",
    "        center -- tuple of horizontal and vertical coordinates of ellipse center\n",
    "        theta -- angle (in degrees) of ellipse rotation\n",
    "        isLog -- if True, transform logarithm for \n",
    "    \"\"\"\n",
    "    # Calculate points on ellipse\n",
    "    phi = np.linspace(0, 2 * np.pi, 200)\n",
    "    vals = np.matrix([width * np.cos(phi), height * np.sin(phi)])\n",
    "\n",
    "    # Set up rotation matrix\n",
    "    theta_rad = np.deg2rad(theta)\n",
    "    c = np.cos(theta_rad)\n",
    "    s = np.sin(theta_rad)\n",
    "    R = np.matrix([[c, -s], [s, c]])\n",
    "\n",
    "    # Rotate and shift ellipse\n",
    "    vals_trafo = (R * vals + np.matrix([center]).T).A\n",
    "\n",
    "    # If logarithmic data were given, raise with basis 10\n",
    "    if isLog:\n",
    "        vals_trafo = np.power(10, vals_trafo)\n",
    "\n",
    "    # Ensure that first and last point are numerically equal\n",
    "    vals_trafo[:,-1] = vals_trafo[:,0]\n",
    "\n",
    "    return vals_trafo[0,:], vals_trafo[1,:]\n",
    "\n",
    "def eliminate_outliers(iData, thresh_chsq=3.6, thresh_delta=1e-2):\n",
    "    \"\"\"\n",
    "    Eliminates outliers based on heuristic thresholds.\n",
    "    A logical index array is returned, in which outliers are marked\n",
    "    as False and traces to keep as True.\n",
    "\n",
    "    Input:\n",
    "        iData -- index of dataset in D and R\n",
    "        thresh_chsq -- chisquare threshold: sort out traces above threshold\n",
    "        thresh_delta -- delta threshold: sort out traces below threshold\n",
    "\n",
    "    Returns:\n",
    "        idx -- logical array with as many elements as traces\n",
    "    \"\"\"\n",
    "    # Obtain chi-square values\n",
    "    cqr = R[iData]['red']['chisq']\n",
    "    cqg = R[iData]['green']['chisq']\n",
    "    if cqr.shape != cqg.shape:\n",
    "        raise ValueError(\"Incompatible shapes of RFP and GFP data.\")\n",
    "\n",
    "    # Get index based on chi-square thresholding\n",
    "    rawr = D[iData]['rfp']\n",
    "    rawg = D[iData]['gfp']\n",
    "    qmr = np.log10(cqr / rawr.max(axis=0))\n",
    "    qmg = np.log10(cqg / rawg.max(axis=0))\n",
    "    idx = (qmr < thresh_chsq) & (qmg < thresh_chsq)\n",
    "\n",
    "    # Sort out outliers based on small delta value\n",
    "    deltr = R[iData]['red']['params']['deltr'].values\n",
    "    deltg = R[iData]['green']['params']['deltg'].values\n",
    "    idx &= (deltr > thresh_delta) & (deltg > thresh_delta)\n",
    "\n",
    "    return idx\n",
    "\n",
    "def trimOutliers(x, log=False):\n",
    "    \"\"\"Returns a logical index of all non-outliers.\"\"\"\n",
    "    if log:\n",
    "        x = np.log10(x)\n",
    "    p90 = np.percentile(x, 95)\n",
    "    thrsh = p90 * 1.5\n",
    "    i_in = x <= thrsh\n",
    "    if log:\n",
    "        p05 = np.percentile(x, 5)\n",
    "        i_in = np.logical_and(i_in, x >= p05 - .5 * (p90 - p05))\n",
    "    return i_in\n",
    "\n",
    "def plotParamCorrelations(p_green, p_red, lbl='', param_name='',\n",
    "                          param_unit='a.u.', loglog=False, pdf=None,\n",
    "                          trim=True, fit_diag=False, equal_aspect=False,\n",
    "                          p2_green=None, p2_red=None, desc1=None, desc2=None,\n",
    "                          draw_ellipse=True):\n",
    "    \"\"\"\n",
    "    Plots the parameter time correlations for dataset `ds`\n",
    "    and saves the plot to PDF if `pdf` is a `PdfPages` instance.\n",
    "\n",
    "    Arguments:\n",
    "        p_green -- 1-dim array of parameter values for GFP\n",
    "        p_red -- array of parameter values for RFP, same shape as p_green\n",
    "        lbl -- dataset label for figure title\n",
    "        param_name -- parameter name, used for figure labels\n",
    "        param_unit -- unit of parameter, used in axes labels\n",
    "        loglog -- if True, scale both axes logarithmically\n",
    "        pdf -- PdfPages object to save figure to (None for no saving)\n",
    "        trim -- if True, trim outliers\n",
    "        fit_diag -- if True, fit and plot diagonal (slope 1) to data\n",
    "        equal_aspect -- if True, use same scaling for x- and y-axis\n",
    "        p2_green -- 1-dim array of second dataset for GFP\n",
    "        p2_red --1-dim array of second dataset for RFP, same shape as p2_green\n",
    "        desc1 -- short label for first dataset\n",
    "        desc2 -- short label for second dataset\n",
    "        draw_ellipse -- if True, draw a confidence ellipse\n",
    "    \"\"\"\n",
    "    # Define constants\n",
    "    col1 = \"#1f77b4\"\n",
    "    col2 = \"#ff7f0e\"\n",
    "\n",
    "    if desc1 is not None:\n",
    "        descc1 = \"{}: \".format(desc1)\n",
    "    else:\n",
    "        descc1 = \"\"\n",
    "    if desc2 is not None:\n",
    "        descc2 = \"{}: \".format(desc2)\n",
    "    else:\n",
    "        descc2 = \"\"\n",
    "\n",
    "    # Test for second dataset\n",
    "    hasSecond = False\n",
    "    if p2_green is not None and p2_red is not None:\n",
    "        hasSecond = True\n",
    "    \n",
    "    # Trim outliers\n",
    "    out_lbl = \"\"\n",
    "    out_lbl2 = \"\"\n",
    "    if trim:\n",
    "        ing = trimOutliers(p_green, log=loglog)\n",
    "        inr = trimOutliers(p_red, log=loglog)\n",
    "        if hasSecond:\n",
    "            ing2 = trimOutliers(p2_green, log=loglog)\n",
    "            inr2 = trimOutliers(p2_red, log=loglog)\n",
    "\n",
    "            # Get combined value ranges\n",
    "            maxg = max(p_green[ing].max(), p2_green[ing2].max())\n",
    "            ming = min(p_green[ing].min(), p2_green[ing2].min())\n",
    "            maxr = max(p_red[inr].max(), p2_red[inr2].max())\n",
    "            minr = min(p_red[inr].min(), p2_red[inr2].min())\n",
    "\n",
    "            # Trim outliers outside of combined ranges\n",
    "            in_all = np.all(np.stack((p_green >= ming, p_green <= maxg,\n",
    "                                      p_red >= minr, p_red <= maxr)), axis=0)\n",
    "            p_green = p_green[in_all]\n",
    "            p_red = p_red[in_all]\n",
    "            in_all2 = np.all(np.stack((p2_green >= ming, p2_green <= maxg,\n",
    "                                      p2_red >= minr, p2_red <= maxr)), axis=0)\n",
    "            p2_green = p2_green[in_all2]\n",
    "            p2_red = p2_red[in_all2]\n",
    "\n",
    "            # Set outlier labels\n",
    "            n_out = np.sum(~in_all)\n",
    "            out_lbl = \" ($+${} outlier{})\".format(n_out, '' if n_out == 1 else 's')\n",
    "            n_out2 = np.sum(~in_all2)\n",
    "            out_lbl2 = \" ($+${} outlier{})\".format(n_out2, '' if n_out2 == 1 else 's')\n",
    "\n",
    "        elif not np.all(ing) or not np.all(inr):\n",
    "            in_all = np.logical_and(ing, inr)\n",
    "            n_out = np.sum(np.logical_not(in_all))\n",
    "            p_green = p_green[in_all]\n",
    "            p_red = p_red[in_all]\n",
    "            out_lbl = \"\\n($+${} outlier{})\".format(n_out, '' if n_out == 1 else 's')\n",
    "\n",
    "    # Make figure and axes\n",
    "    f, ax = plt.subplots(figsize=(4, 4))\n",
    "\n",
    "    # Plot data\n",
    "    ax.plot(p_green, p_red, '.', color=col1, ms=2, mec='none', zorder=3)\n",
    "    if hasSecond:\n",
    "        ax.plot(p2_green, p2_red, '.', color=col2, ms=2, mec='none', zorder=3)\n",
    "\n",
    "    # Fix and get axes limits\n",
    "    if loglog:\n",
    "        ax.set_xscale('log')\n",
    "        ax.set_yscale('log')\n",
    "    ax.autoscale(enable=False)\n",
    "    if equal_aspect:\n",
    "        ax.set_aspect('equal')\n",
    "    xlim = np.array(ax.get_xlim())\n",
    "    ylim = np.array(ax.get_ylim())\n",
    "\n",
    "    # Fit and plot line\n",
    "    if fit_diag:\n",
    "        offset = so.least_squares(lambda p: p_green - p_red + p, np.zeros(1)).x\n",
    "        diag_x_vals = np.linspace(xlim.min(), xlim.max(), 200)\n",
    "        diag_y_vals = diag_x_vals + offset\n",
    "        ax.plot(diag_x_vals, diag_y_vals, '-', color=col1, lw=.75, zorder=2)\n",
    "\n",
    "    # Draw confidence ellipse\n",
    "    if draw_ellipse:\n",
    "        ellinfo = do_PCA(p_green, p_red, isLog=loglog)\n",
    "        ell_x, ell_y = construct_ellipse(width=ellinfo['width'],\n",
    "                        height=ellinfo['height'], center=ellinfo['center'],\n",
    "                        theta=ellinfo['angle'], isLog=loglog)\n",
    "        ax.plot(ell_x, ell_y, '-', lw=.5, color=col1, zorder=4)\n",
    "        ctr = ellinfo['center']\n",
    "        if loglog:\n",
    "            ctr = np.power(10, ctr)\n",
    "        ax.plot(*ctr, marker='+', mec=col1, mew=.75, ms=5, zorder=4)\n",
    "        if hasSecond:\n",
    "            ellinfo = do_PCA(p2_green, p2_red, isLog=loglog)\n",
    "            ell_x, ell_y = construct_ellipse(width=ellinfo['width'],\n",
    "                            height=ellinfo['height'], center=ellinfo['center'],\n",
    "                            theta=ellinfo['angle'], isLog=loglog)\n",
    "            ax.plot(ell_x, ell_y, '-', lw=.5, color=col2, zorder=4)\n",
    "            ctr = ellinfo['center']\n",
    "            if loglog:\n",
    "                ctr = np.power(10, ctr)\n",
    "            ax.plot(*ctr, marker='+', mec=col2, mew=.75, ms=5, zorder=4)\n",
    "\n",
    "    # Draw diagonal\n",
    "    diag_lim = (max(xlim.min(), ylim.min()), min(xlim.max(), ylim.max()))\n",
    "    ax.plot(diag_lim, diag_lim, '--k', lw=.5, zorder=1)\n",
    "\n",
    "    # Show legend\n",
    "    if hasSecond:\n",
    "        dummy1 = mlin.Line2D([], [], marker='.', ls='', color=col1, label=desc1)\n",
    "        dummy2 = mlin.Line2D([], [], marker='.', ls='', color=col2, label=desc2)\n",
    "        ax.legend(handles=(dummy1, dummy2),\n",
    "                  loc=('lower left' if fit_diag else 'lower right'),\n",
    "                  borderpad=.1, borderaxespad=.1, labelspacing=.2, handletextpad=.5)\n",
    "        #ax.legend(handles=(dummy1, dummy2), ncol=2, mode='expand',\n",
    "        #         loc='lower left', bbox_to_anchor=(0, 1, 1, 0.1))\n",
    "\n",
    "    # Write offset\n",
    "    if fit_diag:\n",
    "        ax.text(xlim.max()-.1, ylim.min()+.1,\n",
    "                \"$t_\\mathrm{{r}}-t_\\mathrm{{g}}={:.2f}$\".format(\n",
    "                    offset.item()), {'color': col1, 'ha': 'right', 'va': 'bottom'})\n",
    "\n",
    "    # Write statistics information\n",
    "    stat_text = \"{}{} cells\".format(descc1, p_green.size) + out_lbl\n",
    "    if hasSecond:\n",
    "        stat_text += \"\\n{}{} cells\".format(descc2, p2_green.size) + out_lbl2\n",
    "    ax.text(xlim.min() * 1.01, ylim.max() * .99, stat_text,\n",
    "           {'color': 'k', 'ha': 'left', 'va': 'top'})\n",
    "\n",
    "    # Axes formatting\n",
    "    ax.set_xlabel(\"{} for GFP [{}]\".format(param_name, param_unit), color='g')\n",
    "    ax.set_ylabel(\"{} for RFP [{}]\".format(param_name, param_unit), color='r')\n",
    "    ax.set_title('GFP/RFP correlation of {}\\n{}'.format(param_name, lbl))\n",
    "    f.tight_layout(pad=0)\n",
    "\n",
    "    # Show, optionally save, and close the figure\n",
    "    plt.show(f)\n",
    "    if pdf is not None:\n",
    "        pdf.savefig(f)\n",
    "    f.clear()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Onset time correlations for separate files\n",
    "with PdfPages(getOutpath(\"onset_correlations.pdf\")) as pdf:\n",
    "    for i in range(len(R)):\n",
    "        # Load relevant data\n",
    "        if want_clean:\n",
    "            idx = eliminate_outliers(i)\n",
    "        else:\n",
    "            idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "        t_green = R[i]['green']['params']['tg'].values[idx]\n",
    "        t_red = R[i]['red']['params']['tr'].values[idx]\n",
    "        lbl = getDataLabel(D[i])\n",
    "\n",
    "        #plotOnsetCorrelations(t_green, t_red, lbl, pdf)\n",
    "        plotParamCorrelations(t_green, t_red, lbl=lbl, param_name='Onset time',\n",
    "                param_unit='h', fit_diag=True, equal_aspect=True, pdf=pdf,\n",
    "                loglog=False, trim=want_trim, draw_ellipse=want_ellipse)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Correlate onsets of merged datasets\n",
    "onset_sorted = []\n",
    "samples = set()\n",
    "conditions = [\"control\", \"siRNA\"]\n",
    "\n",
    "# Write relevant data in list\n",
    "for i in range(len(R)):\n",
    "    if want_clean:\n",
    "        idx = eliminate_outliers(i)\n",
    "    else:\n",
    "        idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "    onset_sorted.append({\n",
    "        'green': R[i]['green']['params']['tg'].values[idx],\n",
    "        'red': R[i]['red']['params']['tr'].values[idx],\n",
    "        'sample': D[i]['sample'],\n",
    "        'condition': D[i]['condition']\n",
    "    })\n",
    "\n",
    "    samples.add(onset_sorted[i]['sample'])\n",
    "    if onset_sorted[i]['condition'] not in conditions:\n",
    "        conditions.append(onset_sorted[i]['condition'])\n",
    "\n",
    "# Iterate over relevant data list\n",
    "with PdfPages(getOutpath(\"onset_correlations_merged.pdf\")) as pdf:\n",
    "    for sample in samples:\n",
    "        for condition in conditions:\n",
    "            # Get the indices for this property combination\n",
    "            idx = [i for i in range(len(onset_sorted))\n",
    "                   if onset_sorted[i]['sample'] == sample\n",
    "                   and onset_sorted[i]['condition'] == condition]\n",
    "            #if len(idx) == 0:\n",
    "            #    continue\n",
    "\n",
    "            # Merge corresponding datasets\n",
    "            t_red = np.concatenate([onset_sorted[i]['red'] for i in idx])\n",
    "            t_green = np.concatenate([onset_sorted[i]['green'] for i in idx])\n",
    "            lbl = \"{} ({})\".format(sample, condition)\n",
    "\n",
    "            # Plot merged datasets\n",
    "            plotParamCorrelations(t_green, t_red, lbl=lbl, param_name='Onset time',\n",
    "                    param_unit='h', fit_diag=True, equal_aspect=True, pdf=pdf,\n",
    "                    loglog=want_loglog, trim=want_trim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlate onsets of merged datasets in one figure\n",
    "onset_sorted = []\n",
    "samples = set()\n",
    "conditions = [\"control\", \"siRNA\"]\n",
    "\n",
    "# Write relevant data in list\n",
    "for i in range(len(R)):\n",
    "    if want_clean:\n",
    "        idx = eliminate_outliers(i)\n",
    "    else:\n",
    "        idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "    onset_sorted.append({\n",
    "        'green': R[i]['green']['params']['tg'].values[idx],\n",
    "        'red': R[i]['red']['params']['tr'].values[idx],\n",
    "        'sample': D[i]['sample'],\n",
    "        'condition': D[i]['condition']\n",
    "    })\n",
    "\n",
    "    samples.add(onset_sorted[i]['sample'])\n",
    "    if onset_sorted[i]['condition'] not in conditions:\n",
    "        conditions.append(onset_sorted[i]['condition'])\n",
    "\n",
    "# Iterate over relevant data list\n",
    "with PdfPages(getOutpath(\"onset_correlations_merged_log.pdf\")) as pdf:\n",
    "    for sample in samples:\n",
    "        lbl = \"{}\".format(sample)\n",
    "        p_red = []\n",
    "        p_green = []\n",
    "        desc = []\n",
    "\n",
    "        for condition in conditions:\n",
    "            # Get the indices for this property combination\n",
    "            idx = [i for i in range(len(onset_sorted))\n",
    "                   if onset_sorted[i]['sample'] == sample\n",
    "                   and onset_sorted[i]['condition'] == condition]\n",
    "            if len(idx) == 0:\n",
    "                continue\n",
    "\n",
    "            # Merge corresponding datasets\n",
    "            p_red.append(np.concatenate([onset_sorted[i]['red'] for i in idx]))\n",
    "            p_green.append(np.concatenate([onset_sorted[i]['green'] for i in idx]))\n",
    "            desc.append(condition)\n",
    "\n",
    "        # Plot merged datasets\n",
    "        plotParamCorrelations(p_green[0], p_red[0], lbl=lbl, param_name='Onset time',\n",
    "                param_unit='h', fit_diag=True, equal_aspect=True, pdf=pdf,\n",
    "                loglog=False, trim=want_trim,\n",
    "                p2_green=p_green[1], p2_red=p_red[1], desc1=desc[0], desc2=desc[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# m_ktl for separate files\n",
    "with PdfPages(getOutpath(\"m-ktl_correlations.pdf\")) as pdf:\n",
    "    for i in range(len(R)):\n",
    "        # Load relevant data\n",
    "        if want_clean:\n",
    "            idx = eliminate_outliers(i)\n",
    "        else:\n",
    "            idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "        p_green = R[i]['green']['params']['m_ktl'].values[idx]\n",
    "        p_red = R[i]['red']['params']['m_ktl'].values[idx]\n",
    "        lbl = getDataLabel(D[i])\n",
    "\n",
    "        plotParamCorrelations(p_green, p_red, lbl=lbl, pdf=pdf,\n",
    "                              param_name=r\"$mk_\\mathrm{tl}$\",\n",
    "                              param_unit=r\"$\\mathregular{h^{-1}}$\",\n",
    "                              loglog=True, trim=want_trim)\n",
    "\n",
    "# delta for separate files\n",
    "with PdfPages(getOutpath(\"delta_correlations.pdf\")) as pdf:\n",
    "    for i in range(len(R)):\n",
    "        # Load relevant data\n",
    "        if want_clean:\n",
    "            idx = eliminate_outliers(i)\n",
    "        else:\n",
    "            idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "        p_green = R[i]['green']['params']['deltg'].values[idx]\n",
    "        p_red = R[i]['red']['params']['deltr'].values[idx]\n",
    "        lbl = getDataLabel(D[i])\n",
    "\n",
    "        plotParamCorrelations(p_green, p_red, lbl=lbl, pdf=pdf,\n",
    "                              param_name=r\"$\\delta$\",\n",
    "                              param_unit=r\"$\\mathregular{h^{-1}}$\",\n",
    "                              loglog=want_loglog, trim=want_trim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Correlate m_ktl and delta of merged datasets\n",
    "mktl_sorted = []\n",
    "delta_sorted = []\n",
    "samples = set()\n",
    "conditions = [\"control\", \"siRNA\"]\n",
    "\n",
    "# Write relevant data in list\n",
    "for i in range(len(R)):\n",
    "    sample = D[i]['sample']\n",
    "    condition = D[i]['condition']\n",
    "\n",
    "    samples.add(sample)\n",
    "    if condition not in conditions:\n",
    "        conditions.append(condition)\n",
    "\n",
    "    if want_clean:\n",
    "        idx = eliminate_outliers(i)\n",
    "    else:\n",
    "        idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "    mktl_sorted.append({\n",
    "        'green': R[i]['green']['params']['m_ktl'].values[idx],\n",
    "        'red': R[i]['red']['params']['m_ktl'].values[idx],\n",
    "        'sample': sample,\n",
    "        'condition': condition\n",
    "    })\n",
    "    delta_sorted.append({\n",
    "        'green': R[i]['green']['params']['deltg'].values[idx],\n",
    "        'red': R[i]['red']['params']['deltr'].values[idx],\n",
    "        'sample': sample,\n",
    "        'condition': condition\n",
    "    })\n",
    "\n",
    "# Iterate over merged m_ktl list\n",
    "with PdfPages(getOutpath(\"m-ktl_correlations_merged.pdf\")) as pdf:\n",
    "    for sample in samples:\n",
    "        for condition in conditions:\n",
    "            # Get the indices for this property combination\n",
    "            idx = [i for i in range(len(mktl_sorted))\n",
    "                   if mktl_sorted[i]['sample'] == sample\n",
    "                   and mktl_sorted[i]['condition'] == condition]\n",
    "            #if len(idx) == 0:\n",
    "            #    continue\n",
    "\n",
    "            # Merge corresponding datasets\n",
    "            p_red = np.concatenate([mktl_sorted[i]['red'] for i in idx])\n",
    "            p_green = np.concatenate([mktl_sorted[i]['green'] for i in idx])\n",
    "            lbl = \"{} ({})\".format(sample, condition)\n",
    "\n",
    "            # Plot merged datasets\n",
    "            plotParamCorrelations(p_green, p_red, lbl=lbl, pdf=pdf,\n",
    "                    param_name=r\"$mk_\\mathrm{tl}$\",\n",
    "                    param_unit=r\"$\\mathregular{h^{-1}}$\",\n",
    "                    loglog=want_loglog, trim=want_trim)\n",
    "\n",
    "# Iterate over merged delta list\n",
    "with PdfPages(getOutpath(\"delta_correlations_merged.pdf\")) as pdf:\n",
    "    for sample in samples:\n",
    "        for condition in conditions:\n",
    "            # Get the indices for this property combination\n",
    "            idx = [i for i in range(len(delta_sorted))\n",
    "                   if delta_sorted[i]['sample'] == sample\n",
    "                   and delta_sorted[i]['condition'] == condition]\n",
    "            #if len(idx) == 0:\n",
    "            #    continue\n",
    "\n",
    "            # Merge corresponding datasets\n",
    "            p_red = np.concatenate([delta_sorted[i]['red'] for i in idx])\n",
    "            p_green = np.concatenate([delta_sorted[i]['green'] for i in idx])\n",
    "            lbl = \"{} ({})\".format(sample, condition)\n",
    "\n",
    "            # Plot merged datasets\n",
    "            plotParamCorrelations(p_green, p_red, lbl=lbl, pdf=pdf,\n",
    "                    param_name=r\"$\\delta$\",\n",
    "                    param_unit=r\"$\\mathregular{h^{-1}}$\",\n",
    "                    loglog=want_loglog, trim=want_trim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Correlate m_ktl and delta of merged datasets in one figure\n",
    "mktl_sorted = []\n",
    "delta_sorted = []\n",
    "samples = set()\n",
    "conditions = [\"control\", \"siRNA\"]\n",
    "\n",
    "# Write relevant data in list\n",
    "for i in range(len(R)):\n",
    "    sample = D[i]['sample']\n",
    "    condition = D[i]['condition']\n",
    "\n",
    "    samples.add(sample)\n",
    "    if condition not in conditions:\n",
    "        conditions.append(condition)\n",
    "\n",
    "    if want_clean:\n",
    "        idx = eliminate_outliers(i)\n",
    "    else:\n",
    "        idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "    mktl_sorted.append({\n",
    "        'green': R[i]['green']['params']['m_ktl'].values[idx],\n",
    "        'red': R[i]['red']['params']['m_ktl'].values[idx],\n",
    "        'sample': sample,\n",
    "        'condition': condition\n",
    "    })\n",
    "    delta_sorted.append({\n",
    "        'green': R[i]['green']['params']['deltg'].values[idx],\n",
    "        'red': R[i]['red']['params']['deltr'].values[idx],\n",
    "        'sample': sample,\n",
    "        'condition': condition\n",
    "    })\n",
    "\n",
    "# Iterate over merged m_ktl list\n",
    "with PdfPages(getOutpath(\"m-ktl_correlations_merged_log.pdf\")) as pdf:\n",
    "    for sample in samples:\n",
    "        lbl = \"{}\".format(sample)\n",
    "        p_red = []\n",
    "        p_green = []\n",
    "        desc = []\n",
    "\n",
    "        for condition in conditions:\n",
    "            # Get the indices for this property combination\n",
    "            idx = [i for i in range(len(mktl_sorted))\n",
    "                   if mktl_sorted[i]['sample'] == sample\n",
    "                   and mktl_sorted[i]['condition'] == condition]\n",
    "            #if len(idx) == 0:\n",
    "            #    continue\n",
    "\n",
    "            # Merge corresponding datasets\n",
    "            p_red.append(np.concatenate([mktl_sorted[i]['red'] for i in idx]))\n",
    "            p_green.append(np.concatenate([mktl_sorted[i]['green'] for i in idx]))\n",
    "            desc.append(condition)\n",
    "\n",
    "        # Plot merged datasets\n",
    "        plotParamCorrelations(p_green[0], p_red[0], lbl=lbl, pdf=pdf,\n",
    "                param_name=r\"$mk_\\mathrm{tl}$\",\n",
    "                param_unit=r\"$\\mathregular{h^{-1}}$\",\n",
    "                loglog=want_loglog,  trim=want_trim,\n",
    "                equal_aspect=True,\n",
    "                p2_green=p_green[1], p2_red=p_red[1], desc1=desc[0], desc2=desc[1])\n",
    "\n",
    "# Iterate over merged delta list\n",
    "with PdfPages(getOutpath(\"delta_correlations_merged_log.pdf\")) as pdf:\n",
    "    for sample in samples:\n",
    "        lbl = \"{}\".format(sample)\n",
    "        p_red = []\n",
    "        p_green = []\n",
    "        desc = []\n",
    "\n",
    "        for condition in conditions:\n",
    "            # Get the indices for this property combination\n",
    "            idx = [i for i in range(len(delta_sorted))\n",
    "                   if delta_sorted[i]['sample'] == sample\n",
    "                   and delta_sorted[i]['condition'] == condition]\n",
    "\n",
    "            # Merge corresponding datasets\n",
    "            p_red.append(np.concatenate([delta_sorted[i]['red'] for i in idx]))\n",
    "            p_green.append(np.concatenate([delta_sorted[i]['green'] for i in idx]))\n",
    "            desc.append(condition)\n",
    "\n",
    "        # Plot merged datasets\n",
    "        plotParamCorrelations(p_green[0], p_red[0], lbl=lbl, pdf=pdf,\n",
    "                param_name=r\"$\\delta$\",\n",
    "                param_unit=r\"$\\mathregular{h^{-1}}$\",\n",
    "                loglog=want_loglog,  trim=want_trim,\n",
    "                equal_aspect=True,\n",
    "                p2_green=p_green[1], p2_red=p_red[1], desc1=desc[0], desc2=desc[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up merged parameter lists\n",
    "# Stores all parameters in a data structure accessible as\n",
    "# parameters[<sample>][<condition>][<parameter>][<color>]\n",
    "par_names = (\"onset\", \"mktl\", \"delta\")\n",
    "colors = (\"red\", \"green\")\n",
    "parameters = {}\n",
    "\n",
    "# Write relevant data in list\n",
    "for i, r in enumerate(R):\n",
    "    if want_clean:\n",
    "        idx = eliminate_outliers(i)\n",
    "    else:\n",
    "        idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "    sample = D[i]['sample']\n",
    "    condition = D[i]['condition']\n",
    "\n",
    "    if sample not in parameters:\n",
    "        parameters[sample] = {}\n",
    "    if condition not in parameters[sample]:\n",
    "        parameters[sample][condition] = {}\n",
    "    for p in par_names:\n",
    "        if p not in parameters[sample][condition]:\n",
    "            parameters[sample][condition][p] = {}\n",
    "        for c in colors:\n",
    "            if c not in parameters[sample][condition][p]:\n",
    "                parameters[sample][condition][p][c] = []\n",
    "\n",
    "    parameters[sample][condition][\"onset\"][\"red\"] = np.concatenate(\n",
    "        (parameters[sample][condition][\"onset\"][\"red\"],\n",
    "         r['red']['params']['tr'].values[idx]))\n",
    "    parameters[sample][condition][\"onset\"][\"green\"] = np.concatenate(\n",
    "        (parameters[sample][condition][\"onset\"][\"green\"],\n",
    "         r['green']['params']['tg'].values[idx]))\n",
    "    parameters[sample][condition][\"mktl\"][\"red\"] = np.concatenate(\n",
    "        (parameters[sample][condition][\"mktl\"][\"red\"],\n",
    "         r['red']['params']['m_ktl'].values[idx]))\n",
    "    parameters[sample][condition][\"mktl\"][\"green\"] = np.concatenate(\n",
    "        (parameters[sample][condition][\"mktl\"][\"green\"],\n",
    "         r['green']['params']['m_ktl'].values[idx]))\n",
    "    parameters[sample][condition][\"delta\"][\"red\"] = np.concatenate(\n",
    "        (parameters[sample][condition][\"delta\"][\"red\"],\n",
    "         r['red']['params']['deltr'].values[idx]))\n",
    "    parameters[sample][condition][\"delta\"][\"green\"] = np.concatenate(\n",
    "        (parameters[sample][condition][\"delta\"][\"green\"],\n",
    "         r['green']['params']['deltg'].values[idx]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Population analysis\n",
    "\n",
    "Goal: test whether control and siRNA are equally distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ranksums(X, Y, roundDigits=None, debug=False):\n",
    "    \"\"\"Calculates the ranksums for the Wilcoxon-Mann-Whitney-U test.\n",
    "    For degenerate values, the mean rank is used.\n",
    "\n",
    "    Input:\n",
    "        X -- 1-dimensional array of first sample (need not be sorted)\n",
    "        Y -- 2-dimensional array of second sample (need not be sorted)\n",
    "        roundDigits -- number of digits to round to; no rounding if None\n",
    "        debug -- if True, the rank lists are printed\n",
    "\n",
    "    Returns:\n",
    "        tuple (ranksum of X, ranksum of Y)\n",
    "    \"\"\"\n",
    "    # Flatten, round and sort X and Y\n",
    "    if len(X.shape) > 1:\n",
    "        X = X.flatten()\n",
    "    if len(Y.shape) > 1:\n",
    "        Y = Y.flatten()\n",
    "    if roundDigits is not None:\n",
    "        X = X.round(decimals=roundDigits)\n",
    "        Y = Y.round(decimals=roundDigits)\n",
    "    XY = np.concatenate((X,Y))\n",
    "    XY.sort()\n",
    "\n",
    "    # Compute ranks\n",
    "    R = {}\n",
    "    iR = 0\n",
    "    xy_uq, xy_ct = np.unique(XY, return_counts=True)\n",
    "    for i, c in enumerate(xy_ct):\n",
    "        iR_old = iR\n",
    "        iR += c\n",
    "        R[xy_uq[i]] = np.arange(iR, iR_old, -1).sum() / c\n",
    "\n",
    "    # Compute rank sums\n",
    "    if debug:\n",
    "        xnew = []\n",
    "        ynew = []\n",
    "\n",
    "    sumX = 0\n",
    "    for x in X:\n",
    "        sumX += R[x]\n",
    "        if debug:\n",
    "            xnew.append(R[x])\n",
    "\n",
    "    sumY = 0\n",
    "    for y in Y:\n",
    "        sumY += R[y]\n",
    "        if debug:\n",
    "            ynew.append(R[y])\n",
    "\n",
    "    if debug:\n",
    "        print(\"X: {}\\nY:{}\".format(str(xnew), str(ynew)))\n",
    "\n",
    "    return sumX, sumY\n",
    "\n",
    "def criticalValue(U, n1, n2):\n",
    "    \"\"\"Computes the critical value for the Wilcoxon-Mann-Whitney test.\n",
    "    The critical value is estimated using an estimation formula\n",
    "    assuming large n1 and n2.\n",
    "    The null hypothesis cannot be rejected when the critical value\n",
    "    is inside the desired range of a gaussian distribution.\n",
    "\n",
    "    Input:\n",
    "        U -- the smaller Mann-Whitney-U test statistic of the two samples\n",
    "        n1 -- the number of elements in one of the samples\n",
    "        n2 -- the number of elements in the other sample\n",
    "\n",
    "    Returns:\n",
    "        Estimated critical value\n",
    "    \"\"\"\n",
    "    return (U - n1 * n2 / 2) / np.sqrt(n1 * n2 * (n1 + n2 + 1) / 12)\n",
    "\n",
    "def WilcoxonMannWhitney(X, Y, alpha=.05, roundDigits=None, verbose=False):\n",
    "    \"\"\"Wilcoxon-Mann-Whitney-U test\n",
    "\n",
    "    Input:\n",
    "        X -- first sample as 1-dim array\n",
    "        Y -- second sample as 1-dim array\n",
    "        alpha -- significance level\n",
    "        roundDigits -- if not None, X and Y are rounded; corresponds\n",
    "                       to the \"decimals\" parameter of ndarray.round\n",
    "        verbose -- if True, print detailed information\n",
    "\n",
    "    Returns:\n",
    "        True if null hypothesis cannot be rejected (X and Y are from\n",
    "        same population), False otherwise\n",
    "    \"\"\"\n",
    "    # Obtain ranksums\n",
    "    sumX, sumY = ranksums(X, Y, roundDigits=roundDigits)\n",
    "    nX = X.size\n",
    "    nY = Y.size\n",
    "\n",
    "    # Find critical value\n",
    "    Ux = sumX - nX * (nX + 1) / 2\n",
    "    Uy = sumY - nY * (nY + 1) / 2\n",
    "    Z = criticalValue(min(Ux, Uy), nX, nY)\n",
    "\n",
    "    # Get interval of non-critical values\n",
    "    intZ = ss.norm.ppf((alpha/2, 1-alpha/2))\n",
    "\n",
    "    # Get p-value\n",
    "    sgnZ = -1 if Z < 0 else 1\n",
    "    p = ss.norm.cdf(-sgnZ * Z) + 1 - ss.norm.cdf(sgnZ * Z)\n",
    "\n",
    "    # Get test result\n",
    "    isAccept = False\n",
    "    if Z > intZ[0] and Z < intZ[1]:\n",
    "        isAccept = True\n",
    "\n",
    "    # DEBUG\n",
    "    if verbose:\n",
    "        print(\"\"\"Ranksums: x={:6.0f}, y={:6.0f}\n",
    "U:        x={:6.3f}, y={:6.3f}\n",
    "Z={:.3f} {} [{:.3f}, {:.3f}]\n",
    "p-value={:.3f}, 𝛼={:.3f}\n",
    "Null hypothesis is {}.\"\"\".format(sumX, sumY, Ux, Uy, Z,\n",
    "                                 '∈' if isAccept else '∉',\n",
    "                                 intZ[0], intZ[1], p, alpha, isAccept))\n",
    "\n",
    "    return isAccept, p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array((1,2,3,4,4,6))\n",
    "Y = np.array((2,2,3,4,5,6))\n",
    "ranksums(X, Y, debug=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = np.array((1,2,3,4,4,6,1,1))\n",
    "Y = np.array((2,2,3,4,5,6,7,10))\n",
    "WilcoxonMannWhitney(X, Y, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "WilcoxonMannWhitney(parameters['A549']['control']['mktl']['green'],\n",
    "                    parameters['A549']['siRNA']['mktl']['green'],\n",
    "                    roundDigits=3, alpha=.01, verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def BrownForsythe(Y, alpha=.05, verbose=False):\n",
    "    \"\"\"Brown-Forsythe test\"\"\"\n",
    "    k = len(Y)\n",
    "    N = 0\n",
    "    Ni = np.zeros(k)\n",
    "    Zij = []\n",
    "    Zi = []\n",
    "    denominator = 0\n",
    "\n",
    "    for i, y in enumerate(Y):\n",
    "        Ni[i] = y.size\n",
    "        N += Ni[i]\n",
    "        my = np.median(y)\n",
    "        Zij.insert(i, abs(y - my))\n",
    "        Zi.insert(i, np.mean(Zij[i]))\n",
    "        denominator += np.sum((Zij[i] - Zi[i])**2)\n",
    "\n",
    "    Z = np.mean(np.concatenate(Zij))\n",
    "    numerator = np.sum(Ni * (Zi - Z)**2)\n",
    "\n",
    "    W = (N - k) / (k - 1) * numerator / denominator\n",
    "    maxW = ss.f.ppf((1 - alpha,), k - 1, N - k)[0]\n",
    "    p = 1 - ss.f.cdf(W, k - 1, N - k)\n",
    "\n",
    "    isAccept = False\n",
    "    if W <= maxW:\n",
    "        isAccept = True\n",
    "\n",
    "    if verbose:\n",
    "        print(\"W={:.3f} {} {:.3f}\\np-value={:.3f}, 𝛼={:.3f}\\nNull hypothesis is {}.\".format(\n",
    "            W, \"≤\" if isAccept else \"≰\", maxW, p, alpha, isAccept))\n",
    "\n",
    "    return isAccept"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "BrownForsythe((parameters['A549']['control']['mktl']['green'],\n",
    "              parameters['A549']['siRNA']['mktl']['green']),\n",
    "             verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "BrownForsythe((parameters['A549']['control']['mktl']['green'],\n",
    "              parameters['A549']['siRNA']['mktl']['green']),\n",
    "             verbose=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "samples = (\"Huh7\", \"A549\")\n",
    "par_names = (\"onset\", \"mktl\", \"delta\")\n",
    "colors = (\"red\", \"green\")\n",
    "alpha = .05\n",
    "\n",
    "for sample in samples:\n",
    "    for par_name in par_names:\n",
    "        for color in colors:\n",
    "            print(\"********************\\n{:9s}: {:5s}\\n{:9s}: {:5s}\\n{:9s}: {:5s}\\n\".format(\n",
    "                    \"Sample\", sample, \"Parameter\", par_name, \"Color\", color))\n",
    "            print(\"Wilcoxon-Mann-Whitney-U test:\")\n",
    "            WilcoxonMannWhitney(parameters[sample]['control'][par_name][color],\n",
    "                                parameters[sample][ 'siRNA' ][par_name][color],\n",
    "                                roundDigits=3, alpha=alpha, verbose=True)\n",
    "            print(\"\\nBrown-Forsythe test:\")\n",
    "            BrownForsythe((parameters[sample]['control'][par_name][color],\n",
    "                           parameters[sample][ 'siRNA' ][par_name][color]),\n",
    "                          alpha=alpha, verbose=True)\n",
    "            print(\"********************\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = (\"Huh7\",)\n",
    "par_names = (\"mktl\",)\n",
    "colors = (\"red\", \"green\")\n",
    "alpha = .05\n",
    "\n",
    "for sample in samples:\n",
    "    for par_name in par_names:\n",
    "        for color in colors:\n",
    "            print(\"********************\\n{:9s}: {:5s}\\n{:9s}: {:5s}\\n{:9s}: {:5s}\\n\".format(\n",
    "                    \"Sample\", sample, \"Parameter\", par_name, \"Color\", color))\n",
    "            print(\"Wilcoxon-Mann-Whitney-U test:\")\n",
    "            WilcoxonMannWhitney(parameters[sample]['control'][par_name][color],\n",
    "                                parameters[sample][ 'siRNA' ][par_name][color],\n",
    "                                roundDigits=3, alpha=alpha, verbose=True)\n",
    "            print(\"\\nBrown-Forsythe test:\")\n",
    "            BrownForsythe((parameters[sample]['control'][par_name][color],\n",
    "                           parameters[sample][ 'siRNA' ][par_name][color]),\n",
    "                          alpha=alpha, verbose=True)\n",
    "            print(\"********************\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Probe populations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compare delta distributions\n",
    "d_h = parameters[\"Huh7\"][\"siRNA\"][\"delta\"][\"green\"]\n",
    "d_a = parameters[\"A549\"][\"siRNA\"][\"delta\"][\"green\"]\n",
    "\n",
    "d_h = np.log10(d_h)\n",
    "d_a = np.log10(d_a)\n",
    "\n",
    "f, ax = plt.subplots()\n",
    "_ = ax.hist(d_h, bins=150)\n",
    "\n",
    "ax.set_xlabel(\"$\\log(\\delta)$\")\n",
    "ax.set_ylabel(\"Occurrences\")\n",
    "plt.show(f)\n",
    "plt.close(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "db = DBSCAN(eps=0.055, min_samples=4)\n",
    "lbls = db.fit(d_h.reshape((-1,1)))\n",
    "\n",
    "\n",
    "values = []\n",
    "for l in np.unique(db.labels_):\n",
    "    values.append(d_h[db.labels_ == l])\n",
    "\n",
    "colm = np.concatenate((((0,0,0),), matplotlib.cm.tab10.colors), axis=0)\n",
    "\n",
    "f, ax = plt.subplots()\n",
    "_ = ax.hist(values, bins=100, histtype='barstacked', log=True, color=colm[:len(values),:])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "d_h.reshape((1,-1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "matplotlib.cm.tab10.colors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colm = np.concatenate((((0,0,0),), matplotlib.cm.tab10.colors), axis=0)\n",
    "colm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Playground for parameter correlations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make list of deltas and mktls\n",
    "delta_sorted = []\n",
    "mktl_sorted = []\n",
    "onset_sorted = []\n",
    "samples = set()\n",
    "conditions = [\"control\", \"siRNA\"]\n",
    "\n",
    "# Write relevant data in list\n",
    "for i in range(len(R)):\n",
    "    sample = D[i]['sample']\n",
    "    condition = D[i]['condition']\n",
    "\n",
    "    samples.add(sample)\n",
    "    if condition not in conditions:\n",
    "        conditions.append(condition)\n",
    "\n",
    "    if want_clean:\n",
    "        idx = eliminate_outliers(i)\n",
    "    else:\n",
    "        idx = np.ones((D[i]['nTraces']), dtype=np.bool_)\n",
    "\n",
    "    delta_sorted.append({\n",
    "        'green': R[i]['green']['params']['deltg'].values[idx],\n",
    "        'red': R[i]['red']['params']['deltr'].values[idx],\n",
    "        'sample': sample,\n",
    "        'condition': condition\n",
    "    })\n",
    "    mktl_sorted.append({\n",
    "        'green': R[i]['green']['params']['m_ktl'].values[idx],\n",
    "        'red': R[i]['red']['params']['m_ktl'].values[idx],\n",
    "        'sample': sample,\n",
    "        'condition': condition\n",
    "    })\n",
    "    onset_sorted.append({\n",
    "        'green': R[i]['green']['params']['tg'].values[idx],\n",
    "        'red': R[i]['red']['params']['tr'].values[idx],\n",
    "        'sample': D[i]['sample'],\n",
    "        'condition': D[i]['condition']\n",
    "    })\n",
    "\n",
    "params = {}\n",
    "for sample in samples:\n",
    "    lbl = \"{}\".format(sample)\n",
    "    params[sample] = {}\n",
    "\n",
    "    for condition in conditions:\n",
    "        params[sample][condition] = {'delta': {}, 'mktl': {}, 't0': {}}\n",
    "\n",
    "        # Get the indices for this property combination\n",
    "        idx = [i for i in range(len(delta_sorted))\n",
    "               if delta_sorted[i]['sample'] == sample\n",
    "               and delta_sorted[i]['condition'] == condition]\n",
    "\n",
    "        # Merge corresponding datasets\n",
    "        params[sample][condition]['delta']['red'] = \\\n",
    "            np.concatenate([delta_sorted[i]['red'] for i in idx])\n",
    "        params[sample][condition]['delta']['green'] = \\\n",
    "            np.concatenate([delta_sorted[i]['green'] for i in idx])\n",
    "\n",
    "        params[sample][condition]['mktl']['red'] = \\\n",
    "            np.concatenate([mktl_sorted[i]['red'] for i in idx])\n",
    "        params[sample][condition]['mktl']['green'] = \\\n",
    "            np.concatenate([mktl_sorted[i]['green'] for i in idx])\n",
    "\n",
    "        params[sample][condition]['t0']['red'] = \\\n",
    "            np.concatenate([onset_sorted[i]['red'] for i in idx])\n",
    "        params[sample][condition]['t0']['green'] = \\\n",
    "            np.concatenate([onset_sorted[i]['green'] for i in idx])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Calculate Pearson correlation coefficient and save to XLSX\n",
    "samples = set(params.keys())\n",
    "conditions = set()\n",
    "par_names = set()\n",
    "colors = ['red', 'green']\n",
    "corr_lin_col_name = \"Pearson (linear)\"\n",
    "corr_log_col_name = \"Perason (logarithmic)\"\n",
    "\n",
    "# Collect keys\n",
    "for k1 in samples:\n",
    "    conditions.update(params[k1].keys())\n",
    "    for k2 in conditions:\n",
    "        par_names.update(params[k1][k2].keys())\n",
    "\n",
    "# Initialize table\n",
    "muli = pd.MultiIndex(names=(\"Sample\", \"Condition\", \"Parameter\"),\n",
    "                     levels=[[], [], []], labels=[[], [], []])\n",
    "par_corr_tab = pd.DataFrame(index=muli,\n",
    "                            columns=(corr_lin_col_name,corr_log_col_name),\n",
    "                            dtype=np.float_)\n",
    "\n",
    "# Populate table\n",
    "for k1, d1 in params.items():\n",
    "    for k2, d2 in d1.items():\n",
    "        for k3, d3 in d2.items():\n",
    "            log_idx = np.logical_and(d3['green'] > 0, d3['red'] > 0)\n",
    "            d3_log_green = np.log10(d3['green'][log_idx])\n",
    "            d3_log_red = np.log10(d3['red'][log_idx])\n",
    "\n",
    "            r_pearson_lin = ss.pearsonr(d3['green'], d3['red'])[0]\n",
    "            r_pearson_log = ss.pearsonr(d3_log_green, d3_log_red)[0]\n",
    "\n",
    "            print(\"{:4s} {:7s} {:5s}: pearson correlation: lin={: 5.3f} log={: 5.3f}\".format(\n",
    "                k1, k2, k3, r_pearson_lin, r_pearson_log))\n",
    "            par_corr_tab.loc[(k1, k2, k3), (corr_lin_col_name, corr_log_col_name)] = \\\n",
    "                (r_pearson_lin, r_pearson_log)\n",
    "\n",
    "# Save to excel file\n",
    "xlsx_file = getOutpath(\"CORR-COEFF\" + \".xlsx\")\n",
    "xlsx_writer = pd.ExcelWriter(xlsx_file, engine='xlsxwriter')\n",
    "par_corr_tab.to_excel(xlsx_writer, sheet_name=\"Correlation\")\n",
    "xlsx_writer.save()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "par_corr_tab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Kolmogorov-Smirnov test\n",
    "sample = \"Huh7\"\n",
    "res_red = ss.ks_2samp(params[sample][\"control\"][\"mktl\"][\"red\"],\n",
    "                      params[sample][\"siRNA\"][\"mktl\"][\"red\"])\n",
    "res_green = ss.ks_2samp(params[sample][\"control\"][\"mktl\"][\"green\"],\n",
    "                        params[sample][\"siRNA\"][\"mktl\"][\"green\"])\n",
    "print(res_red.pvalue)\n",
    "print(res_green.pvalue)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Print delta values conditionally\n",
    "sample = 'Huh7'\n",
    "condition = 'control'\n",
    "\n",
    "for i_delta, delta in enumerate(delta_sorted):\n",
    "    # Filter values\n",
    "    if not delta['sample'] == sample or not delta['condition'] == condition:\n",
    "        continue\n",
    "\n",
    "    #idx = eliminate_outliers(i_delta)\n",
    "    measurement = D[i_delta]['measurement']\n",
    "    deltas_r = R[i_delta]['red']['params']['deltr'].values#[idx]\n",
    "    deltas_g = R[i_delta]['green']['params']['deltg'].values#[idx]\n",
    "\n",
    "    for j_d in range(deltas_r.size):\n",
    "        if deltas_r[j_d] < 1e-2 or deltas_g[j_d] < 1e-2:\n",
    "            print(\"{:16s}({:3d}): {:8.5f} {:8.5f}\".format(\n",
    "                measurement, j_d, deltas_r[j_d], deltas_g[j_d]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\chi^2$ analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "i = 4 # 4 or 9\n",
    "r = R[i]\n",
    "d = D[i]\n",
    "print(getDataLabel(d))\n",
    "\n",
    "# Raw traces\n",
    "dr = d['rfp']\n",
    "dg = d['gfp']\n",
    "\n",
    "# Fitted traces\n",
    "\n",
    "\n",
    "# Chi-square (sum of squared residuals)\n",
    "cqr = r['red']['chisq']\n",
    "cqg = r['green']['chisq']\n",
    "\n",
    "# Logarithm of chi-square\n",
    "lcqr = np.log10(cqr)\n",
    "lcqg = np.log10(cqg)\n",
    "\n",
    "# Maxima of raw traces\n",
    "mr = dr.max(axis=0)\n",
    "mg = dg.max(axis=0)\n",
    "\n",
    "# Relative chi-square\n",
    "qmr = np.log10(cqr / mr)\n",
    "qmg = np.log10(cqg / mg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate fits\n",
    "pr = r['red']['params']\n",
    "pg = r['green']['params']\n",
    "\n",
    "fr = np.empty_like(dr)\n",
    "fg = np.empty_like(dg)\n",
    "\n",
    "for j in range(d['nTraces']):\n",
    "    fr[:,j] = red_p.eval(**pr.iloc[j,:].to_dict(), t=d['t'])\n",
    "    fg[:,j] = green_p.eval(**pg.iloc[j,:].to_dict(), t=d['t'])\n",
    "\n",
    "lnrr = np.log10(np.sum(((fr - dr) / mr)**2, axis=0))\n",
    "lnrg = np.log10(np.sum(((fg - dg) / mg)**2, axis=0))"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for j, val in enumerate(lnrr):\n",
    "    print(\"{:03d}: {:10g} {:10g}\".format(j, val, lnrg[j]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, 2)\n",
    "ax[0].hist(qmr, color='r', bins=80)\n",
    "ax[1].hist(qmg, color='g', bins=80)\n",
    "plt.show(f)\n",
    "plt.close(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "for j in range(D[i]['nTraces']):\n",
    "    #print(\"{:03d}: {:15f} {:15f}\".format(j, lcqr[j], lcqg[j]))\n",
    "    print(\"{:03d}: {:15f} {:15f}\".format(j, qmr[j], qmg[j]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.flatnonzero(qmr >= 3.6))\n",
    "print(np.flatnonzero(qmg >= 3.6))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "j = 166\n",
    "plt.plot(d['t'], fr[:,j], d['t'], dr[:,j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "j = 237\n",
    "plt.plot(d['t'], fg[:,j], d['t'], dg[:,j])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(1, 2)\n",
    "ax[0].plot(mr, lcqr, '.r', ms=1.5)\n",
    "ax[1].plot(mg, lcqg, '.g', ms=1.5)\n",
    "\n",
    "plt.show(f)\n",
    "plt.close(f)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "for i, r in enumerate(R):\n",
    "    cqr = r['red']['chisq']\n",
    "    cqg = r['green']['chisq']\n",
    "\n",
    "    f, ax = plt.subplots()\n",
    "\n",
    "    ax.boxplot((np.log10(cqr), np.log10(cqg)), labels=('red', 'green'))\n",
    "\n",
    "    plt.show(f)\n",
    "    plt.close(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameter distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot violin distributions of the data sets (both separate and combined)\n",
    "pn_both = ('m_ktl',)\n",
    "#pn_red = ('tr', 'kmr', 'betr', 'deltr', 'offr')\n",
    "#pn_green = ('tg', 'kmg', 'betg', 'deltg', 'offg')\n",
    "pn_red = ('tr', 'deltr', 'offr')\n",
    "pn_green = ('tg', 'deltg', 'offg')\n",
    "\n",
    "grid = (2, len(pn_both)+max(len(pn_red), len(pn_green)))\n",
    "\n",
    "with PdfPages(getOutpath('parameter_distributions.pdf')) as pdf:\n",
    "    for ds in range(len(D)):\n",
    "\n",
    "        par_kde = {}\n",
    "        fit_types = []\n",
    "\n",
    "        # Check for separate fit\n",
    "        if 'red' in R[ds] and 'green' in R[ds]:\n",
    "            hasSeparate = True\n",
    "            fit_types += ['red', 'green']\n",
    "        else:\n",
    "            hasSeparate = False\n",
    "\n",
    "        # Check for combined fit\n",
    "        par_kde_combined = {}\n",
    "        if 'combined' in R[ds]:\n",
    "            hasCombined = True\n",
    "            fit_types += ['combined']\n",
    "        else:\n",
    "            hasCombined = False\n",
    "\n",
    "        # Calculate parameter distributions\n",
    "        for t in fit_types:\n",
    "            par_kde[t] = parameter_KDE(R[ds][t]['params'])\n",
    "\n",
    "        # Plot parameter distributions\n",
    "        for typeName, hasType in zip(('separate', 'combined'), (hasSeparate, hasCombined)):\n",
    "            if not hasType:\n",
    "                continue\n",
    "\n",
    "            fig = plt.figure()\n",
    "            gs = GridSpec(grid[0], grid[1])\n",
    "\n",
    "            if typeName == 'combined':\n",
    "                # Combined fit; define specific settings\n",
    "                pn_green_temp = pn_green\n",
    "                pn_red_temp = pn_red\n",
    "                offset_both = len(pn_both)\n",
    "                kde_label_green = 'combined'\n",
    "                kde_label_red = 'combined'\n",
    "\n",
    "                # Plot combined parameters\n",
    "                for pi, label in enumerate(pn_both):\n",
    "                    ax = plt.subplot(gs.new_subplotspec((pi, 0), rowspan=2))\n",
    "                    data = par_kde['combined'][label]\n",
    "                    clr_face = '#0000ff55'\n",
    "                    #clr_edge = '#000099ff'\n",
    "                    plot_kde(ax, data, label, clr_face)\n",
    "            else:\n",
    "                # Separate fit; define specific settings\n",
    "                pn_green_temp = pn_both + pn_green\n",
    "                pn_red_temp = pn_both + pn_red\n",
    "                offset_both = 0\n",
    "                kde_label_green = 'green'\n",
    "                kde_label_red = 'red'\n",
    "\n",
    "            # Plot green parameters\n",
    "            for pi, par_label in enumerate(pn_green_temp):\n",
    "                ax = plt.subplot(gs.new_subplotspec((0, pi+offset_both)))\n",
    "                data = par_kde[kde_label_green][par_label]\n",
    "                clr_face = '#00ff0055'\n",
    "                #clr_edge = '#009900ff'\n",
    "                plot_kde(ax, data, par_label, clr_face)\n",
    "\n",
    "            # Plot red parameters\n",
    "            for pi, par_label in enumerate(pn_red_temp):\n",
    "                ax = plt.subplot(gs.new_subplotspec((1, pi+offset_both)))\n",
    "                data = par_kde[kde_label_red][par_label]\n",
    "                clr_face = '#ff000055'\n",
    "                #clr_edge = '#990000ff'\n",
    "                plot_kde(ax, data, par_label, clr_face)\n",
    "\n",
    "            # Show and close figure\n",
    "            fig.suptitle(getDataLabel(D[ds]) + \" (\" + typeName + \" fit)\")\n",
    "            fig.tight_layout(pad=0, rect=(0, 0, 1, .93))\n",
    "            pdf.savefig(fig, bbox_inches='tight')\n",
    "            plt.show(fig)\n",
    "            plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(par_kde['green'].keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Playground\n",
    "This section contains code that was/is used for developing ideas."
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "# Plot parameter correlations\n",
    "\n",
    "# Get parameters to be correlated\n",
    "par_cor = (('tr', 'tg'), ('m_ktl', 'm_ktl'), ('kmr', 'kmg'),\n",
    "           ('betr', 'betg'), ('deltr', 'deltg'), ('offr', 'offg'))\n",
    "\n",
    "with PdfPages(os.path.join(getOutpath(), '{:s}_both_parameter_correlations.pdf'.format(getTimeStamp()))) as pdf:\n",
    "    for pr, pg in par_cor:\n",
    "        # Get results of data sets\n",
    "        r1 = R[0]\n",
    "        r2 = R[1]\n",
    "\n",
    "        # Get parameter values\n",
    "        val1r = r1['red']['params'].loc[:,pr].values\n",
    "        val1g = r1['green']['params'].loc[:,pg].values\n",
    "        val2r = r2['red']['params'].loc[:,pr].values\n",
    "        val2g = r2['green']['params'].loc[:,pg].values\n",
    "\n",
    "        # Sort out outliers\n",
    "        idx1 = np.ones(np.size(val1r), dtype=np.bool_)\n",
    "        idx2 = np.ones(np.size(val2r), dtype=np.bool_)\n",
    "        #isr = valr.argsort()[-2:]\n",
    "        #isg = valg.argsort()[-2:]\n",
    "        #if valr[isr[0]] < 0.9 * valr[isr[1]]:\n",
    "        #    idx1[isr[1]] = False\n",
    "        #if valg[isg[0]] < 0.9 * valg[isg[1]]:\n",
    "        #    idx1[isg[1]] = False\n",
    "\n",
    "        # Plot\n",
    "        fig = plt.figure(figsize=(4.6,4))\n",
    "        tit = \"{} [{}]\\nCorrelation {} – {}\".format(\n",
    "            D[0]['sample'], D[0]['measurement'], pr, pg)\n",
    "        ax1 = fig.add_subplot(111)\n",
    "        ax1.set_xscale('log')\n",
    "        ax1.set_yscale('log')\n",
    "        ax1.plot(val1r[idx1], val1g[idx1], '.', label=D[0]['condition'])\n",
    "        ax1.plot(val2r[idx2], val2g[idx2], '.', label=D[1]['condition'])\n",
    "\n",
    "        ax1.set_autoscale_on(False)\n",
    "        lmt = np.array([ax1.get_xlim(), ax1.get_ylim()])\n",
    "        diag = (lmt[:,0].max(), lmt[:,1].min())\n",
    "        ax1.plot(diag, diag, '-k')\n",
    "\n",
    "        ax1.set_xlabel(pr, color='r')\n",
    "        ax1.set_ylabel(pg, color='g')\n",
    "        ax1.set_title(tit)\n",
    "        ax1.legend()\n",
    "\n",
    "        fig.tight_layout(pad=0)\n",
    "        plt.show(fig)\n",
    "        pdf.savefig(fig)\n",
    "        plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot parameter correlations\n",
    "\n",
    "# Get parameters to be correlated\n",
    "par_cor = (('tr', 'tg'), ('m_ktl', 'm_ktl'), ('kmr', 'kmg'),\n",
    "           ('betr', 'betg'), ('deltr', 'deltg'), ('offr', 'offg'))\n",
    "\n",
    "for i, r in enumerate(R):\n",
    "    with PdfPages(os.path.join(getOutpath(), '{:s}_{}_parameter_correlations.pdf'.format(getTimeStamp(), getDataLabel(D[i], True)))) as pdf:\n",
    "        for pr, pg in par_cor:\n",
    "            # Get parameter values\n",
    "            valr = r['red']['params'].loc[:,pr].values\n",
    "            valg = r['green']['params'].loc[:,pg].values\n",
    "\n",
    "            # Sort out outliers\n",
    "            idx1 = np.ones(np.size(valr), dtype=np.bool_)\n",
    "            #isr = valr.argsort()[-2:]\n",
    "            #isg = valg.argsort()[-2:]\n",
    "            #if valr[isr[0]] < 0.9 * valr[isr[1]]:\n",
    "            #    idx1[isr[1]] = False\n",
    "            #if valg[isg[0]] < 0.9 * valg[isg[1]]:\n",
    "            #    idx1[isg[1]] = False\n",
    "\n",
    "            # Plot\n",
    "            fig = plt.figure(figsize=(12,4))\n",
    "            tit = \"{}\\nCorrelation {} – {}\".format(getDataLabel(D[i]), pr, pg)\n",
    "            ax1 = fig.add_subplot(131)\n",
    "            ax1.set_xscale('log')\n",
    "            ax1.set_yscale('log')\n",
    "            ax1.plot(valr[idx1], valg[idx1], '.')\n",
    "\n",
    "            ax1.set_autoscale_on(False)\n",
    "            lmt = np.array([ax1.get_xlim(), ax1.get_ylim()])\n",
    "            diag = (lmt[:,0].max(), lmt[:,1].min())\n",
    "            ax1.plot(diag, diag, '-k')\n",
    "\n",
    "            ax1.set_xlabel(pr, color='r')\n",
    "            ax1.set_ylabel(pg, color='g')\n",
    "            ax1.set_title(tit)\n",
    "\n",
    "            fig.tight_layout(pad=0)\n",
    "            plt.show(fig)\n",
    "            pdf.savefig(fig)\n",
    "            plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot the parameter distributions for the datasets\n",
    "ds_keys = list(R.keys())\n",
    "ds_keys.sort()\n",
    "params = R[ds_keys[0]]['combined']['params'].columns\n",
    "grid = (len(params), len(ds_keys))\n",
    "i_col = 0\n",
    "\n",
    "pdffile = os.path.join(getOutpath(), '{:s}_parameters.pdf'.format(getTimeStamp()))\n",
    "with PdfPages(pdffile) as pdf:\n",
    "    fig = plt.figure()\n",
    "    fig.set_figheight(grid[0] * .8 * fig.get_figheight())\n",
    "    fig.set_figwidth(grid[1] * .8 * fig.get_figwidth())\n",
    "\n",
    "    for ds in ds_keys:\n",
    "        i_row = 0\n",
    "        for p in params:\n",
    "            ax = plt.subplot2grid(grid, (i_row, i_col))\n",
    "            ax.hist(R[ds]['combined']['params'][p], bins=100)\n",
    "            if i_row == grid[0] - 1:\n",
    "                ax.set_xlabel('Value [a.u.]')\n",
    "            if i_col == 0:\n",
    "                ax.set_ylabel('Occurrences [#]')\n",
    "            ax.set_title('{:s}: {:s}'.format(ds, p))\n",
    "            i_row += 1\n",
    "        i_col += 1\n",
    "\n",
    "    pdf.savefig(fig)\n",
    "    plt.show(fig)\n",
    "    plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Plot onset time correlations\n",
    "pdffile = os.path.join(getOutpath(), '{:s}_onset_correlations.pdf'.format(getTimeStamp()))\n",
    "with PdfPages(pdffile) as pdf:\n",
    "    for k in R.keys():\n",
    "        fig = plt.figure()\n",
    "        plt.plot([0, 30], [0, 30], 'k-')\n",
    "        plt.plot(R[k]['combined']['params']['tr'], R[k]['combined']['params']['tg'], '.')\n",
    "        plt.xlabel('Onset RFP [h]')\n",
    "        plt.ylabel('Onset GFP [h]')\n",
    "        plt.title(k)\n",
    "        pdf.savefig(fig)\n",
    "        plt.show()\n",
    "        plt.close()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Degradation rate ratio\n",
    "def plotHistograms(maxH):\n",
    "    Rkeys = sorted(R.keys())\n",
    "    for ds in Rkeys:\n",
    "        #deltg = R[ds]['green']['params']['deltg']\n",
    "        #deltr = R[ds]['red']['params']['deltr']\n",
    "        deltg = R[ds]['combined']['params']['deltg']\n",
    "        deltr = R[ds]['combined']['params']['deltr']\n",
    "        quot = deltg / deltr\n",
    "\n",
    "        fig = plt.figure()\n",
    "        plt.hist(quot, bins=150, range=(0, maxH))\n",
    "        plt.title(ds)\n",
    "        plt.xlabel('$\\delta_\\mathrm{green} / \\delta_\\mathrm{red}$ [a.u.]')\n",
    "        plt.ylabel('Occurrences [#]')\n",
    "        plt.show(fig)\n",
    "        plt.close(fig)\n",
    "\n",
    "wdg.interact(plotHistograms, maxH=wdg.IntSlider(\n",
    "    value=100, min=0, max=1000, step=10, description='Histogram maximum', continuous_update=False));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Fit distribution to degradation rate quotient histograms\n",
    "def gamma(x, p=2, b=1, s=10):\n",
    "    return s * b**p * x**(p-1) * np.exp(-b * x) / sc.special.gamma(p)\n",
    "\n",
    "def gamma2(x, p1=1.9, p2=2.1, b1=0.9, b2=1.1, s1=10, s2=10):\n",
    "    return gamma(x, p1, b1, s1) + gamma(x, p2, b2, s2)\n",
    "\n",
    "def weibull(x, lmbd=.2, k=2, s=10):\n",
    "    return s * lmbd * k * (lmbd * x)**(k - 1) * np.exp(- (lmbd * x)**k)\n",
    "\n",
    "def weibull2(x, lmbd1=.15, lmbd2=.25, k1=1.9, k2=2.1, s1=10, s2=10):\n",
    "    return weibull(x, lmbd=lmbd1, k=k1, s=s1) + weibull(x, lmbd=lmbd2, k=k2, s=s2)\n",
    "\n",
    "# Define models\n",
    "model_gamma = lm.Model(gamma)\n",
    "model_gamma.set_param_hint(name='p', min=.01)\n",
    "model_gamma.set_param_hint(name='b', min=.01)\n",
    "model_gamma.set_param_hint(name='s', min=1)\n",
    "\n",
    "model_gamma2 = lm.Model(gamma2)\n",
    "model_gamma2.set_param_hint(name='p1', min=.01)\n",
    "model_gamma2.set_param_hint(name='p2', min=.01)\n",
    "model_gamma2.set_param_hint(name='b1', min=.01)\n",
    "model_gamma2.set_param_hint(name='b2', min=.01)\n",
    "model_gamma2.set_param_hint(name='s1', min=1)\n",
    "model_gamma2.set_param_hint(name='s2', min=1)\n",
    "\n",
    "model_weibull = lm.Model(weibull)\n",
    "model_weibull.set_param_hint(name='lmbd', min=.001)\n",
    "model_weibull.set_param_hint(name='k', min=.001, max=5)\n",
    "model_weibull.set_param_hint(name='s', min=1)\n",
    "\n",
    "model_weibull2 = lm.Model(weibull2)\n",
    "model_weibull2.set_param_hint(name='lmbd1', min=.001)\n",
    "model_weibull2.set_param_hint(name='lmbd2', min=.001)\n",
    "model_weibull2.set_param_hint(name='k1', min=.001, max=5)\n",
    "model_weibull2.set_param_hint(name='k2', min=.001, max=5)\n",
    "model_weibull2.set_param_hint(name='s1', min=1)\n",
    "model_weibull2.set_param_hint(name='s2', min=1)\n",
    "\n",
    "maxH = 40\n",
    "\n",
    "with PdfPages(os.path.join(getOutpath(), '{:s}_degradation_distribution.pdf'.format(getTimeStamp()))) as pdf:\n",
    "    for ds in sorted(R.keys()):\n",
    "        # Calculate degradation rate quotient\n",
    "        deltg = R[ds]['combined']['params']['deltg']\n",
    "        deltr = R[ds]['combined']['params']['deltr']\n",
    "        quot = deltg / deltr\n",
    "\n",
    "        # Create histogram\n",
    "        fig = plt.figure()\n",
    "        ax = fig.add_subplot(1, 2, 1)\n",
    "        hist_val, hist_edg = ax.hist(quot, bins=70, range=(0, maxH), label='Histogram')[:2]\n",
    "        hist_ctr = (hist_edg[:-1] + hist_edg[1:]) / 2\n",
    "\n",
    "        # Fit models\n",
    "        result_g = model_gamma.fit(hist_val, x=hist_ctr)\n",
    "        result_g2 = model_gamma2.fit(hist_val, x=hist_ctr)\n",
    "        result_w = model_weibull.fit(hist_val, x=hist_ctr)\n",
    "        result_w2 = model_weibull2.fit(hist_val, x=hist_ctr)\n",
    "\n",
    "        # Select models\n",
    "        #print('gamma: {}'.format(result_g.chisqr))\n",
    "        #print('gamma2: {}'.format(result_g2.chisqr))\n",
    "        #print('weibull: {}'.format(result_w.chisqr))\n",
    "        #print('weibull2: {}'.format(result_w2.chisqr))\n",
    "\n",
    "        if result_g2.chisqr < .7 * result_g.chisqr:\n",
    "            res_g = result_g2\n",
    "            name_g = 'gamma2'\n",
    "        else:\n",
    "            res_g = result_g\n",
    "            name_g = 'gamma'\n",
    "\n",
    "        if result_w2.chisqr < .7 * result_w.chisqr:\n",
    "            res_w = result_w2\n",
    "            name_w = 'weibull2'\n",
    "        else:\n",
    "            res_w = result_w\n",
    "            name_w = 'weibull'\n",
    "\n",
    "        # Plot models\n",
    "        x = np.linspace(.1, 5, 100)\n",
    "        ax.plot(hist_ctr, res_g.best_fit, '-', label=name_g, color='orange')\n",
    "        ax.plot(hist_ctr, res_w.best_fit, '-', label=name_w, color='magenta')\n",
    "        ax.legend()\n",
    "        ax.set_xlabel('$\\delta_\\mathrm{green} / \\delta_\\mathrm{red}$ [a.u.]')\n",
    "        ax.set_ylabel('Counts [#]')\n",
    "        ax.set_title(ds)\n",
    "\n",
    "        # Print fit reports\n",
    "        rep = res_g.fit_report(show_correl=False) + '\\n' + res_w.fit_report(show_correl=False)\n",
    "        ax = fig.add_subplot(1, 2, 2)\n",
    "        ax.set_axis_off()\n",
    "        ax.text(0, 1, rep, ha='left', va='top', family='monospace', size=5.5)\n",
    "\n",
    "        # Display, save and close figure\n",
    "        plt.show(fig)\n",
    "        pdf.savefig(fig)\n",
    "        plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Scatter plot of degradation rates\n",
    "Rkeys = sorted(R.keys())\n",
    "for ds in Rkeys:\n",
    "    deltg = R[ds]['combined']['params']['deltg']\n",
    "    deltr = R[ds]['combined']['params']['deltr']\n",
    "\n",
    "    fig = plt.figure()\n",
    "    h = plt.plot(deltg, deltr, '.')\n",
    "    plt.title(ds)\n",
    "    plt.xlabel('$\\delta_\\mathrm{green}$ [a.u.]')\n",
    "    plt.ylabel('$\\delta_\\mathrm{red}$ [a.u.]')\n",
    "    plt.xscale('log')\n",
    "    plt.yscale('log')\n",
    "    plt.show(fig)\n",
    "    plt.close(fig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print medians of deltas per dataset\n",
    "for i, r in enumerate(R):\n",
    "    print('[{:31s}] deltg: {:5.3f}, deltr: {:5.3f}'.format(\n",
    "        getDataLabel(D[i]),\n",
    "        np.median(r[\"green\"][\"params\"][\"deltg\"].values),\n",
    "        np.median(r[\"red\"][\"params\"][\"deltr\"].values)))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}