profile-fake-data.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from itertools import product\n",
    "import os, sys\n",
    "import time\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "\n",
    "from flowstock import Akagi, ForwardSimulator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HiddenPrints:\n",
    "    def __enter__(self):\n",
    "        self._original_stdout = sys.stdout\n",
    "        sys.stdout = open(os.devnull, 'w')\n",
    "\n",
    "    def __exit__(self, exc_type, exc_val, exc_tb):\n",
    "        sys.stdout.close()\n",
    "        sys.stdout = self._original_stdout"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_params(num_cells, seed=0):\n",
    "\n",
    "    np.random.seed(seed)\n",
    "    \n",
    "    n = int(num_cells**0.5)\n",
    "    x = np.linspace(-1, 1, n)\n",
    "    y = np.linspace(-1, 1, n)\n",
    "    xv, yv = np.meshgrid(x, y, sparse=False, indexing='ij')\n",
    "\n",
    "    loc = []\n",
    "    for i in range(n):\n",
    "        for j in range(n):\n",
    "            loc.append([x[i], y[j]])\n",
    "\n",
    "    loc = np.array(loc)\n",
    "\n",
    "    rad = 0.8\n",
    "    N_init = np.exp(-((loc**2).sum(axis=1)**0.5 - rad)**2)\n",
    "    N_init = (N_init * 1000 * N_init.max()).astype(int)\n",
    "    N_init.reshape((n, n))\n",
    "\n",
    "    s = np.exp((-4 * loc**2).sum(axis=1))\n",
    "    s.reshape((n, n))\n",
    "\n",
    "    pi = N_init / N_init.max() / 10\n",
    "    pi.reshape((n, n))\n",
    "\n",
    "    # Array of distances between cells\n",
    "    d = np.zeros((num_cells, num_cells))\n",
    "    for i in range(num_cells):\n",
    "        for j in range(num_cells):\n",
    "            d[i, j] = ((loc[i] - loc[j])**2).sum()**0.5\n",
    "\n",
    "    beta = 1\n",
    "\n",
    "    K = 1.5\n",
    "    \n",
    "    return N_init, d, K, pi, s, beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def generate_fake_data(N_init, d, K, pi, s, beta, seed=0):\n",
    "\n",
    "    simulator = ForwardSimulator(pi, s, beta, d, K)\n",
    "    simulator.noise_amplitude = 0.0\n",
    "\n",
    "    # num_steps = 11\n",
    "    num_steps = 2\n",
    "    np.random.seed(seed)\n",
    "    N, M_true = simulator.simulate(N_init, num_steps)\n",
    "\n",
    "    return N, M_true"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def estimate(N, M, d, K, eps, lamda):\n",
    "    a = Akagi(N, d, K)\n",
    "    a.lamda = lamda\n",
    "    \n",
    "    result = a.exact_inference(eps)\n",
    "    \n",
    "    return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Precompile numba expressions\n",
    "def precompile():\n",
    "    num_cells = 121\n",
    "    params = generate_params(num_cells)\n",
    "    _, d, K, pi_true, s_true, beta_true = params\n",
    "\n",
    "    N_true, M_true = generate_fake_data(*params)\n",
    "\n",
    "    with HiddenPrints():\n",
    "        M_est, pi_est, s_est, beta_est = estimate(N_true*1.0, M_true, d, K, 0.1, 10.0)\n",
    "\n",
    "%time precompile()\n",
    "%time precompile()\n",
    "%time precompile()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def nae(tr, es):\n",
    "    return np.abs((tr - es)).sum() / np.abs(tr).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "\n",
    "columns = ['num regions', 'eps', 'lamda', 'scale', 'M_nae', 'pi_nae', 's_nae', 'beta_nae', 'time (s)']\n",
    "df = pd.DataFrame(columns=columns)\n",
    "\n",
    "for (num_cells, seed) in product(np.arange(6, 14)**2, range(3)):\n",
    "\n",
    "    print(\"Generate params, \", num_cells, \" cells\")\n",
    "    params = generate_params(num_cells, seed=seed)\n",
    "    _, d, K, pi_true, s_true, beta_true = params\n",
    "\n",
    "    print(\"Generate fake data\")\n",
    "    N_true, M_true = generate_fake_data(*params, seed=seed)\n",
    "\n",
    "    for (eps, lamda, scale) in product(np.logspace(-1, -4, 4), np.logspace(1, 2, 2), np.logspace(0, 3, 1)):\n",
    "        \n",
    "        lamda = float(lamda / scale)\n",
    "        \n",
    "        print(\"Estimate M from fake N, eps={eps}, lamda={lamda}\".format(eps=eps, lamda=lamda))\n",
    "        with HiddenPrints():\n",
    "            t_0 = time.time()\n",
    "            M_est, pi_est, s_est, beta_est = estimate(N_true*scale, M_true, d, K, eps, lamda)\n",
    "            t_1 = time.time()\n",
    "\n",
    "        M_est /= scale\n",
    "        \n",
    "        df.loc[len(df)] = [\n",
    "            num_cells,\n",
    "            eps, lamda, scale,\n",
    "            nae(M_true, M_est),\n",
    "            nae(pi_true, pi_est), nae(s_true/s_true.max(), s_est), nae(beta_true, beta_est), \n",
    "            t_1 - t_0,\n",
    "        ]"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "scrolled": true
   },
   "source": [
    "estimate(N_true, M_true, d, K, eps, float(lamda/scale))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['eps'] = df['eps'].astype('category')\n",
    "df['lamda'] = df['lamda'].astype('category')\n",
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['time (s)'].min()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df['time (s)'].max()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(2,1, figsize = (6, 10))\n",
    "\n",
    "# NAE\n",
    "sns.lineplot(\n",
    "    data=df,\n",
    "    x='num regions',\n",
    "    y='M_nae',\n",
    "    hue='eps',\n",
    "    style='eps',\n",
    "    ax=ax1,\n",
    ")\n",
    "\n",
    "ax1.set_ylim(0,)\n",
    "\n",
    "# Time\n",
    "ln = sns.lineplot(\n",
    "    data=df,\n",
    "    x='num regions',\n",
    "    y='time (s)',\n",
    "    hue='eps',\n",
    "    style='eps',\n",
    "    ax=ax2,\n",
    ")\n",
    "\n",
    "# ax2.set_ylim(0,)\n",
    "ax2.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.distplot(\n",
    "    df['time (s)'],\n",
    "    kde=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax_eps, ax_lamda) = plt.subplots(2, 1, figsize=(6, 10))\n",
    "\n",
    "sns.lineplot(\n",
    "    data=df,\n",
    "    x='eps',\n",
    "    y='M_nae',\n",
    "    ax=ax_eps,\n",
    ")\n",
    "ax_eps.set_xscale('log')\n",
    "ax_eps.set_ylim(0,)\n",
    "\n",
    "plot = sns.lineplot(\n",
    "    data=df,\n",
    "    x='lamda',\n",
    "    y='M_nae',\n",
    "    ax=ax_lamda,\n",
    ")\n",
    "ax_lamda.set_xscale('log')\n",
    "ax_lamda.set_ylim(0,)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot = sns.lineplot(\n",
    "    data=df,\n",
    "    x='lamda',\n",
    "    y='time (s)',\n",
    ")\n",
    "plot.axes.set_xscale('log')\n",
    "plot.axes.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot = sns.lineplot(\n",
    "    data=df,\n",
    "    x='eps',\n",
    "    y='time (s)',\n",
    ")\n",
    "plot.axes.set_xscale('log')\n",
    "plot.axes.set_yscale('log')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(\n",
    "    np.abs((pi_est-pi_true)).reshape(13, 13))\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(\n",
    "    np.abs(s_est-s_true/s_true.max()).reshape(13, 13)\n",
    ")\n",
    "plt.colorbar()"
   ]
  }
 ],
 "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}