{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Eigenvalue Trajectories (Bos 2016)\n\nHere we calculate the eigenvalue trajectories of the :cite:t:`potjans2014` \nmicrocircuit model including modifications made in :cite:t:`bos2016`.\n\nThis example reproduces Fig 4. 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\nfrom collections import defaultdict\n\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": [
        "Calculate all necessary quantities and finally the eigenvalues of the \neffective connectivity matrix. To add further details to the plot, we\nhere also compute the sensitivity measure for all eigenmodes\nand 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)\neigenvalues = np.linalg.eig(\n    nnmt.lif.exp.effective_connectivity(microcircuit))[0]\nresorted_eigenvalues, new_indices = (\n    nnmt.lif.exp._match_eigenvalues_across_frequencies(eigenvalues))\nsensitivity_dict = nnmt.lif.exp.sensitivity_measure_all_eigenmodes(\n    network=microcircuit)\n# calculate the power spectra\npower_spectra = nnmt.lif.exp.power_spectra(microcircuit).T"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Identify the indices of the eigenvalues that should be plotted to \nreproduce Fig. 4 in Bos 2016.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for i in range(8):\n    i = str(i)\n    print(sensitivity_dict[i]['critical_frequency'])\n    \n# Here, we manually identified the following indices to correspond to the \n# markers shown in the publication.\neigenvalues_to_be_plotted = [str(i) for i in [0, 1, 2, 3, 5, 6]]\nprint(f'Eigenvalues to be plotted: {eigenvalues_to_be_plotted}')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "For comparison with the complete circuit, we here alter the indegree matrix \nsuch that it just consist of the individual isolated layers\nand calculate the corresponding eigenvalue spectra and power spectra.\nNote: This isolated layers are treated at the same working point as the whole\ncircuit.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "isolated_layers_results = defaultdict(str)\nfor i, layer in enumerate(['23', '4', '5', '6']):\n    print(f'Modify connectivity to obtain isolated {layer}.')\n    isolated_layers_results[layer] = defaultdict(str)\n\n    microcircuit_isolated_layers = 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')\n        \n    if reduce_frequency_resolution:\n        microcircuit_isolated_layers.change_parameters(\n            changed_analysis_params={'df': 5},\n            overwrite=True)\n        derived_analysis_params = (\n            microcircuit_isolated_layers._calculate_dependent_analysis_parameters())\n        microcircuit_isolated_layers.analysis_params.update(\n            derived_analysis_params)\n    \n    isolated_layers_results[layer]['network'] = microcircuit_isolated_layers\n    \n    # calculate working point for exponentially shape post synaptic currents\n    nnmt.lif.exp.working_point(\n        isolated_layers_results[layer]['network'], method='taylor')\n\n    reducing_matrix = np.zeros((8,8))\n    reducing_matrix[2*i:2*i+2, 2*i:2*i+2] = np.ones([2,2])\n    # for layer i, set indegree matrix such that it is isolated\n    isolated_layers_results[layer]['network'].network_params.update(\n        K = full_indegree_matrix*reducing_matrix)\n    \n    print('connectivity changed...')\n    # calculate the transfer function\n    nnmt.lif.exp.transfer_function(\n        isolated_layers_results[layer]['network'], method='taylor')\n    # calculate the delay distribution matrix\n    nnmt.network_properties.delay_dist_matrix(\n        isolated_layers_results[layer]['network'])\n    eigenvalue_spectra_layer = np.linalg.eig(\n            nnmt.lif.exp.effective_connectivity(\n                isolated_layers_results[layer]['network']))[0]\n    # calculate the power spectra\n    power_spectra_layer = nnmt.lif.exp.power_spectra(\n        isolated_layers_results[layer]['network']).T\n    \n    resorted_eigenvalue_spectra_layer, new_indices_layer = (\n        nnmt.lif.exp._match_eigenvalues_across_frequencies(\n        eigenvalue_spectra_layer))\n    \n    isolated_layers_results[layer]['eigenvalue_spectra'] = (\n        eigenvalue_spectra_layer)\n    isolated_layers_results[layer]['power_spectra'] = power_spectra_layer"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Calculate the sensitivity measure dictionary for each isolated layer.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "for i, layer in enumerate(['23', '4', '5', '6']):\n    layer = isolated_layers_results[layer]\n    layer['sensitivity_dict'] = (\n        nnmt.lif.exp.sensitivity_measure_all_eigenmodes(layer['network']))\n\n# identify which eigenvalues should be plotted to reproduce Fig.4 of Bos 2016\nfor i, layer in enumerate(['23', '4', '5', '6']):\n    for j in range(2):\n        j = str(j)\n        print(layer, j,\n              isolated_layers_results[layer][\n                  'sensitivity_dict'][j]['critical_frequency'])\n\nlayer_eigenvalue_tuples_to_be_plotted = [('23', '1'), ('4', '1'),\n                                         ('5', '1'), ('6', '1'),\n                                         ('23', '0'), ('5', '0')]"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally, we plot 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(1, 3, figure=fig)\n\n# Panel A\ngsA = gridspec.GridSpecFromSubplotSpec(2, 1,\n                                       subplot_spec=grid_specification[0])\n# Panel B\ngsB = gridspec.GridSpecFromSubplotSpec(2, 1,\n                                       subplot_spec=grid_specification[1])\n# Panel C\ngsC = gridspec.GridSpecFromSubplotSpec(2, 1,\n                                       subplot_spec=grid_specification[2])\n\n### Panel A\n# top\nax = fig.add_subplot(gsA[0])\nN = resorted_eigenvalues.shape[0] # N frequencies\ndc = 1/float(N)\nfor i in range(0, N, 3):\n    ax.plot(resorted_eigenvalues[i].real, \n            resorted_eigenvalues[i].imag, '.',\n                color=(0.9-0.9*i*dc, 0.9-0.9*i*dc, 0.9-0.9*i*dc),\n                markersize=1.0, zorder=1)\nax.scatter(1,0, s=15, color='r')\nax.set_ylim([-4, 6.5])\nax.set_xlim([-11.5, 2])\nax.set_xticks([-9, -6, -3, 0])\nax.set_yticks([-3, 0, 3, 6])\nax.set_ylabel('Im($\\lambda(\\omega)$)')\n\n# bottom\nax = fig.add_subplot(gsA[1])\nfor i in range(0, N, 3):\n    ax.plot(resorted_eigenvalues[i].real, \n            resorted_eigenvalues[i].imag, '.',\n                color=(0.9-0.9*i*dc, 0.9-0.9*i*dc, 0.9-0.9*i*dc),\n                markersize=1.0, zorder=1)\n# frequencies where eigenvalue trajectory is closest to one\nfmaxs = [np.round(sensitivity_dict[i]['critical_frequency'],1) \n         for i in eigenvalues_to_be_plotted]\nmarkers = ['<', '>', '^', 'v', 'o', '+']\nfor i, eig_index in enumerate(eigenvalues_to_be_plotted):\n    eigc = sensitivity_dict[eig_index]['critical_eigenvalue']\n    ax.plot(eigc.real, eigc.imag, markers[i], color='black',#colors_array[i],\n                mew=1, ms=4, label=str(fmaxs[i])+'Hz')\nax.legend(bbox_to_anchor=(-0.35, -0.9, 1.6, 0.5), loc='center', \n                ncol=3, mode=\"expand\", borderaxespad=3.5, fontsize=7,\n                numpoints = 1)\nax.scatter(1,0, s=5, color='r')\n# box = ax.get_position()\n# ax.set_position([box.x0, box.y0-box.height, box.width, box.height*2])\nax.set_xlabel('Re($\\lambda(\\omega)$)')\nax.set_ylabel('Im($\\lambda(\\omega)$)')\nax.set_ylim([-0.3, 0.5])\nax.set_xlim([0.1, 1.1])\nax.set_yticks([-0.2, 0, 0.2, 0.4])\nax.set_xticks([0.2, 0.4, 0.6, 0.8, 1.0])\n\n\n### Panel B\n\ndef get_parameter_plot():\n    colors = [[] for i in range(4)]\n    colors[1] = [(0.0, 0.7, 0.0), (0.0, 1.0, 0.0)]\n    colors[3] = [(0.0, 0.0, 0.4), (0.0, 0.0, 1.0)]\n    colors[0] = [(0.7, 0.0, 0.0), (1.0, 0.0, 0.0)]\n    colors[2] = [(0.5, 0.0, 0.5), (1.0, 0.0, 1.0)]\n    return colors\n\ncolors = get_parameter_plot()\n\ndef get_color(i, layer):\n    cont_colors = [(1.0-i*dc, 0.0, 0.0), (0.0, 1.0-i*dc, 0.0), \n                    (1.0-i*dc, 0.0, 1.0-i*dc), (0.0, 0.0, 1.0-i*dc)]\n    index = ['23', '4', '5',  '6'].index(layer)\n    return cont_colors[index]\n\n# top\nax = fig.add_subplot(gsB[0])\nN = resorted_eigenvalues.shape[1] # N populations\ndc = 1/float(N)\nfor i, layer in enumerate(['23', '4', '5', '6']):\n    # [:, 2*layer:2*layer+2] in the original code serves to \n    # plot only the non-zero eigenspectra \n    for j in range(N):\n        plt.scatter(\n            isolated_layers_results[layer]['eigenvalue_spectra'][:, j].real,\n            isolated_layers_results[layer]['eigenvalue_spectra'][:, j].imag,\n            color=get_color(j, layer),\n            s=0.5, zorder=1)\nax.scatter(1,0, s=15, color='r')\nax.set_ylim([-4, 6.5])\nax.set_xlim([-11.5, 2])\nax.set_xticks([-9, -6, -3, 0])\nax.set_yticks([-3, 0, 3, 6])\nax.set_ylabel('Im($\\lambda(\\omega)$)')\n\n# bottom\nax = fig.add_subplot(gsB[1])\nfor i, layer in enumerate(['23', '4', '5', '6']):\n    for j in range(N):\n        plt.scatter(\n            isolated_layers_results[layer]['eigenvalue_spectra'][:, j].real,\n            isolated_layers_results[layer]['eigenvalue_spectra'][:, j].imag,\n            color=get_color(j, layer),\n            s=0.5, zorder=1)\n        \n# frequencies where eigenvalue trajectory is closest to one\nfmaxs = [np.round(\n    isolated_layers_results[i[0]][\n        'sensitivity_dict'][i[1]]['critical_frequency'],1) for i \n         in layer_eigenvalue_tuples_to_be_plotted]\n\nmarkers = ['<', '>', '^', 'v', 'o', '+']\nfor i, (layer, eig_index) in enumerate(layer_eigenvalue_tuples_to_be_plotted):\n    eigc = isolated_layers_results[layer][\n        'sensitivity_dict'][eig_index]['critical_eigenvalue']\n    ax.plot(eigc.real, eigc.imag, markers[i], color='black',#colors_array[i],\n                mew=1, ms=4, label=str(fmaxs[i])+'Hz')\nax.legend(bbox_to_anchor=(-0.35, -0.9, 1.6, 0.5), loc='center', \n                ncol=3, mode=\"expand\", borderaxespad=3.5, fontsize=7,\n                numpoints = 1)\nax.scatter(1,0, s=5, color='r')\nax.set_xlabel('Re($\\lambda(\\omega)$)')\nax.set_ylabel('Im($\\lambda(\\omega)$)')\nax.set_ylim([-0.3, 0.5])\nax.set_xlim([0.1, 1.1])\nax.set_yticks([-0.2, 0, 0.2, 0.4])\nax.set_xticks([0.2, 0.4, 0.6, 0.8, 1.0])\n\n### panel C ##\nfreqs_isolated_layers = microcircuit_isolated_layers.analysis_params[\n    'omegas']/(2*np.pi)\nfrequencies = microcircuit.analysis_params['omegas']/(2*np.pi)\n\n# loop across layers 23 and 4\nfor i, layer in enumerate(['23', '4']):\n    ax = fig.add_subplot(gsC[i])\n\n    # loop across excitatory and inhibitory\n    for j in [0,1]:\n        ax.plot(freqs_isolated_layers,\n                isolated_layers_results[layer]['power_spectra'][j+2*i], \n                color=colors[i][j])\n        ax.plot(frequencies, power_spectra[2*i+j], \n                color='black', linestyle='dashed')\n        \n    ax.set_yscale('log')\n    ax.set_xticks([100, 200, 300])\n    ax.set_ylabel('power')\n    ax.set_yticks([1e-2, 1e-4])\n\n    if i == 0:\n        ax.set_ylim([5*1e-6, 5*1e-2])\n    else:\n        ax.set_ylim([2*1e-5, 4*1e-1])\n        ax.set_xlabel('frequency $f$(1/$s$)')\n        \nplt.savefig('figures/eigenvalue_trajectories_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.13"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}