create_sample_model_output.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample script to show each models output"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import seaborn as sns\n",
    "import os\n",
    "from pysb.logging import setup_logger\n",
    "from pysb.simulator import CudaSSASimulator, ScipyOdeSimulator, OpenCLSSASimulator, StochKitSimulator, BngSimulator\n",
    "from pysb.examples.earm_1_0 import model as earm_model\n",
    "from pysb.examples.schloegl import model as schoelgl_model\n",
    "from pysb.examples.kinase_cascade import model as kinase_model\n",
    "from pysb.examples.michment import model as michment_model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysb.bng import generate_equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "generate_equations(kinase_model)\n",
    "generate_equations(earm_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_model_stats(model):\n",
    "    generate_equations(model)\n",
    "    return {\n",
    "        'model_name': model.name,\n",
    "        'n_species': len(model.species),\n",
    "        'n_reactions': len(model.reactions),\n",
    "        'n_parameters': len(model.parameters)\n",
    "        \n",
    "    }"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_model_stats(kinase_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_model_stats(earm_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_model_stats(michment_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "get_model_stats(schoelgl_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tspan = np.linspace(0, 10, 101)\n",
    "sim = OpenCLSSASimulator(michment_model, tspan=tspan, verbose=True)\n",
    "traj = sim.run(number_sim=100).dataframe\n",
    "traj.reset_index(inplace=True)\n",
    "\n",
    "# ode = ScipyOdeSimulator(michment_model, tspan=tspan).run()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(traj.head(10))\n",
    "obs = ['E_free','S_free', 'ES_complex', 'Product']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_model(model, tspan, n_sim=1000, obs=None):\n",
    "    sim = OpenCLSSASimulator(model, tspan=tspan, verbose=True, precision=np.float64)\n",
    "\n",
    "    ode = ScipyOdeSimulator(model, tspan=tspan).run()\n",
    "\n",
    "#     sim = CUDASimulator(model, tspan=tspan, verbose=True, )\n",
    "    traj = sim.run(tspan=tspan, number_sim=n_sim)\n",
    "    plot(traj, ode, tspan, obs, obs)\n",
    "    return\n",
    "    for n, i in enumerate(model.species):\n",
    "        name = '__s{}'.format(n)\n",
    "        print(name)\n",
    "        plot(traj, ode, tspan, name, i)\n",
    "\n",
    "def plot(traj, ode, tspan, obs, title):\n",
    "    x = traj.dataframe[obs].unstack(0).values\n",
    "\n",
    "    plt.figure()\n",
    "    plt.title(title)\n",
    "    # create line traces\n",
    "    plt.plot(tspan, x, '0.5', lw=2, alpha=0.25)  # individual trajectories\n",
    "    plt.plot(tspan, x.mean(1), 'k-*', lw=3, label=\"Mean\")\n",
    "    plt.plot(tspan, x.min(1), 'b--', lw=3, label=\"Minimum\")\n",
    "    plt.plot(tspan, x.max(1), 'r--', lw=3, label=\"Maximum\")\n",
    "\n",
    "    # adding ODE solution to plot\n",
    "    plt.plot(tspan, ode.dataframe[name], 'g--', lw=3, label=\"ODE\")\n",
    "\n",
    "    plt.xlabel('Time(s)')\n",
    "    plt.ylabel('Number of molecules')\n",
    "    if not os.path.exists(\"figures\"):\n",
    "        os.mkdir('figures')\n",
    "    plt.savefig('figures/{}.png'.format(obs), bbox_inches='tight')\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "tspan = np.linspace(0, 8000, 101)\n",
    "name = 'tBid_total'\n",
    "# name = 'tBid_total'\n",
    "run_model(earm_model, tspan, n_sim=1000, obs=name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "tspan = np.linspace(0, 300, 101)\n",
    "name = 'ppERK'\n",
    "run_model(kinase_model, tspan, 100, obs=name)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tspan = np.linspace(0, 100, 101)\n",
    "name = 'X_total'\n",
    "\n",
    "run_model(scholgl_model, tspan, n_sim=2**12, obs=name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "def save_output(model, tspan, n_sim=1000):\n",
    "    name = model.name.split('.')[-1]\n",
    "    \n",
    "    def save_traj(traj, sim):\n",
    "        s_name = \"{name}_{n_sim}_{sim}.csv.gz\".format(name=name, n_sim=n_sim, sim=sim)\n",
    "        traj.to_csv(s_name, compression='gzip')\n",
    "        \n",
    "    print(\"Running OpenCLSimulator\")\n",
    "    traj = OpenCLSimulator(model, platform='NVIDIA',device='gpu'\n",
    "                          ).run(tspan=tspan, number_sim=n_sim).dataframe\n",
    "    save_traj(traj, 'opencl')\n",
    "    \n",
    "    print(\"Running CUDASimulator\")\n",
    "    traj = CUDASimulator(model).run(tspan=tspan, number_sim=n_sim).dataframe\n",
    "    save_traj(traj, 'cuda')\n",
    "    \n",
    "    print(\"Running BngSimulator\")\n",
    "    traj = BngSimulator(model).run(tspan=tspan, n_runs=n_sim).dataframe\n",
    "    save_traj(traj, 'bng')\n",
    "    \n",
    "    print(\"Running StochKitSimulator\")\n",
    "    traj = StochKitSimulator(model).run(tspan=tspan, n_runs=n_sim).dataframe\n",
    "    save_traj(traj, 'stochkit')\n",
    "    \n",
    "def load_traj(model, n_sim=1000):\n",
    "    \n",
    "    name = model.name.split('.')[-1]\n",
    "    bng_file = \"{name}_{n_sim}_{sim}.csv.gz\".format(name=name, n_sim=n_sim, sim='bng')\n",
    "    stochkit_file = \"{name}_{n_sim}_{sim}.csv.gz\".format(name=name, n_sim=n_sim, sim='stochkit')\n",
    "    cuda_file = \"{name}_{n_sim}_{sim}.csv.gz\".format(name=name, n_sim=n_sim, sim='cuda')\n",
    "    opencl_file = \"{name}_{n_sim}_{sim}.csv.gz\".format(name=name, n_sim=n_sim, sim='opencl')\n",
    "    \n",
    "    bng_data = pd.read_csv(bng_file)\n",
    "    stochkit_data = pd.read_csv(stochkit_file)\n",
    "    cuda_data = pd.read_csv(cuda_file)\n",
    "    opencl_data = pd.read_csv(opencl_file)\n",
    "    bng_data['simulator'] = 'bng'\n",
    "    stochkit_data['simulator'] = 'stochkit'\n",
    "    cuda_data['simulator'] = 'CUDASimulator'\n",
    "    opencl_data['simulator'] = 'OpenCLSimulator'\n",
    "\n",
    "    df = pd.concat([bng_data, stochkit_data, cuda_data, opencl_data], ignore_index=True)\n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_plot(df, obs):\n",
    "    plt.figure(figsize=(6,6))\n",
    "    ax = sns.lineplot(\n",
    "        x=\"time\", y=obs,\n",
    "        markers=True, dashes=True,\n",
    "        ci='sd', \n",
    "        estimator=np.average,\n",
    "        hue=\"simulator\", \n",
    "        data=df\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tspan = np.linspace(0, 100, 301)\n",
    "scholgl_model.parameters['X_0'].value = 250\n",
    "save_output(scholgl_model, tspan, n_sim=1000)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}