Figure-4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "98b14a52-4e2e-45cd-a988-8b2b17004e4a",
   "metadata": {},
   "source": [
    "# Figure 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9434ac2-45f4-4118-bdf2-5c43e2cccdcb",
   "metadata": {},
   "source": [
    "---\n",
    "If you're running this on **Google Colab**, then uncomment and run the next two cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f474dabe-2362-43d8-8f8d-5631b5b94b5b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !git clone https://github.com/Mark-Kramer/Golden-Framework.git"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7028880a-efc3-42b0-82b8-6f4b2af2f77e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import sys\n",
    "# sys.path.insert(0,'/content/Golden-Framework')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ad1eab1-eeed-4425-993d-a4ea3a858241",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a874dd6-6fcf-46bf-8166-cc313d760c9b",
   "metadata": {},
   "source": [
    "## Load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cbe0b839-03fb-492c-b150-b2f167bf27eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.io as io\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.backends.backend_pdf import PdfPages"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73292613-8edb-4ab4-921d-139fef77d231",
   "metadata": {},
   "source": [
    "## Load functions for simulations and plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85b528db-f366-47cd-9443-0453c39f0070",
   "metadata": {},
   "outputs": [],
   "source": [
    "from Golden_Ratio_functions import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "229fb3d8-c492-4b7c-a18b-93fa155ca414",
   "metadata": {},
   "source": [
    "## Figure 4B: Two ensembles & constant gain only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1125a7d0-66b1-4337-8857-2e35724237fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio   = (1+np.sqrt(5))/2\n",
    "f_k     = ratio**[4,4,4,4, 6,6,6,6]\n",
    "network = {\"f_k\":f_k, \"k_perturbed\":4}\n",
    "gain    = {\"f_S\":0, \"g_S\":0}\n",
    "\n",
    "x,t = run_damped_coupled_harmonic_oscillators(network, gain)\n",
    "\n",
    "f = plt.figure(figsize=(2, 3), dpi=80)\n",
    "plot_model_traces(t,x,network,gain)\n",
    "#f.savefig(\"./PDFs/Figure-4B.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0aeac1a9-c849-478d-8dbd-e9570351b648",
   "metadata": {},
   "source": [
    "## Figure 4C: Two ensembles & sinusoidal gain"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a507d350-0523-4714-997d-533c2b046cdd",
   "metadata": {},
   "outputs": [],
   "source": [
    "ratio   = (1+np.sqrt(5))/2\n",
    "f_k     = ratio**[4,4,4,4, 6,6,6,6]\n",
    "network = {\"f_k\":f_k, \"k_perturbed\":4}\n",
    "gain    = {\"f_S\":ratio**5, \"g_S\":50}\n",
    "\n",
    "x,t = run_damped_coupled_harmonic_oscillators(network, gain)\n",
    "\n",
    "f = plt.figure(figsize=(2, 3), dpi=80)\n",
    "plot_model_traces(t,x,network,gain)\n",
    "#f.savefig(\"./PDFs/Figure-4C.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15e30ec1-0283-4356-a3f0-07d7f985ae6d",
   "metadata": {},
   "source": [
    "## Figure 4D (Part 1): Compute target response versus sinusoidal gain frequency, save it.\n",
    "This part is slow, and iterates over a range of `f_S`. Results saved in a `.mat` file, available on GitHub."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f3fc9c8-780e-43e4-b46c-5f8c71b4d21a",
   "metadata": {},
   "outputs": [],
   "source": [
    "golden    = (1+np.sqrt(5))/2;\n",
    "f_k       = ratio**[4,4,4,4, 6,6,6,6]\n",
    "network   = {\"f_k\":f_k, \"k_perturbed\":4}\n",
    "\n",
    "f_S       = np.arange(1,50,0.1)\n",
    "g_S       = 50\n",
    "A         = np.zeros([16,f_S.size])\n",
    "\n",
    "for i in np.arange(f_S.size):\n",
    "\n",
    "    print([\"%.1f\"% f_S[i]][0],end=' ')\n",
    "    \n",
    "    gain    = {\"f_S\":f_S[i], \"g_S\":g_S}\n",
    "    x,t = run_damped_coupled_harmonic_oscillators(network, gain)\n",
    "    \n",
    "    amp   = np.abs(hilbert(x));\n",
    "    i0    = np.where( (t>=0) & (t<=1.5) );\n",
    "    A[:,i]= np.mean(amp[:,i0[0]],1)\n",
    "\n",
    "res = {\"A\":A, \"f_S\":f_S, \"g_S\":g_S, \"f_k\":f_k, \"network\":network, \"gain\":gain}\n",
    "io.savemat('Figure-4D.mat', res)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aa482f6b-a781-4395-ad75-8444aedf2964",
   "metadata": {},
   "source": [
    "## Figure 4D (Part 2): Load target response versus sinusoidal gain frequency, plot it.\n",
    "\n",
    "This part is fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "85710591-15b5-4d87-a3de-777428b45dbc",
   "metadata": {},
   "outputs": [],
   "source": [
    "res       = io.loadmat(\"Figure-4D.mat\")\n",
    "# Note: if you're running on Google Colab, then load the .mat file like this:\n",
    "# res = io.loadmat(\"/content/Golden-Framework/Figure-4D.mat\");\n",
    "f         = plt.figure(figsize=(5, 3), dpi=80)\n",
    "ratio     = (1+np.sqrt(5))/2\n",
    "f_k_ticks = ratio**np.arange(2,9+1)\n",
    "plot_gain_traces_two_groups(res, f_k_ticks)\n",
    "#f.savefig(\"./PDFs/Figure-4D.pdf\", bbox_inches='tight')"
   ]
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}