Result summary.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caa81a92",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pds\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0c22abb",
   "metadata": {},
   "outputs": [],
   "source": [
    "method_dict = {\"vi\": \"D-CODE\", \"diff\": \"SR-T\", \"spline\": \"SR-S\", \"gp\": \"SR-G\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "249e983a",
   "metadata": {},
   "outputs": [],
   "source": [
    "val_dict = {\n",
    "    \"noise\": \"sigma\",\n",
    "    \"freq\": \"del_t\",\n",
    "    \"n\": \"n\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8bb315d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_list = [\"GompertzODE\", \"LogisticODE\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4189eda7",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_df(df, x_val=\"sigma\"):\n",
    "    for method in method_dict.keys():\n",
    "        df_sub = df[df.method == method]\n",
    "        df_sub = df_sub.dropna()\n",
    "\n",
    "        plt.fill_between(\n",
    "            df_sub[x_val],\n",
    "            df_sub.rate - df_sub.rate_sd,\n",
    "            df_sub.rate + df_sub.rate_sd,\n",
    "            alpha=0.3,\n",
    "        )\n",
    "\n",
    "        plt.plot(df_sub[x_val], df_sub.rate, \"o-\", label=method_dict[method])\n",
    "        plt.ylim(-0.05, 1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "71d13c12",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "plt.style.use(\"tableau-colorblind10\")\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "plt.rcParams[\"font.size\"] = \"13\"\n",
    "\n",
    "counter = 1\n",
    "for i in range(len(ode_list)):\n",
    "    ode = ode_list[i]\n",
    "    for val_key, x_val in val_dict.items():\n",
    "        print(ode, val_key, x_val)\n",
    "        df = pds.read_csv(\"results/{}-{}.txt\".format(ode, val_key), header=None)\n",
    "\n",
    "        df.columns = [\n",
    "            \"ode\",\n",
    "            \"freq\",\n",
    "            \"n\",\n",
    "            \"sigma\",\n",
    "            \"method\",\n",
    "            \"rate\",\n",
    "            \"rate_sd\",\n",
    "            \"ks\",\n",
    "            \"ks_sd\",\n",
    "        ]\n",
    "        df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "        df = df.sort_values([\"method\", x_val])\n",
    "\n",
    "        plot_conf = 230 + counter\n",
    "        plt.subplot(plot_conf)\n",
    "        plot_df(df, x_val=x_val)\n",
    "\n",
    "        if counter == 1 or counter == 4:\n",
    "            plt.ylabel(\"Success Prob.\", size=16)\n",
    "\n",
    "        if counter == 1:\n",
    "            plt.title(\"Varying noise level $\\sigma_R$\")\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 2:\n",
    "            plt.title(\"Gompertz Model \\n Varying step size $\\Delta t$\")\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 3:\n",
    "            plt.title(\"Varying sample size $N$\")\n",
    "        elif counter == 5:\n",
    "            plt.title(\"Generalized Logistic Model\")\n",
    "\n",
    "        if counter == 4:\n",
    "            plt.xlabel(r\"$\\sigma_R$\", size=16)\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 5:\n",
    "            plt.xlabel(r\"$\\Delta t$\", size=16)\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 6:\n",
    "            plt.xlabel(r\"$N$\", size=16)\n",
    "\n",
    "        counter += 1\n",
    "\n",
    "\n",
    "plt.legend(title=\"Methods\", bbox_to_anchor=(1.05, 1), loc=\"upper left\")\n",
    "plt.tight_layout(pad=0.2)\n",
    "plt.savefig(\"growth_results.png\", dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bece3c8e",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "e612c3a4",
   "metadata": {},
   "source": [
    "## Selkov Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ce75e4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "ode = \"SelkovODE\"\n",
    "\n",
    "val_key = list(val_dict.keys())[0]\n",
    "x_val = val_dict[val_key]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "caddd892",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pds.read_csv(\"results/{}-{}.txt\".format(ode, val_key), header=None)\n",
    "\n",
    "df.columns = [\"ode\", \"freq\", \"n\", \"sigma\", \"method\", \"rate\", \"rate_sd\", \"ks\", \"ks_sd\"]\n",
    "df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "df = df.sort_values([\"method\", x_val])\n",
    "df[\"x_id\"] = 0\n",
    "df0 = df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "583c4bf9",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pds.read_csv(\"results/{}-{}-1.txt\".format(ode, val_key), header=None)\n",
    "\n",
    "df.columns = [\"ode\", \"freq\", \"n\", \"sigma\", \"method\", \"rate\", \"rate_sd\", \"ks\", \"ks_sd\"]\n",
    "df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "df = df.sort_values([\"method\", x_val])\n",
    "df[\"x_id\"] = 1\n",
    "df1 = df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9387e01e",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pds.read_csv(\"results/{}-{}-param.txt\".format(ode, val_key), header=None)\n",
    "\n",
    "df.columns = [\n",
    "    \"ode\",\n",
    "    \"freq\",\n",
    "    \"n\",\n",
    "    \"sigma\",\n",
    "    \"method\",\n",
    "    \"sigma_rmse\",\n",
    "    \"sigma_sd\",\n",
    "    \"rho_rmse\",\n",
    "    \"rho_sd\",\n",
    "]\n",
    "df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "df = df.sort_values([\"method\", x_val])\n",
    "df[\"x_id\"] = 0\n",
    "df0_param = df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09e14f44",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pds.read_csv(\"results/{}-{}-param-1.txt\".format(ode, val_key), header=None)\n",
    "\n",
    "df.columns = [\n",
    "    \"ode\",\n",
    "    \"freq\",\n",
    "    \"n\",\n",
    "    \"sigma\",\n",
    "    \"method\",\n",
    "    \"sigma_rmse\",\n",
    "    \"sigma_sd\",\n",
    "    \"rho_rmse\",\n",
    "    \"rho_sd\",\n",
    "]\n",
    "df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "df = df.sort_values([\"method\", x_val])\n",
    "df[\"x_id\"] = 1\n",
    "df1_param = df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "956031dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pds.concat([df0, df1])\n",
    "df_param = pds.concat([df0_param, df1_param])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bbca74f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pds.merge(df, df_param)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9d4ebc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_rate = pds.pivot_table(df, values=\"rate\", index=[\"x_id\", \"method\"], columns=\"sigma\")\n",
    "tbl_rate[\"val\"] = \"rate\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ed70f299",
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_rate_sd = pds.pivot_table(\n",
    "    df, values=\"rate_sd\", index=[\"x_id\", \"method\"], columns=\"sigma\"\n",
    ")\n",
    "tbl_rate_sd[\"val\"] = \"rate_sd\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "74117c39",
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_sigma_rmse = pds.pivot_table(\n",
    "    df, values=\"sigma_rmse\", index=[\"x_id\", \"method\"], columns=\"sigma\"\n",
    ")\n",
    "tbl_sigma_rmse[\"val\"] = \"sigma_rmse\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "039a1ec8",
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_sigma_sd = pds.pivot_table(\n",
    "    df, values=\"sigma_sd\", index=[\"x_id\", \"method\"], columns=\"sigma\"\n",
    ")\n",
    "tbl_sigma_sd[\"val\"] = \"sigma_sd\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57f0286e",
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_sigma_ks = pds.pivot_table(\n",
    "    df, values=\"ks\", index=[\"x_id\", \"method\"], columns=\"sigma\"\n",
    ")\n",
    "tbl_sigma_ks[\"val\"] = \"ks\"\n",
    "tbl_sigma_ks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6887512",
   "metadata": {},
   "outputs": [],
   "source": [
    "tbl_sigma_ks_sd = pds.pivot_table(\n",
    "    df, values=\"ks_sd\", index=[\"x_id\", \"method\"], columns=\"sigma\"\n",
    ")\n",
    "tbl_sigma_ks_sd[\"val\"] = \"ks_sd\"\n",
    "tbl_sigma_ks_sd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "747d235b",
   "metadata": {},
   "outputs": [],
   "source": [
    "selkov_table = pds.concat([tbl_rate, tbl_rate_sd, tbl_sigma_rmse, tbl_sigma_sd])\n",
    "selkov_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d47fa0dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "tt = selkov_table.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1ae0f5f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "tt[tt[\"method\"] == \"gp\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2d3e6ab",
   "metadata": {},
   "outputs": [],
   "source": [
    "tt[tt[\"method\"] == \"diff\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0846ec69",
   "metadata": {},
   "outputs": [],
   "source": [
    "selkov_table.to_csv(\"Selkov_results.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5138f8b6",
   "metadata": {},
   "source": [
    "## Lorenz results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "771376ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_df_ax(df, ax, x_val=\"sigma\"):\n",
    "    for method in method_dict.keys():\n",
    "        df_sub = df[df.method == method]\n",
    "        df_sub = df_sub.dropna()\n",
    "\n",
    "        ax.fill_between(\n",
    "            df_sub[x_val],\n",
    "            df_sub.rate - df_sub.rate_sd,\n",
    "            df_sub.rate + df_sub.rate_sd,\n",
    "            alpha=0.3,\n",
    "        )\n",
    "\n",
    "        ax.plot(df_sub[x_val], df_sub.rate, \"o-\", label=method_dict[method])\n",
    "        ax.set_ylim(-0.05, 1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "261120db",
   "metadata": {},
   "outputs": [],
   "source": [
    "## this part is the placeholder\n",
    "\n",
    "\n",
    "def lorenz(x, y, z, s=10, r=28, b=2.667):\n",
    "    \"\"\"\n",
    "    Given:\n",
    "       x, y, z: a point of interest in three dimensional space\n",
    "       s, r, b: parameters defining the lorenz attractor\n",
    "    Returns:\n",
    "       x_dot, y_dot, z_dot: values of the lorenz attractor's partial\n",
    "           derivatives at the point x, y, z\n",
    "    \"\"\"\n",
    "    x_dot = s * (y - x)\n",
    "    y_dot = r * x - y - x * z\n",
    "    z_dot = x * y - b * z\n",
    "    return x_dot, y_dot, z_dot\n",
    "\n",
    "\n",
    "dt = 0.01\n",
    "num_steps = 10000\n",
    "\n",
    "# Need one more for the initial values\n",
    "xs = np.empty(num_steps + 1)\n",
    "ys = np.empty(num_steps + 1)\n",
    "zs = np.empty(num_steps + 1)\n",
    "\n",
    "# Set initial values\n",
    "xs[0], ys[0], zs[0] = (0.0, 1.0, 1.05)\n",
    "\n",
    "# Step through \"time\", calculating the partial derivatives at the current point\n",
    "# and using them to estimate the next point\n",
    "for i in range(num_steps):\n",
    "    x_dot, y_dot, z_dot = lorenz(xs[i], ys[i], zs[i])\n",
    "    xs[i + 1] = xs[i] + (x_dot * dt)\n",
    "    ys[i + 1] = ys[i] + (y_dot * dt)\n",
    "    zs[i + 1] = zs[i] + (z_dot * dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "07a7c113",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pickle\n",
    "\n",
    "with open(\"results/Lorenz_traj.pkl\", \"rb\") as f:\n",
    "    diff_dict = pickle.load(f)\n",
    "\n",
    "with open(\"results/Lorenz_vi_traj.pkl\", \"rb\") as f:\n",
    "    vi_dict = pickle.load(f)\n",
    "\n",
    "with open(\"results/Lorenz_true_traj.pkl\", \"rb\") as f:\n",
    "    true_dict = pickle.load(f)\n",
    "\n",
    "with open(\"results/Lorenz_node_traj2.pkl\", \"rb\") as f:\n",
    "    node_dict = pickle.load(f)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "72d6bf98",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_trac(ax, xs, ys, zs, title, lw=0.5):\n",
    "    elev = 5.0\n",
    "    azim = 120.0\n",
    "    ax.view_init(elev, azim)\n",
    "    ax.w_xaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "    ax.w_yaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "    ax.w_zaxis.set_pane_color((1.0, 1.0, 1.0, 1.0))\n",
    "\n",
    "    ax.plot(xs, ys, zs, lw=lw)\n",
    "    ax.set_title(title)\n",
    "\n",
    "\n",
    "#     ax.set_xticks([])\n",
    "#     ax.set_yticks([])\n",
    "#     ax.set_zticks([])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "482cf86d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "\n",
    "plt.figure(figsize=(12, 6))\n",
    "plt.style.use(\"tableau-colorblind10\")\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "plt.rcParams[\"font.size\"] = \"13\"\n",
    "\n",
    "\n",
    "gs = gridspec.GridSpec(3, 12)\n",
    "ax1a = plt.subplot(gs[0, :4])\n",
    "ax1b = plt.subplot(gs[0, 4:8])\n",
    "ax1c = plt.subplot(gs[0, 8:])\n",
    "\n",
    "# ax1a = plt.subplot(gs[0, :3])\n",
    "# ax1b = plt.subplot(gs[0, 3:6])\n",
    "# ax1c = plt.subplot(gs[0, 6:9])\n",
    "\n",
    "ax2a = plt.subplot(gs[1:, :3], projection=\"3d\")\n",
    "ax2b = plt.subplot(gs[1:, 3:6], projection=\"3d\")\n",
    "ax2c = plt.subplot(gs[1:, 6:9], projection=\"3d\")\n",
    "ax2d = plt.subplot(gs[1:, 9:], projection=\"3d\")\n",
    "\n",
    "\n",
    "for i, ax in enumerate(plt.gcf().axes):\n",
    "\n",
    "    if i < 3:\n",
    "        x_id = i\n",
    "\n",
    "        if x_id == 0:\n",
    "            df = pds.read_csv(\"results/Lorenz-noise.txt\", header=None)\n",
    "        else:\n",
    "            df = pds.read_csv(\"results/Lorenz-noise-{}.txt\".format(x_id), header=None)\n",
    "\n",
    "        df.columns = [\n",
    "            \"ode\",\n",
    "            \"freq\",\n",
    "            \"n\",\n",
    "            \"sigma\",\n",
    "            \"method\",\n",
    "            \"rate\",\n",
    "            \"rate_sd\",\n",
    "            \"ks\",\n",
    "            \"ks_sd\",\n",
    "        ]\n",
    "        df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "        df = df.sort_values([\"method\", \"sigma\"])\n",
    "        plot_df_ax(df, ax)\n",
    "        ax.set_xlabel(\"Noise $\\sigma_R$\")\n",
    "        if i == 0:\n",
    "            ax.set_title(\"Success Prob. $\\dot{x}_1(t)$\")\n",
    "        elif i == 1:\n",
    "            ax.set_title(\"Success Prob. $\\dot{x}_2(t)$\")\n",
    "        else:\n",
    "            ax.set_title(\"Success Prob. $\\dot{x}_3(t)$\")\n",
    "\n",
    "            ax.legend(bbox_to_anchor=(1.005, 1), loc=\"upper left\", fontsize=10)\n",
    "    #             ax.legend(loc='center left', fontsize=10)\n",
    "    else:\n",
    "\n",
    "        if i == 3:\n",
    "            plot_trac(ax, true_dict[\"x\"], true_dict[\"y\"], true_dict[\"z\"], str(i))\n",
    "            ax.set_title(\"Ground truth\")\n",
    "        elif i == 4:\n",
    "            plot_trac(ax, vi_dict[\"x\"], vi_dict[\"y\"], vi_dict[\"z\"], str(i), lw=0.3)\n",
    "            ax.set_title(\"D-CODE\")\n",
    "        elif i == 5:\n",
    "            plot_trac(\n",
    "                ax, diff_dict[\"x\"], diff_dict[\"y\"], diff_dict[\"z\"], str(i), lw=0.8\n",
    "            )\n",
    "            ax.set_title(\"SR-T\")\n",
    "        else:\n",
    "            plot_trac(\n",
    "                ax, node_dict[\"x\"], node_dict[\"y\"], node_dict[\"z\"], str(i), lw=2.0\n",
    "            )\n",
    "            ax.set_title(\"Neural ODE\")\n",
    "\n",
    "        ax.set_zlim(0, 50)\n",
    "        ax.set_xlim(-25, 25)\n",
    "        ax.set_ylim(-25, 25)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "plt.savefig(\"lorenz.png\", dpi=200)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "998fba0d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "db6539c9",
   "metadata": {},
   "source": [
    "## sensitivity plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8e6405d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_df2(df_sub, x_val=\"n_basis\"):\n",
    "\n",
    "    plt.fill_between(\n",
    "        df_sub[x_val],\n",
    "        df_sub.rate - df_sub.rate_sd,\n",
    "        df_sub.rate + df_sub.rate_sd,\n",
    "        alpha=0.3,\n",
    "    )\n",
    "\n",
    "    plt.plot(df_sub[x_val], df_sub.rate, \"o-\")\n",
    "    plt.ylim(-0.05, 1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44469395",
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_list = [\"GompertzODE\", \"Lorenz\"]\n",
    "bas_list = [\"sine\", \"cubic\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c81d0d92",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "plt.style.use(\"tableau-colorblind10\")\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "plt.rcParams[\"font.size\"] = \"13\"\n",
    "\n",
    "counter = 1\n",
    "for i in range(len(ode_list)):\n",
    "    ode = ode_list[i]\n",
    "\n",
    "    for bas in bas_list:\n",
    "\n",
    "        df = pds.read_csv(\"results/sensitivity_{}.txt\".format(ode), header=None)\n",
    "        df.columns = [\"ode\", \"basis\", \"n_basis\", \"N\", \"rate\", \"rate_sd\", \"ks\", \"ks_sd\"]\n",
    "        df = df.sort_values([\"basis\", \"n_basis\"])\n",
    "\n",
    "        df_sub = df[df[\"basis\"] == bas]\n",
    "        df_sub = df_sub.dropna()\n",
    "\n",
    "        plot_conf = 220 + counter\n",
    "        plt.subplot(plot_conf)\n",
    "        plot_df2(df_sub)\n",
    "\n",
    "        if counter > 2:\n",
    "            plt.xlabel(\"Number of basis\", size=16)\n",
    "\n",
    "        #         if counter == 1 or counter == 4:\n",
    "        #             plt.ylabel('Recovery Rate', size=16)\n",
    "        plt.title(\"{} - {}\".format(ode, bas))\n",
    "\n",
    "        #         if counter == 1:\n",
    "        #             plt.title('Varying noise level $\\sigma_R$')\n",
    "        #             plt.xscale('log')\n",
    "        #         elif counter == 2:\n",
    "        #             plt.title('Gompertz Model \\n Varying step size $\\Delta t$')\n",
    "        #         elif counter == 3:\n",
    "        #             plt.title('Varying sample size $N$')\n",
    "        #         elif counter == 4:\n",
    "        #             plt.title('Generalized Logistic Model')\n",
    "\n",
    "        #         if counter == 4:\n",
    "        #             plt.xlabel(r'$\\sigma_R$', size=16)\n",
    "        #             plt.xscale('log')\n",
    "        #         elif counter == 5:\n",
    "        #             plt.xlabel(r'$\\Delta t$', size=16)\n",
    "        #         elif counter == 6:\n",
    "        #             plt.xlabel(r'$N$', size=16)\n",
    "\n",
    "        counter += 1\n",
    "\n",
    "\n",
    "plt.tight_layout(pad=0.2)\n",
    "plt.savefig(\"sensitivity_results.png\", dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3c2c746a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e113185a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1c5a71fb",
   "metadata": {},
   "source": [
    "## objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb55c8eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "method_dict = {\"vi\": \"D-CODE\", \"diff\": \"SR-T\", \"spline\": \"SR-S\", \"gp\": \"SR-G\"}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ee84c5d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "val_dict = {\n",
    "    \"noise\": \"sigma\",\n",
    "    \"freq\": \"del_t\",\n",
    "    \"n\": \"n\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "392d4e77",
   "metadata": {},
   "outputs": [],
   "source": [
    "ode_list = [\"GompertzODE\", \"LogisticODE\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "65833b2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_df(df, x_val=\"sigma\"):\n",
    "    for method in method_dict.keys():\n",
    "        df_sub = df[df.method == method]\n",
    "        df_sub = df_sub.dropna()\n",
    "\n",
    "        #         if x_val == 'sigma':\n",
    "        #             df_sub = df_sub[df_sub[x_val] < 0.6]\n",
    "\n",
    "        plt.fill_between(\n",
    "            df_sub[x_val], df_sub.ks - df_sub.ks_sd, df_sub.ks + df_sub.ks_sd, alpha=0.3\n",
    "        )\n",
    "\n",
    "        plt.plot(df_sub[x_val], df_sub.ks, \"o-\", label=method_dict[method])\n",
    "\n",
    "\n",
    "#         plt.ylim([-0.05, None])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1b2fab75",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12, 4))\n",
    "plt.style.use(\"tableau-colorblind10\")\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "plt.rcParams[\"font.size\"] = \"13\"\n",
    "\n",
    "counter = 1\n",
    "for i in range(len(ode_list)):\n",
    "    ode = ode_list[i]\n",
    "    for val_key, x_val in val_dict.items():\n",
    "        df = pds.read_csv(\"results/{}-{}.txt\".format(ode, val_key), header=None)\n",
    "\n",
    "        df.columns = [\n",
    "            \"ode\",\n",
    "            \"freq\",\n",
    "            \"n\",\n",
    "            \"sigma\",\n",
    "            \"method\",\n",
    "            \"rate\",\n",
    "            \"rate_sd\",\n",
    "            \"ks\",\n",
    "            \"ks_sd\",\n",
    "        ]\n",
    "        df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "        df = df.sort_values([\"method\", x_val])\n",
    "\n",
    "        plot_conf = 230 + counter\n",
    "        plt.subplot(plot_conf)\n",
    "        plot_df(df, x_val=x_val)\n",
    "\n",
    "        if counter == 1 or counter == 4:\n",
    "            plt.ylabel(\"Objective $d_x$\", size=16)\n",
    "\n",
    "        if counter == 1:\n",
    "            plt.title(\"Varying noise level $\\sigma_R$\")\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 2:\n",
    "            plt.title(\"Gompertz Model \\n Varying step size $\\Delta t$\")\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 3:\n",
    "            plt.title(\"Varying sample size $N$\")\n",
    "        elif counter == 5:\n",
    "            plt.title(\"Generalized Logistic Model\")\n",
    "\n",
    "        if counter == 4:\n",
    "            plt.xlabel(r\"$\\sigma_R$\", size=16)\n",
    "            plt.xscale(\"log\")\n",
    "        elif counter == 5:\n",
    "            plt.xlabel(r\"$\\Delta t$\", size=16)\n",
    "        elif counter == 6:\n",
    "            plt.xlabel(r\"$N$\", size=16)\n",
    "\n",
    "        counter += 1\n",
    "\n",
    "\n",
    "plt.legend(title=\"Methods\", bbox_to_anchor=(1.05, 1), loc=\"upper left\")\n",
    "plt.tight_layout(pad=0.2)\n",
    "plt.savefig(\"growth_results_obj.png\", dpi=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c43a57a8",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b49955c",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "38687cbc",
   "metadata": {},
   "source": [
    "## Lorenz objective"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8158a255",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_df_ax2(df, ax, x_val=\"sigma\"):\n",
    "    for method in method_dict.keys():\n",
    "        df_sub = df[df.method == method]\n",
    "        df_sub = df_sub.dropna()\n",
    "\n",
    "        ax.fill_between(\n",
    "            df_sub[x_val], df_sub.ks - df_sub.ks_sd, df_sub.ks + df_sub.ks_sd, alpha=0.3\n",
    "        )\n",
    "\n",
    "        ax.plot(df_sub[x_val], df_sub.ks, \"o-\", label=method_dict[method])\n",
    "\n",
    "\n",
    "#         ax.set_ylim(-0.05, 1.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bc28f0d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "\n",
    "\n",
    "plt.figure(figsize=(12, 2.5))\n",
    "plt.style.use(\"tableau-colorblind10\")\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "plt.rcParams[\"font.size\"] = \"13\"\n",
    "\n",
    "\n",
    "# gs = gridspec.GridSpec(3, 9)\n",
    "ax1a = plt.subplot(1, 3, 1)\n",
    "ax1b = plt.subplot(1, 3, 2)\n",
    "ax1c = plt.subplot(1, 3, 3)\n",
    "\n",
    "# ax1a = plt.subplot(gs[0, :3])\n",
    "# ax1b = plt.subplot(gs[0, 3:6])\n",
    "# ax1c = plt.subplot(gs[0, 6:9])\n",
    "\n",
    "# ax2a = plt.subplot(gs[1:, :3], projection='3d')\n",
    "# ax2b = plt.subplot(gs[1:, 3:6], projection='3d')\n",
    "# ax2c = plt.subplot(gs[1:, 6:9], projection='3d')\n",
    "# ax2d = plt.subplot(gs[1:, 9:], projection='3d')\n",
    "\n",
    "\n",
    "for i, ax in enumerate(plt.gcf().axes):\n",
    "\n",
    "    if i < 3:\n",
    "        x_id = i\n",
    "\n",
    "        if x_id == 0:\n",
    "            df = pds.read_csv(\"results/Lorenz-noise.txt\", header=None)\n",
    "        else:\n",
    "            df = pds.read_csv(\"results/Lorenz-noise-{}.txt\".format(x_id), header=None)\n",
    "\n",
    "        df.columns = [\n",
    "            \"ode\",\n",
    "            \"freq\",\n",
    "            \"n\",\n",
    "            \"sigma\",\n",
    "            \"method\",\n",
    "            \"rate\",\n",
    "            \"rate_sd\",\n",
    "            \"ks\",\n",
    "            \"ks_sd\",\n",
    "        ]\n",
    "        df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "        df = df.sort_values([\"method\", \"sigma\"])\n",
    "        plot_df_ax2(df, ax)\n",
    "        ax.set_xlabel(\"Noise $\\sigma_R$\")\n",
    "        if i == 0:\n",
    "            ax.set_title(\"Objective $d_x$ for $\\dot{x}_1(t)$\")\n",
    "        elif i == 1:\n",
    "            ax.set_title(\"Objective $d_x$ for $\\dot{x}_2(t)$\")\n",
    "        else:\n",
    "            ax.set_title(\"Objective $d_x$ for $\\dot{x}_3(t)$\")\n",
    "\n",
    "            ax.legend(bbox_to_anchor=(1.005, 1), loc=\"upper left\", fontsize=10)\n",
    "#             ax.legend(loc='center left', fontsize=10)\n",
    "\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "\n",
    "plt.savefig(\"lorenz_objective.png\", dpi=200)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "364f0656",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "265f0dbf",
   "metadata": {},
   "source": [
    "## Fraction ODE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "83d63650",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "plt.figure(figsize=(12, 2.5))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.style.use(\"tableau-colorblind10\")\n",
    "colors = plt.rcParams[\"axes.prop_cycle\"].by_key()[\"color\"]\n",
    "plt.rcParams[\"font.size\"] = \"13\"\n",
    "df = pds.read_csv(\"results/FracODE-noise.txt\", header=None)\n",
    "df.columns = [\"ode\", \"freq\", \"n\", \"sigma\", \"method\", \"rate\", \"rate_sd\", \"ks\", \"ks_sd\"]\n",
    "df[\"del_t\"] = 1.0 / df[\"freq\"]\n",
    "df = df.sort_values([\"method\", \"sigma\"])\n",
    "\n",
    "x_val = \"sigma\"\n",
    "for method in method_dict.keys():\n",
    "    df_sub = df[df.method == method]\n",
    "    df_sub = df_sub.dropna()\n",
    "\n",
    "    plt.fill_between(\n",
    "        df_sub[x_val],\n",
    "        df_sub.rate - df_sub.rate_sd,\n",
    "        df_sub.rate + df_sub.rate_sd,\n",
    "        alpha=0.3,\n",
    "    )\n",
    "\n",
    "    plt.plot(df_sub[x_val], df_sub.rate, \"o-\", label=method_dict[method])\n",
    "    plt.ylim(-0.05, 1.05)\n",
    "plt.title(\"Discover Prob.\")\n",
    "plt.xlabel(\"Noise level $\\sigma$\")\n",
    "\n",
    "ax = plt.subplot(1, 2, 2)\n",
    "plot_df_ax2(df, ax)\n",
    "ax.set_title(\"Objective $d_x$\")\n",
    "ax.legend(bbox_to_anchor=(1.005, 1), loc=\"upper left\", fontsize=10)\n",
    "plt.xlabel(\"Noise level $\\sigma$\")\n",
    "plt.savefig(\"frac.png\", dpi=200)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc0100b3",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9a66987a",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de74ac31",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}