added demo 2 for network burst inference

Author: rdgao

notebooks/demo-1_automind_inference_workflow.ipynb

@@ -567,7 +567,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
    "outputs": [
     {
@@ -587,6 +587,10 @@
     ")\n",
     "df_samples.insert(loc=0, column=\"params_settings.batch_seed\", value=batch_seed)\n",
     "df_samples.insert(loc=1, column=\"params_settings.random_seed\", value=random_seeds)\n",
+    "\n",
+    "# NOTE: This sets the correct early-stopping condition, i.e., to assess the simulation from 0.1 to 10.1 seconds.\n",
+    "params_dict['params_analysis']['analysis_window'] = [0.1, None] \n",
+    "\n",
     "params_dict_run = data_utils.fill_params_dict(\n",
     "    params_dict, df_samples, posterior.as_dict, n_samples\n",
     ")"

notebooks/demo-2_automind_inference_from_spikes.ipynb

@@ -10,6 +10,8 @@
     "\n",
     "We will walk through how to compute network burst summary statistics such that we can use them for inference, get samples from the trained density estimator, and then simulate and analyze a few of those discovered model configurations to check if they indeed reproduce the target observation.\n",
     "\n",
+    "Many thanks to [Tereza Zuskinova](https://github.com/TerezaZuskinova) at IST Austria for her help in developing and testing this demo.\n",
+    "\n",
     "For more details, visit the [AutoMIND preprint](https://www.biorxiv.org/content/10.1101/2024.08.21.608969v1).\n",
     "\n",
     "---\n",
@@ -114,7 +116,9 @@
    "metadata": {},
    "source": [
     "### Load and process target data from raw spike times\n",
-    "Here, we use another recording from the organoid dataset, and demonstrate how to start from raw spike times to computing summary features amenable to be used for inference."
+    "Here, we use another recording from the organoid dataset, and demonstrate how to start from raw spike times to computing summary features amenable to be used for inference.\n",
+    "\n",
+    "You can download the `.npz` file [here](https://figshare.com/s/3f1467f8fb0f328aed16) in `./observations`."
    ]
   },
   {
@@ -164,28 +168,6 @@
     "spikes_in_list[0]"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": 5,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "\n",
-    "\n",
-    "# data_all = np.load(\"../../../data/MEA_muotri/trujillo/spikes_only.npz\", allow_pickle=True)\n",
-    "# idx = 26\n",
-    "# print(data_all[\"recs\"][idx], data_all[\"t\"][idx][-1])\n",
-    "# spikes_in_list = data_all[\"spikes\"][idx][0]\n",
-    "# t_end = 250.\n",
-    "# rec_name = data_all[\"recs\"][idx]\n",
-    "# # org_spikes_all, org_t_all, org_name_all = (\n",
-    "# #     data_all[\"spikes\"],\n",
-    "# #     data_all[\"t\"],\n",
-    "# #     data_all[\"recs\"],\n",
-    "# # )\n",
-    "# np.savez('./example_raw_data.npz', spikes=spikes_in_list, t_end=t_end, recording_name=rec_name, fs=12500)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -197,13 +179,12 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
    "outputs": [],
    "source": [
-    "spikes_dict = data_utils.convert_spike_array_to_dict(spikes_in_list, fs, 'exc_spikes') # \n",
-    "spikes_dict['t_end'] = t_end\n",
-    "spikes_dict['inh_spikes'] = {}\n",
+    "spikes_dict = data_utils.convert_spike_array_to_dict(spikes_in_list, fs, 'exc_spikes')\n",
+    "spikes_dict['t_end'] = t_end # end recording time, you can take this as the last spike time plus a little more.\n",
     "\n",
     "# Analyze data with the same params_dict as the simulated data. \n",
     "# Note that depending on your own data, some settings may need to be adjusted in params_dict['params_analysis'].\n",
@@ -221,7 +202,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 6,
    "metadata": {},
    "outputs": [
     {
@@ -238,7 +219,7 @@
    "source": [
     "xlims = [0, 250]\n",
     "fig, axs = plt.subplots(2,1,figsize=(2, 1))\n",
-    "plot_utils._plot_raster_pretty(spikes_dict, XL=xlims, every_other=50, fontsize=5, ax=axs[0])\n",
+    "plot_utils._plot_raster_pretty(spikes_dict, plot_inh=False, XL=xlims, every_other=50, fontsize=5, ax=axs[0])\n",
     "plot_utils._plot_rates_pretty(\n",
     "        spikes_analyzed[1],\n",
     "        XL=xlims,        \n",
@@ -267,7 +248,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
@@ -317,7 +298,7 @@
        "0                1.421244  "
       ]
      },
-     "execution_count": 8,
+     "execution_count": 7,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -329,7 +310,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 8,
    "metadata": {},
    "outputs": [
     {
@@ -369,13 +350,13 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 10,
+   "execution_count": 9,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "8c869864d4634528b092c3e95293da68",
+       "model_id": "332026e114db48129bb0a59970f5749b",
        "version_major": 2,
        "version_minor": 0
       },
@@ -1065,7 +1046,7 @@
        "[20 rows x 30 columns]"
       ]
      },
-     "execution_count": 10,
+     "execution_count": 9,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -1108,7 +1089,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 10,
    "metadata": {},
    "outputs": [
     {
@@ -1156,7 +1137,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 11,
    "metadata": {},
    "outputs": [
     {
@@ -1187,15 +1168,15 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 12,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
       "cache non-existent.\n",
-      "20241121-1242|20241121-85|20241121-1629|20241121-5|20241121-375|20241121-1401|20241121-1002|20241121-343|20241121-1975|20241121-1430|20241121-211|20241121-1429|20241121-1791|20241121-360|20241121-1479|20241121-565|20241121-794|20241121-1121|20241121-1938|20241121-1829|Simulations took 158.22 seconds.\n",
+      "20241121-211|20241121-1121|20241121-794|20241121-565|20241121-1938|20241121-1242|20241121-85|20241121-1429|20241121-1975|20241121-1829|20241121-5|20241121-1401|20241121-375|20241121-360|20241121-1629|20241121-343|20241121-1479|20241121-1791|20241121-1430|20241121-1002|Simulations took 158.26 seconds.\n",
       "cache non-existent.\n"
      ]
     }
@@ -1233,7 +1214,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 13,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1254,7 +1235,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 14,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -1268,7 +1249,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 15,
    "metadata": {},
    "outputs": [
     {
@@ -1395,7 +1376,8 @@
     "    )    \n",
     "    ax2.set_ylabel(\"Rate\", labelpad=-5)\n",
     "\n",
-    "    plot_utils._plot_raster_pretty(spikes_dict, XL=xlims, every_other=50, fontsize=5, ax=ax3)\n",
+    "    # plot real data next to it\n",
+    "    plot_utils._plot_raster_pretty(spikes_dict, XL=xlims, plot_inh=False, every_other=50, fontsize=5, ax=ax3)\n",
     "    plot_utils._plot_rates_pretty(\n",
     "            spikes_analyzed[1],\n",
     "            XL=xlims,        \n",