{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# PyNEST Microcircuit: Helper Functions\n\nHelper functions for network construction, simulation and evaluation of the\nmicrocircuit.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "from matplotlib.patches import Polygon\nimport matplotlib.pyplot as plt\nimport os\nimport sys\nimport numpy as np\nif 'DISPLAY' not in os.environ:\n    import matplotlib\n    matplotlib.use('Agg')\n\n\ndef num_synapses_from_conn_probs(conn_probs, popsize1, popsize2):\n    \"\"\"Computes the total number of synapses between two populations from\n    connection probabilities.\n\n    Here it is irrelevant which population is source and which target.\n\n    Parameters\n    ----------\n    conn_probs\n        Matrix of connection probabilities.\n    popsize1\n        Size of first poulation.\n    popsize2\n        Size of second population.\n\n    Returns\n    -------\n    num_synapses\n        Matrix of synapse numbers.\n    \"\"\"\n    prod = np.outer(popsize1, popsize2)\n    num_synapses = np.log(1. - conn_probs) / np.log((prod - 1.) / prod)\n    return num_synapses\n\n\ndef postsynaptic_potential_to_current(C_m, tau_m, tau_syn):\n    \"\"\" Computes a factor to convert postsynaptic potentials to currents.\n\n    The time course of the postsynaptic potential ``v`` is computed as\n    :math: `v(t)=(i*h)(t)`\n    with the exponential postsynaptic current\n    :math:`i(t)=J\\mathrm{e}^{-t/\\tau_\\mathrm{syn}}\\Theta (t)`,\n    the voltage impulse response\n    :math:`h(t)=\\frac{1}{\\tau_\\mathrm{m}}\\mathrm{e}^{-t/\\tau_\\mathrm{m}}\\Theta (t)`,\n    and\n    :math:`\\Theta(t)=1` if :math:`t\\geq 0` and zero otherwise.\n\n    The ``PSP`` is considered as the maximum of ``v``, i.e., it is\n    computed by setting the derivative of ``v(t)`` to zero.\n    The expression for the time point at which ``v`` reaches its maximum\n    can be found in Eq. 5 of [1]_.\n\n    The amplitude of the postsynaptic current ``J`` corresponds to the\n    synaptic weight ``PSC``.\n\n    References\n    ----------\n    .. [1] Hanuschkin A, Kunkel S, Helias M, Morrison A and Diesmann M (2010)\n           A general and efficient method for incorporating precise spike times\n           in globally time-driven simulations.\n           Front. Neuroinform. 4:113.\n           DOI: `10.3389/fninf.2010.00113 <https://doi.org/10.3389/fninf.2010.00113>`__.\n\n    Parameters\n    ----------\n    C_m\n        Membrane capacitance (in pF).\n    tau_m\n        Membrane time constant (in ms).\n    tau_syn\n        Synaptic time constant (in ms).\n\n    Returns\n    -------\n    PSC_over_PSP\n        Conversion factor to be multiplied to a `PSP` (in mV) to obtain a `PSC`\n        (in pA).\n\n    \"\"\"\n    sub = 1. / (tau_syn - tau_m)\n    pre = tau_m * tau_syn / C_m * sub\n    frac = (tau_m / tau_syn) ** sub\n\n    PSC_over_PSP = 1. / (pre * (frac**tau_m - frac**tau_syn))\n    return PSC_over_PSP\n\n\ndef dc_input_compensating_poisson(bg_rate, K_ext, tau_syn, PSC_ext):\n    \"\"\" Computes DC input if no Poisson input is provided to the microcircuit.\n\n    Parameters\n    ----------\n    bg_rate\n        Rate of external Poisson generators (in spikes/s).\n    K_ext\n        External indegrees.\n    tau_syn\n        Synaptic time constant (in ms).\n    PSC_ext\n        Weight of external connections (in pA).\n\n    Returns\n    -------\n    DC\n        DC input (in pA) which compensates lacking Poisson input.\n    \"\"\"\n    DC = bg_rate * K_ext * PSC_ext * tau_syn * 0.001\n    return DC\n\n\ndef adjust_weights_and_input_to_synapse_scaling(\n        full_num_neurons,\n        full_num_synapses,\n        K_scaling,\n        mean_PSC_matrix,\n        PSC_ext,\n        tau_syn,\n        full_mean_rates,\n        DC_amp,\n        poisson_input,\n        bg_rate,\n        K_ext):\n    \"\"\" Adjusts weights and external input to scaling of indegrees.\n\n    The recurrent and external weights are adjusted to the scaling\n    of the indegrees. Extra DC input is added to compensate for the\n    scaling in order to preserve the mean and variance of the input.\n\n    Parameters\n    ----------\n    full_num_neurons\n        Total numbers of neurons.\n    full_num_synapses\n        Total numbers of synapses.\n    K_scaling\n        Scaling factor for indegrees.\n    mean_PSC_matrix\n        Weight matrix (in pA).\n    PSC_ext\n        External weight (in pA).\n    tau_syn\n        Synaptic time constant (in ms).\n    full_mean_rates\n        Firing rates of the full network (in spikes/s).\n    DC_amp\n        DC input current (in pA).\n    poisson_input\n        True if Poisson input is used.\n    bg_rate\n        Firing rate of Poisson generators (in spikes/s).\n    K_ext\n        External indegrees.\n\n    Returns\n    -------\n    PSC_matrix_new\n        Adjusted weight matrix (in pA).\n    PSC_ext_new\n        Adjusted external weight (in pA).\n    DC_amp_new\n        Adjusted DC input (in pA).\n\n    \"\"\"\n    PSC_matrix_new = mean_PSC_matrix / np.sqrt(K_scaling)\n    PSC_ext_new = PSC_ext / np.sqrt(K_scaling)\n\n    # recurrent input of full network\n    indegree_matrix = \\\n        full_num_synapses / full_num_neurons[:, np.newaxis]\n    input_rec = np.sum(mean_PSC_matrix * indegree_matrix * full_mean_rates,\n                       axis=1)\n\n    DC_amp_new = DC_amp \\\n        + 0.001 * tau_syn * (1. - np.sqrt(K_scaling)) * input_rec\n\n    if poisson_input:\n        input_ext = PSC_ext * K_ext * bg_rate\n        DC_amp_new += 0.001 * tau_syn * (1. - np.sqrt(K_scaling)) * input_ext\n    return PSC_matrix_new, PSC_ext_new, DC_amp_new\n\n\ndef plot_raster(path, name, begin, end, N_scaling):\n    \"\"\" Creates a spike raster plot of the network activity.\n\n    Parameters\n    -----------\n    path\n        Path where the spike times are stored.\n    name\n        Name of the spike recorder.\n    begin\n        Time point (in ms) to start plotting spikes (included).\n    end\n        Time point (in ms) to stop plotting spikes (included).\n    N_scaling\n        Scaling factor for number of neurons.\n\n    Returns\n    -------\n    None\n\n    \"\"\"\n    fs = 18  # fontsize\n    ylabels = ['L2/3', 'L4', 'L5', 'L6']\n    color_list = np.tile(['#595289', '#af143c'], 4)\n\n    sd_names, node_ids, data = __load_spike_times(path, name, begin, end)\n    last_node_id = node_ids[-1, -1]\n    mod_node_ids = np.abs(node_ids - last_node_id) + 1\n\n    label_pos = [(mod_node_ids[i, 0] + mod_node_ids[i + 1, 1]) /\n                 2. for i in np.arange(0, 8, 2)]\n\n    stp = 1\n    if N_scaling > 0.1:\n        stp = int(10. * N_scaling)\n        print('  Only spikes of neurons in steps of {} are shown.'.format(stp))\n\n    plt.figure(figsize=(8, 6))\n    for i, n in enumerate(sd_names):\n        times = data[i]['time_ms']\n        neurons = np.abs(data[i]['sender'] - last_node_id) + 1\n        plt.plot(times[::stp], neurons[::stp], '.', color=color_list[i])\n    plt.xlabel('time [ms]', fontsize=fs)\n    plt.xticks(fontsize=fs)\n    plt.yticks(label_pos, ylabels, fontsize=fs)\n    plt.savefig(os.path.join(path, 'raster_plot.png'), dpi=300)\n\n\ndef firing_rates(path, name, begin, end):\n    \"\"\" Computes mean and standard deviation of firing rates per population.\n\n    The firing rate of each neuron in each population is computed and stored\n    in a .dat file in the directory of the spike recorders. The mean firing\n    rate and its standard deviation are printed out for each population.\n\n    Parameters\n    -----------\n    path\n        Path where the spike times are stored.\n    name\n        Name of the spike recorder.\n    begin\n        Time point (in ms) to start calculating the firing rates (included).\n    end\n        Time point (in ms) to stop calculating the firing rates (included).\n\n    Returns\n    -------\n    None\n\n    \"\"\"\n    sd_names, node_ids, data = __load_spike_times(path, name, begin, end)\n    all_mean_rates = []\n    all_std_rates = []\n    for i, n in enumerate(sd_names):\n        senders = data[i]['sender']\n        # 1 more bin than node ids per population\n        bins = np.arange(node_ids[i, 0], node_ids[i, 1] + 2)\n        spike_count_per_neuron, _ = np.histogram(senders, bins=bins)\n        rate_per_neuron = spike_count_per_neuron * 1000. / (end - begin)\n        np.savetxt(os.path.join(path, ('rate' + str(i) + '.dat')),\n                   rate_per_neuron)\n        # zeros are included\n        all_mean_rates.append(np.mean(rate_per_neuron))\n        all_std_rates.append(np.std(rate_per_neuron))\n    print('Mean rates: {} spikes/s'.format(np.around(all_mean_rates, decimals=3)))\n    print('Standard deviation of rates: {} spikes/s'.format(\n        np.around(all_std_rates, decimals=3)))\n\n\ndef boxplot(path, populations):\n    \"\"\" Creates a boxblot of the firing rates of all populations.\n\n    To create the boxplot, the firing rates of each neuron in each population\n    need to be computed with the function ``firing_rate()``.\n\n    Parameters\n    -----------\n    path\n        Path where the firing rates are stored.\n    populations\n        Names of neuronal populations.\n\n    Returns\n    -------\n    None\n\n    \"\"\"\n    fs = 18\n    pop_names = [string.replace('23', '2/3') for string in populations]\n    label_pos = list(range(len(populations), 0, -1))\n    color_list = ['#af143c', '#595289']\n    medianprops = dict(linestyle='-', linewidth=2.5, color='black')\n    meanprops = dict(linestyle='--', linewidth=2.5, color='lightgray')\n\n    rates_per_neuron_rev = []\n    for i in np.arange(len(populations))[::-1]:\n        rates_per_neuron_rev.append(\n            np.loadtxt(os.path.join(path, ('rate' + str(i) + '.dat'))))\n\n    plt.figure(figsize=(8, 6))\n    bp = plt.boxplot(rates_per_neuron_rev, 0, 'rs', 0, medianprops=medianprops,\n                     meanprops=meanprops, meanline=True, showmeans=True)\n    plt.setp(bp['boxes'], color='black')\n    plt.setp(bp['whiskers'], color='black')\n    plt.setp(bp['fliers'], color='red', marker='+')\n\n    # boxcolors\n    for i in np.arange(len(populations)):\n        boxX = []\n        boxY = []\n        box = bp['boxes'][i]\n        for j in list(range(5)):\n            boxX.append(box.get_xdata()[j])\n            boxY.append(box.get_ydata()[j])\n        boxCoords = list(zip(boxX, boxY))\n        k = i % 2\n        boxPolygon = Polygon(boxCoords, facecolor=color_list[k])\n        plt.gca().add_patch(boxPolygon)\n    plt.xlabel('firing rate [spikes/s]', fontsize=fs)\n    plt.yticks(label_pos, pop_names, fontsize=fs)\n    plt.xticks(fontsize=fs)\n    plt.savefig(os.path.join(path, 'box_plot.png'), dpi=300)\n\n\ndef __gather_metadata(path, name):\n    \"\"\" Reads names and ids of spike recorders and first and last ids of\n    neurons in each population.\n\n    If the simulation was run on several threads or MPI-processes, one name per\n    spike recorder per MPI-process/thread is extracted.\n\n    Parameters\n    ------------\n    path\n        Path where the spike recorder files are stored.\n    name\n        Name of the spike recorder, typically ``spike_recorder``.\n\n    Returns\n    -------\n    sd_files\n        Names of all files written by spike recorders.\n    sd_names\n        Names of all spike recorders.\n    node_ids\n        Lowest and highest id of nodes in each population.\n\n    \"\"\"\n    # load filenames\n    sd_files = []\n    sd_names = []\n    for fn in sorted(os.listdir(path)):\n        if fn.startswith(name):\n            sd_files.append(fn)\n            # spike recorder name and its ID\n            fnsplit = '-'.join(fn.split('-')[:-1])\n            if fnsplit not in sd_names:\n                sd_names.append(fnsplit)\n\n    # load node IDs\n    node_idfile = open(path + 'population_nodeids.dat', 'r')\n    node_ids = []\n    for l in node_idfile:\n        node_ids.append(l.split())\n    node_ids = np.array(node_ids, dtype='i4')\n    return sd_files, sd_names, node_ids\n\n\ndef __load_spike_times(path, name, begin, end):\n    \"\"\" Loads spike times of each spike recorder.\n\n    Parameters\n    ----------\n    path\n        Path where the files with the spike times are stored.\n    name\n        Name of the spike recorder.\n    begin\n        Time point (in ms) to start loading spike times (included).\n    end\n        Time point (in ms) to stop loading spike times (included).\n\n    Returns\n    -------\n    data\n        Dictionary containing spike times in the interval from ``begin``\n        to ``end``.\n\n    \"\"\"\n    sd_files, sd_names, node_ids = __gather_metadata(path, name)\n    data = {}\n    dtype = {'names': ('sender', 'time_ms'),  # as in header\n             'formats': ('i4', 'f8')}\n    for i, name in enumerate(sd_names):\n        data_i_raw = np.array([[]], dtype=dtype)\n        for j, f in enumerate(sd_files):\n            if name in f:\n                # skip header while loading\n                ld = np.loadtxt(os.path.join(path, f), skiprows=3, dtype=dtype)\n                data_i_raw = np.append(data_i_raw, ld)\n\n        data_i_raw = np.sort(data_i_raw, order='time_ms')\n        # begin and end are included if they exist\n        low = np.searchsorted(data_i_raw['time_ms'], v=begin, side='left')\n        high = np.searchsorted(data_i_raw['time_ms'], v=end, side='right')\n        data[i] = data_i_raw[low:high]\n    return sd_names, node_ids, data"
      ]
    }
  ],
  "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}