Post-Processing and Pipelining (from the Tutorial)¶
Here you find an example of post-processing.
It consists of a main script main.py for the three phases pre-processing, run phase and post-processing of a single neuron simulation and a analysis.py file giving an example of a potential data analysis encompassing plotting the results. Moreover, there exists a pipeline.py file to crunch all first three phases into a single function.
A detail explanation of the example can be found in the Tutorial section.
Download: main.py
Download: analysis.py
Download: pipeline.py
Main¶
__author__ = 'robert'
import numpy as np
import pandas as pd
import logging
import os # For path names working under Linux and Windows
from pypet import Environment, cartesian_product
def run_neuron(traj):
"""Runs a simulation of a model neuron.
:param traj:
Container with all parameters.
:return:
An estimate of the firing rate of the neuron
"""
# Extract all parameters from `traj`
V_init = traj.par.neuron.V_init
I = traj.par.neuron.I
tau_V = traj.par.neuron.tau_V
tau_ref = traj.par.neuron.tau_ref
dt = traj.par.simulation.dt
duration = traj.par.simulation.duration
steps = int(duration / float(dt))
# Create some containers for the Euler integration
V_array = np.zeros(steps)
V_array[0] = V_init
spiketimes = [] # List to collect all times of action potentials
# Do the Euler integration:
print('Starting Euler Integration')
for step in range(1, steps):
if V_array[step-1] >= 1:
# The membrane potential crossed the threshold and we mark this as
# an action potential
V_array[step] = 0
spiketimes.append((step-1)*dt)
elif spiketimes and step * dt - spiketimes[-1] <= tau_ref:
# We are in the refractory period, so we simply clamp the voltage
# to 0
V_array[step] = 0
else:
# Euler Integration step:
dV = -1/tau_V * V_array[step-1] + I
V_array[step] = V_array[step-1] + dV*dt
print('Finished Euler Integration')
# Add the voltage trace and spike times
traj.f_add_result('neuron.$', V=V_array, nspikes=len(spiketimes),
comment='Contains the development of the membrane potential over time '
'as well as the number of spikes.')
# This result will be renamed to `traj.results.neuron.run_XXXXXXXX`.
# And finally we return the estimate of the firing rate
return len(spiketimes) / float(traj.par.simulation.duration) *1000
# *1000 since we have defined duration in terms of milliseconds
def neuron_postproc(traj, result_list):
"""Postprocessing, sorts computed firing rates into a table
:param traj:
Container for results and parameters
:param result_list:
List of tuples, where first entry is the run index and second is the actual
result of the corresponding run.
:return:
"""
# Let's create a pandas DataFrame to sort the computed firing rate according to the
# parameters. We could have also used a 2D numpy array.
# But a pandas DataFrame has the advantage that we can index into directly with
# the parameter values without translating these into integer indices.
I_range = traj.par.neuron.f_get('I').f_get_range()
ref_range = traj.par.neuron.f_get('tau_ref').f_get_range()
I_index = sorted(set(I_range))
ref_index = sorted(set(ref_range))
rates_frame = pd.DataFrame(columns=ref_index, index=I_index)
# This frame is basically a two dimensional table that we can index with our
# parameters
# Now iterate over the results. The result list is a list of tuples, with the
# run index at first position and our result at the second
for result_tuple in result_list:
run_idx = result_tuple[0]
firing_rates = result_tuple[1]
I_val = I_range[run_idx]
ref_val = ref_range[run_idx]
rates_frame.loc[I_val, ref_val] = firing_rates # Put the firing rate into the
# data frame
# Finally we going to store our new firing rate table into the trajectory
traj.f_add_result('summary.firing_rates', rates_frame=rates_frame,
comment='Contains a pandas data frame with all firing rates.')
def add_parameters(traj):
"""Adds all parameters to `traj`"""
print('Adding Parameters')
traj.f_add_parameter('neuron.V_init', 0.0,
comment='The initial condition for the '
'membrane potential')
traj.f_add_parameter('neuron.I', 0.0,
comment='The externally applied current.')
traj.f_add_parameter('neuron.tau_V', 10.0,
comment='The membrane time constant in milliseconds')
traj.f_add_parameter('neuron.tau_ref', 5.0,
comment='The refractory period in milliseconds '
'where the membrane potnetial '
'is clamped.')
traj.f_add_parameter('simulation.duration', 1000.0,
comment='The duration of the experiment in '
'milliseconds.')
traj.f_add_parameter('simulation.dt', 0.1,
comment='The step size of an Euler integration step.')
def add_exploration(traj):
"""Explores different values of `I` and `tau_ref`."""
print('Adding exploration of I and tau_ref')
explore_dict = {'neuron.I': np.arange(0, 1.01, 0.01).tolist(),
'neuron.tau_ref': [5.0, 7.5, 10.0]}
explore_dict = cartesian_product(explore_dict, ('neuron.tau_ref', 'neuron.I'))
# The second argument, the tuple, specifies the order of the cartesian product,
# The variable on the right most side changes fastest and defines the
# 'inner for-loop' of the cartesian product
traj.f_explore(explore_dict)
def main():
filename = os.path.join('hdf5', 'FiringRate.hdf5')
env = Environment(trajectory='FiringRate',
comment='Experiment to measure the firing rate '
'of a leaky integrate and fire neuron. '
'Exploring different input currents, '
'as well as refractory periods',
add_time=False, # We don't want to add the current time to the name,
log_stdout=True,
log_config='DEFAULT',
multiproc=True,
ncores=2, #My laptop has 2 cores ;-)
wrap_mode='QUEUE',
filename=filename,
overwrite_file=True)
traj = env.trajectory
# Add parameters
add_parameters(traj)
# Let's explore
add_exploration(traj)
# Ad the postprocessing function
env.add_postprocessing(neuron_postproc)
# Run the experiment
env.run(run_neuron)
# Finally disable logging and close all log-files
env.disable_logging()
if __name__ == '__main__':
main()
Analysis¶
__author__ = 'robert'
import os
from pypet import Trajectory
import matplotlib.pyplot as plt
def main():
# This time we don't need an environment since we just going to look
# at data in the trajectory
traj = Trajectory('FiringRate', add_time=False)
# Let's load the trajectory from the file
# Only load the parameters, we will load the results on the fly as we need them
filename = os.path.join('hdf5', 'FiringRate.hdf5')
traj.f_load(load_parameters=2, load_derived_parameters=0, load_results=0,
load_other_data=0, filename=filename)
# We'll simply use auto loading so all data will be loaded when needed.
traj.v_auto_load = True
rates_frame = traj.res.summary.firing_rates.rates_frame
# Here we load the data automatically on the fly
plt.figure()
plt.subplot(2,1,1)
#Let's iterate through the columns and plot the different firing rates :
for tau_ref, I_col in rates_frame.iteritems():
plt.plot(I_col.index, I_col, label='Avg. Rate for tau_ref=%s' % str(tau_ref))
# Label the plot
plt.xlabel('I')
plt.ylabel('f[Hz]')
plt.title('Firing as a function of input current `I`')
plt.legend(loc='best')
# Also let's plot an example run, how about run 13 ?
example_run = 13
traj.v_idx = example_run # We make the trajectory behave as a single run container.
# This short statement has two major effects:
# a) all explored parameters are set to the value of run 13,
# b) if there are tree nodes with names other than the current run aka `run_00000013`
# they are simply ignored, if we use the `$` sign or the `crun` statement,
# these are translated into `run_00000013`.
# Get the example data
example_I = traj.I
example_tau_ref = traj.tau_ref
example_V = traj.results.neuron.crun.V # Here crun stands for run_00000013
# We need the time step...
dt = traj.dt
# ...to create an x-axis for the plot
dt_array = [irun * dt for irun in range(len(example_V))]
# And plot the development of V over time,
# Since this is rather repetitive, we only
# plot the first eighth of it.
plt.subplot(2,1,2)
plt.plot(dt_array, example_V)
plt.xlim((0, dt*len(example_V)/8))
# Label the axis
plt.xlabel('t[ms]')
plt.ylabel('V')
plt.title('Example of development of V for I=%s, tau_ref=%s in run %d' %
(str(example_I), str(example_tau_ref), traj.v_idx))
# And let's take a look at it
plt.show()
# Finally revoke the `traj.v_idx=13` statement and set everything back to normal.
# Since our analysis is done here, we could skip that, but it is always a good idea
# to do that.
traj.f_restore_default()
if __name__ == '__main__':
main()
Pipelining¶
Additionally, you can use pipelining.
Since these three steps pre-processing, run-phase, post-processing define a common pipeline, you can actually also make pypet supervise all three steps at once.
You can define a pipeline function, that does the pre-processing and returns the job function plus some optional arguments and the post-processing function with some other optional arguments.
So, you could define the following pipeline function.
The pipeline function has to only accept the trajectory as first argument and
has to return 2 tuples, one for the run function and one for the
post-processing. Since none of our functions takes any other arguments than the trajectory
(and the pos-processing function the result list) we simply return an empty
tuple () for no arguments and an empty dictionary {} for no keyword arguments.
And that’s it, than everything including the pre-processing and addition of parameters is supervised by pypet. Check out the source code below:
__author__ = 'robert'
import logging
import os # For path names working under Windows and Linux
from main import add_parameters, add_exploration, run_neuron, neuron_postproc
from pypet import Environment
def mypipeline(traj):
"""A pipeline function that defines the entire experiment
:param traj:
Container for results and parameters
:return:
Two tuples. First tuple contains the actual run function plus additional
arguments (yet we have none). Second tuple contains the
postprocessing function including additional arguments.
"""
add_parameters(traj)
add_exploration(traj)
return (run_neuron,(),{}), (neuron_postproc,(),{})
def main():
filename = os.path.join('hdf5', 'FiringRate.hdf5')
env = Environment(trajectory='FiringRatePipeline',
comment='Experiment to measure the firing rate '
'of a leaky integrate and fire neuron. '
'Exploring different input currents, '
'as well as refractory periods',
add_time=False, # We don't want to add the current time to the name,
log_stdout=True,
multiproc=True,
ncores=2, #My laptop has 2 cores ;-)
filename=filename,
overwrite_file=True)
env.pipeline(mypipeline)
# Finally disable logging and close all log-files
env.disable_logging()
if __name__ == '__main__':
main()