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