{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Power Spectra (Bos 2016)\n\nHere we calculate the power spectra of the :cite:t:`potjans2014` microcircuit\nmodel including modifications made in :cite:t:`bos2016`.\n\nThis example reproduces Fig. 1E in :cite:t:`bos2016`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nnmt\nimport numpy as np\nimport matplotlib as mpl\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport matplotlib.ticker\n\nplt.style.use('frontiers.mplstyle')\nmpl.rcParams.update({'legend.fontsize': 'medium',  # old: 5.0 was too small\n                     'axes.titlepad': 0.0})"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "First, create an instance of the network model class `Microcircuit`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "microcircuit = nnmt.models.Microcircuit(\n    network_params=\n    '../../tests/fixtures/integration/config/Bos2016_network_params.yaml',\n    analysis_params=\n    '../../tests/fixtures/integration/config/Bos2016_analysis_params.yaml')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The frequency resolution used in the original publication was quite high.\nHere, we reduce the frequency resolution for faster execution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "reduce_frequency_resolution = True\n\nif reduce_frequency_resolution:\n    microcircuit.change_parameters(changed_analysis_params={'df': 5},\n                                overwrite=True)\n    derived_analysis_params = (\n        microcircuit._calculate_dependent_analysis_parameters())\n    microcircuit.analysis_params.update(derived_analysis_params)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calculate all necessary quantities and finally the power spectra.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# calculate working point for exponentially shape post synaptic currents\nnnmt.lif.exp.working_point(microcircuit, method='taylor')\n# calculate the transfer function\nnnmt.lif.exp.transfer_function(microcircuit, method='taylor')\n# calculate the delay distribution matrix\nnnmt.network_properties.delay_dist_matrix(microcircuit)\n# calculate the effective connectivity matrix\nnnmt.lif.exp.effective_connectivity(microcircuit)\n# calculate the power spectra\npower_spectra = nnmt.lif.exp.power_spectra(microcircuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Read the simulated power spectra from the publicated data for comparison.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fix_path = '../../tests/fixtures/integration/data/'\nresult = nnmt.input_output.load_h5(fix_path + \n                                   'Bos2016_publicated_and_converted_data.h5')\nsimulated_power_spectra_1_window = result['fig_microcircuit']['1']\nsimulated_power_spectra_20_window = result['fig_microcircuit']['20']"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting mean-field prediction and simulated results together.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# two column figure, 180 mm wide\nwidth = 180. / 25.4 \nheight = 90. / 25.4\nfig = plt.figure(figsize=(width, height))\ngrid_specification = gridspec.GridSpec(2, 4, figure=fig)\n\nfor layer in [0, 1, 2, 3]:\n    for pop in [0, 1]:\n        j = layer*2+pop\n        \n        ax = fig.add_subplot(grid_specification[pop, layer])\n        \n        ax.plot(simulated_power_spectra_1_window['freq_sim'],\n                   simulated_power_spectra_1_window[f'power{j}'],\n                   color=(0.8, 0.8, 0.8),\n                   label='simulation')\n        ax.plot(simulated_power_spectra_20_window['freq_sim'],\n                   simulated_power_spectra_20_window[f'power{j}'], \n                   color=(0.5, 0.5, 0.5),\n                   label='simulation avg.')\n        ax.plot(microcircuit.analysis_params['omegas']/(2*np.pi),\n                   power_spectra[:, j],\n                   color='black', zorder=2,\n                   label='prediction')\n        \n        ax.set_xlim([10.0, 400.0])\n        ax.set_ylim([1e-6, 1e0])\n        ax.set_yscale('log')\n        \n        population_name = microcircuit.network_params['populations'][j] \n        if population_name == '23E':\n            ax.set_title('2/3E')\n        elif population_name == '23I':\n            ax.set_title('2/3I')\n        else:\n            ax.set_title(population_name)\n            \n            \n        if pop == 1 and layer == 0:\n            ax.set_xlabel(r'frequency $\\omega/2\\pi\\quad(1/\\mathrm{s})$')\n        else:\n            ax.set_xticklabels([])\n        \n        ax.set_xticks([100, 200, 300])\n        y_minor = matplotlib.ticker.LogLocator(\n            base = 10.0, \n            subs = np.arange(1.0, 10.0) * 0.1, \n            numticks = 10)\n        ax.yaxis.set_minor_locator(y_minor)\n        ax.yaxis.set_minor_formatter(matplotlib.ticker.NullFormatter())\n        ax.set_yticks([])\n        if j == 0 or j== 1:\n            ax.set_yticks([1e-5,1e-3,1e-1])\n            \n        if j == 0:\n            ax.set_ylabel(r'power spectrum $|C(\\omega)|\\quad(1/\\mathrm{s}^2)$')\n            ax.legend()\n    \nplt.savefig('figures/power_spectra_Bos2016.png')"
      ]
    }
  ],
  "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.9.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}