Figure-4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "ba424429-ee98-40b5-8d37-d17792c68801",
   "metadata": {},
   "source": [
    "# Figure 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f8406d2-9887-45c4-8eaf-4b948c2f1fbb",
   "metadata": {},
   "source": [
    "NOTE: If you're running this on **Google Colab**, then uncomment and run the next two cells."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ef9ab21-0a1f-43f4-94d8-fab54930e875",
   "metadata": {},
   "outputs": [],
   "source": [
    "# !git clone https://github.com/Mark-Kramer/Aperiodic-Exponent-Model.git"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "026b6f97-c105-4d75-bcb2-e695be43d902",
   "metadata": {},
   "outputs": [],
   "source": [
    "# import sys\n",
    "# sys.path.insert(0,'/content/Aperiodic-Exponent-Model')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3470c5ec-ba91-4c10-b3ff-58789f67e1fb",
   "metadata": {},
   "source": [
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a9d8914-6a7a-4064-bea9-2bf5f7ebe749",
   "metadata": {},
   "source": [
    "## Load packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b187321b-14f3-4cc1-ae68-c89346db50ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.io as io\n",
    "import matplotlib.pyplot as plt\n",
    "from   matplotlib.backends.backend_pdf import PdfPages"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f05d7bf-8269-474c-a01e-9fc1d2dfb2cc",
   "metadata": {},
   "source": [
    "## Load functions for computations and plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d83d19fc-2265-482e-bcd6-4780a0c5733c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from one_over_f_functions import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "751bd552-3820-4ec3-9e8e-3da838e4aded",
   "metadata": {},
   "source": [
    "## Define the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "50e822fe-67a0-451d-9585-501c5db9a07c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate the generative model.\n",
    "def predator_prey(noise,params):\n",
    "    \n",
    "    N = 250000\n",
    "    dt = 0.0002\n",
    "    alpha = params[\"alpha\"]\n",
    "    gamma = params[\"gamma\"]\n",
    "    \n",
    "    x = np.zeros([N])                         # Output variables\n",
    "    y = np.zeros([N])\n",
    "    t = np.arange(N)*dt                       # Time [a.u.]\n",
    "    \n",
    "    x[0] = 0.2+0.2*np.random.rand()           # Initial conditions\n",
    "    y[0] = 0.2+0.2*np.random.rand()\n",
    "    \n",
    "    for n in np.arange(N-1):                  # Simulate the model\n",
    "        x[n+1] = x[n] + dt*(x[n]*(gamma-x[n])/gamma - x[n]*y[n] + noise[\"x\"]*np.random.randn())\n",
    "        y[n+1] = y[n] + dt*(-alpha*y[n]             + x[n]*y[n] + noise[\"y\"]*np.random.randn())\n",
    "        \n",
    "    x = x[50000:]\n",
    "    y = y[50000:]\n",
    "    t = t[50000:]-t[50000]\n",
    "    \n",
    "    return x,y,t"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27943068-c189-420c-8abc-e6f40e69cb7f",
   "metadata": {},
   "source": [
    "## Figure 4A-D: Simulate the model, and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "875b8950-83cb-4a7c-aa42-726a6011b63e",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = {\"x\":[],                           # Set the prey noise,\n",
    "         \"y\":1  }                          # ... and the predator noise.\n",
    "params= {\"alpha\":0.25,                     # Set the parameters.\n",
    "         \"gamma\":0.6}\n",
    "fig, ((ax1, ax2,), (ax3, ax4)) = plt.subplots(2, 2,figsize=(12, 6), dpi=80)\n",
    "\n",
    "noise[\"x\"]=0.0;                            # Prey noise = 0.\n",
    "[x,y,t] = predator_prey(noise,params)      # Simulate model,\n",
    "[S,f]   = compute_spectrum(x,t)            # ... compute the spectrum & fit aperiodic exponent.\n",
    "[aperiodic_exponent, x_linear_fit, y_linear_fit] = estimate_aperiodic_exponent(S,f,[50,250])\n",
    "make_plots(t,x, f,S, x_linear_fit,y_linear_fit,aperiodic_exponent, ax1,ax3);\n",
    "ax1.set_ylim([0.2, 0.3]); ax1.set_xlabel('Time [a.u.]'); ax3.set_xlim([-1,3]); ax3.set_ylim([-25, 0]);\n",
    "\n",
    "noise[\"x\"]=0.01;                           # Prey noise > 0.\n",
    "[x,y,t] = predator_prey(noise,params)      # Simulate model,\n",
    "[S,f]   = compute_spectrum(x,t)            # ... compute the spectrum & fit aperiodic exponent.\n",
    "[aperiodic_exponent, x_linear_fit, y_linear_fit] = estimate_aperiodic_exponent(S,f,[50,250])\n",
    "make_plots(t,x, f,S, x_linear_fit,y_linear_fit,aperiodic_exponent, ax2,ax4);\n",
    "ax2.set_ylim([0.2, 0.3]); ax2.set_xlabel('Time [a.u.]'); ax4.set_xlim([-1,3]); ax4.set_ylim([-25, 0]);\n",
    "\n",
    "#fig.savefig(\"./PDFs/Figure-4ABCD.pdf\", bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3fb871a9-88c7-4f98-b93b-773f64902473",
   "metadata": {},
   "source": [
    "## Figure 4E (Part 1): Iterate over values of noise and plot aperiodic exponent\n",
    "This part is slow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6094ce4-2024-433b-910d-5c271ae70570",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = {\"x\":[],                                    # Prepare the noise variable.\n",
    "         \"y\":1  }                                   # ... fix the predator noise,\n",
    "noise_values = [0,1e-5,2e-5,3e-5,4e-5,5e-5,         # ... list the prey noises.\n",
    "                  1e-4,2e-4,3e-4,4e-4,5e-4,\n",
    "                  6e-4,7e-4,8e-4,9e-4,\n",
    "                  1e-3,1.5e-3,2e-3,3e-3,4e-3,5e-3]\n",
    "params = {\"alpha\":[], \"gamma\":[]}\n",
    "N_replicates = 10                                   # Number of times to repeat simulation.\n",
    "alpha_params = np.arange(0,  0.51,0.1)              # Iterate over parameter alpha,\n",
    "gamma_params = np.arange(0.1,1.01,0.1)              # ... and over parameter gamma.\n",
    "\n",
    "aperiodic_exponents    = np.empty((np.size(alpha_params), np.size(noise_values), np.size(gamma_params)*N_replicates))\n",
    "aperiodic_exponents[:] = np.nan\n",
    "\n",
    "for [j,alpha] in enumerate(alpha_params):           # For each alpha model parameter,\n",
    "    params[\"alpha\"] = alpha\n",
    "\n",
    "    for [n, prey_noise] in enumerate(noise_values): # ... and for each prey noise,\n",
    "        noise[\"x\"]=prey_noise\n",
    "    \n",
    "        counter = 0\n",
    "        for gamma in (gamma_params):                # ... sim model with different gamma parameters,\n",
    "            params[\"gamma\"] = gamma\n",
    "        \n",
    "            print(\"alpha=\",alpha, \", noise_x=\",prey_noise, \", gamma=\",gamma)\n",
    "        \n",
    "            if (1-alpha/gamma)>0:                   # Check physical condition for population > 0.\n",
    "        \n",
    "                print(\"Run it\")\n",
    "\n",
    "                for r in np.arange(N_replicates):   # ... and repeat sim N_replicate times,\n",
    "                    [x,y,t] = predator_prey(noise,params)\n",
    "                    [S,f]  = compute_spectrum(x,t)                     # ... compute the spectrum,\n",
    "                    ae = estimate_aperiodic_exponent(S,f,[50,250])[0]  # ... fit aperiodic exponent,\n",
    "                    aperiodic_exponents[j,n,counter] = ae              # ... and save it.\n",
    "                    counter=counter+1\n",
    "                                                    # Save the results.\n",
    "res = {\"alpha_params\":alpha_params, \"gamma_params\":gamma_params, \"noise_values\":noise_values, \"aperiodic_exponents\":aperiodic_exponents}\n",
    "io.savemat(str('Figure-4E.mat'), res)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "537da9c6-3818-4733-b275-5333e1292b23",
   "metadata": {},
   "source": [
    "## Figure 4E (Part 2): Load the aperiodic exponent, and plot it.\n",
    "\n",
    "This part is fast."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "162171be-5d0c-4502-9f51-916b28ae2d62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note: if you're running on Google Colab, then load the .mat file like this:\n",
    "# res = io.loadmat(\"/content/Aperiodic-Exponent-Model/Figure-4E.mat\");\n",
    "\n",
    "# Or, if not running on Google Colab, then load the .mat file like this:\n",
    "res = io.loadmat(\"Figure-4E.mat\")\n",
    "\n",
    "fig = plt.figure(figsize=(12, 4), dpi=80)\n",
    "alpha_params = res[\"alpha_params\"][0]\n",
    "colors = [\"red\", \"orange\", \"yellow\", \"green\", \"blue\", \"indigo\"]\n",
    "\n",
    "for [k,alpha] in enumerate(alpha_params):\n",
    "    print(\"alpha\", alpha, \"- color:\", colors[k])\n",
    "    res_at_a = {\"noise_values\":res[\"noise_values\"], \"aperiodic_exponents\":res[\"aperiodic_exponents\"][k,:],\n",
    "                \"color\":colors[k]};\n",
    "    make_plot_ae_vs_noise(res_at_a);\n",
    "plt.xlim([0,0.005]);\n",
    "\n",
    "#fig.savefig(\"./PDFs/Figure-4E.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
}