{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Power Spectra of Sub-Circuits (Bos 2016)\n\nHere we calculate the power spectra of subcircuits of the :cite:t:`potjans2014` \nmicrocircuit model including modifications made in :cite:t:`bos2016`.\n\nThis example reproduces Fig. 9 in :cite:t:`bos2016`.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import nnmt\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport matplotlib.gridspec as gridspec\nimport matplotlib.ticker\n\nplt.style.use('frontiers.mplstyle')"
      ]
    },
    {
      "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)\n\nfrequencies = microcircuit.analysis_params['omegas']/(2*np.pi)\nfull_indegree_matrix = microcircuit.network_params['K']"
      ]
    },
    {
      "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')"
      ]
    },
    {
      "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": [
        "Initialize a new network instance for the reduced circuit explaining the \n64 Hz oscillation. First calculate the working point with the full circuit, \nthen reduce the connectivity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "low_gamma_subcircuit = microcircuit.copy()\n\n# construct a matrix with the wanted connections\n# \"The circuit is composed of the connections corresponding to the five largest \n# matrix elements of Z amp (64 Hz) and the eight largest elements \n# of Z freq (64 Hz)\" (Bos 2016)\nreducing_matrix = np.zeros((8,8))\nreducing_matrix[0,0:4] = 1\nreducing_matrix[1,0:2] = 1\nreducing_matrix[2:4,2:4] = 1\nreducing_matrix[3,0] = 1\n\nlow_gamma_subcircuit.network_params.update({\n    'K': full_indegree_matrix*reducing_matrix\n})\n\n# calculate the transfer function\nnnmt.lif.exp.transfer_function(low_gamma_subcircuit, method='taylor')\n# calculate the delay distribution matrix\nnnmt.network_properties.delay_dist_matrix(low_gamma_subcircuit)\n# calculate the effective connectivity matrix\nnnmt.lif.exp.effective_connectivity(low_gamma_subcircuit)\n# calculate the power spectra\nlow_gamma_subcircuit_power_spectra = nnmt.lif.exp.power_spectra(\n    low_gamma_subcircuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Initialize a new network instance to calculate results without connections\nfrom 23E to 4I. First, calculate the working point with the full circuit, \nthen reduce the connectivity.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "without_23E_4I = microcircuit.copy()\n\n# construct a matrix with the wanted connections\nreducing_matrix = np.ones((8,8))\nreducing_matrix[3,0] = 0\n\nwithout_23E_4I.network_params.update({\n    'K': full_indegree_matrix*reducing_matrix\n})\n\n# calculate the transfer function\nnnmt.lif.exp.transfer_function(without_23E_4I, method='taylor')\n# calculate the delay distribution matrix\nnnmt.network_properties.delay_dist_matrix(without_23E_4I)\n# calculate the effective connectivity matrix\nnnmt.lif.exp.effective_connectivity(without_23E_4I)\n# calculate the power spectra\nwithout_23E_4I_power_spectra = nnmt.lif.exp.power_spectra(without_23E_4I)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Plotting everything together.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# two column figure, 180 mm wide\nfig = plt.figure(figsize=(7.08661, 7.08661/2),\n                 constrained_layout=True)\ngrid_specification = gridspec.GridSpec(2, 1, figure=fig)\ncolormap = 'coolwarm'\nlabels = ['23E', '23I', '4E', '4I', '5E', '5I', '6E', '6I']\n\n\ngsA = gridspec.GridSpecFromSubplotSpec(1, 3, \n                                        width_ratios=[1, 1, 1], \n                                        subplot_spec=grid_specification[0])\n\ngsB = gridspec.GridSpecFromSubplotSpec(1, 3, \n                                        width_ratios=[1, 1, 1], \n                                        subplot_spec=grid_specification[1])\n\n\n# reduced circuit for 64 Hz oscillation\nfor gs, pop_idx in zip(gsA, [0,2,3]):\n    ax = fig.add_subplot(gs)\n\n    ax.plot(result['oscillation_origin']['A']['freq_sim'],\n                result['oscillation_origin']['A'][f'power{pop_idx}_sim'], \n                color=(0.8, 0.8, 0.8))\n    ax.plot(result['oscillation_origin']['A']['freq_sim_av'],\n                result['oscillation_origin']['A'][f'power{pop_idx}_sim_av'], \n                color=(0.5, 0.5, 0.5))\n    ax.plot(frequencies, power_spectra[:, pop_idx],\n            color='black', linestyle='dashed', zorder=2)\n    ax.plot(frequencies, low_gamma_subcircuit_power_spectra[:, pop_idx],\n            color='black', zorder=10)\n    ax.set_yscale('log')\n\n    ax.set_title(labels[pop_idx])\n    ax.set_xlabel(r'frequency (1/$s$)')\n\n    ax.set_xticks([20, 40, 60, 80])\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    ax.set_yticks([1e-5,1e-3,1e-1])\n    ax.set_ylabel(r'$|C(\\omega)|$')\n    ax.set_xlim([0.0, 100.0])\n    ax.set_ylim([1e-6, 1e0])\n\n\n# take out 23E -> 4I \nfor gs, pop_idx in zip(gsB, [0,2,3]):\n    ax = fig.add_subplot(gs)\n\n    ax.plot(result['oscillation_origin']['B']['freq_sim'],\n                result['oscillation_origin']['B'][f'power{pop_idx}_sim'], \n                color=(0.8, 0.8, 0.8))\n    ax.plot(result['oscillation_origin']['B']['freq_sim_av'],\n                result['oscillation_origin']['B'][f'power{pop_idx}_sim_av'], \n                color=(0.5, 0.5, 0.5))\n    ax.plot(frequencies, power_spectra[:, pop_idx],\n            color='black', linestyle='dashed', zorder=2)\n    ax.plot(frequencies, without_23E_4I_power_spectra[:, pop_idx],\n            color='black', zorder=10)\n    ax.set_yscale('log')\n\n    ax.set_title(labels[pop_idx])\n    ax.set_xlabel(r'frequency (1/$s$)')\n\n    ax.set_xticks([20, 40, 60, 80])\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    ax.set_yticks([1e-5,1e-3,1e-1])\n    ax.set_ylabel(r'$|C(\\omega)|$')\n    ax.set_xlim([0.0, 100.0])\n    ax.set_ylim([1e-6, 1e0])\n\nplt.savefig('figures/power_spectra_of_subcircuits_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.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}