"""
Functions to produce interactive plots
"""
from __future__ import unicode_literals
from __future__ import print_function
from __future__ import division
from __future__ import absolute_import
from future import standard_library
standard_library.install_aliases()
from netpyne import __gui__
if __gui__:
import matplotlib.pyplot as plt
from matplotlib import mlab
from matplotlib_scalebar import scalebar
from numbers import Number
from .utils import colorList, exception, getSpktSpkid, _showFigure, _saveFigData, getCellsInclude, syncMeasure, _smooth1d
import numpy as np
import pandas as pd
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive raster
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotRaster(include=['allCells'], timeRange=None, maxSpikes=1e8, orderBy='gid', orderInverse=False, popRates=False, spikeHist=False, spikeHistBin=5, syncLines=False, marker='circle', markerSize=3, popColors=None, saveData=None, saveFig=None, showFig=False):
"""Creates an interactive html raster plot of network cells using Bokeh.
Parameters
----------
include : list
Cells to include in the plot.
**Default:**
``['allCells']`` plots all cells
**Options:**
``['all']`` plots all cells and stimulations,
``['allNetStims']`` plots just stimulations,
``['popName1']`` plots a single population,
``['popName1', 'popName2']`` plots multiple populations,
``[120]`` plots a single cell,
``[120, 130]`` plots multiple cells,
``[('popName1', 56)]`` plots a cell from a specific population,
``[('popName1', [0, 1]), ('popName2', [4, 5, 6])]``, plots cells from multiple populations
timeRange : list [start, stop]
Time range to plot.
**Default:**
``None`` plots entire time range
maxSpikes : int
Maximum number of spikes to be plotted.
**Default:** ``1e8``
orderBy : str
Unique numeric cell property by which to order the y-axis.
**Default:** ``'gid'`` orders by cell ID
**Options:**
``'y'`` orders by cell y-location,
``'ynorm'`` orders by cell normalized y-location
orderInverse : bool
Inverts the y-axis order if ``True``.
**Default:** ``False``
popRates : bool
Include population firing rates on plot if ``True``.
**Default:** ``False``
spikeHist : bool
Include spike histogram (spikes/bin) on plot if ``True``.
**Default:** ``False``
spikeHistBin : int
Size of bin in ms to use for spike histogram.
**Default:** ``5``
syncLines : bool
Calculate synchrony measure and plot vertical lines for each spike to evidence synchrony if ``True``.
**Default:** ``False``
marker : str
`Bokeh marker <https://docs.bokeh.org/en/latest/docs/gallery/markers.html>`_ for each spike.
**Default:** ``'circle'``
markerSize : int
Size of Bokeh marker for each spike.
**Default:** ``3``
popColors : dict
Dictionary with custom color (value) used for each population (key).
**Default:** ``None`` uses standard colors
saveData : bool or str
Whether and where to save the data used to generate the plot.
**Default:** ``False``
**Options:** ``True`` autosaves the data,
``'/path/filename.ext'`` saves to a custom path and filename, valid file extensions are ``'.pkl'`` and ``'.json'``
saveFig : bool or str
Whether and where to save the figure.
**Default:** ``False``
**Options:** ``True`` autosaves the figure,
``'/path/filename.html'`` saves to a custom path and filename, only valid file extension is ``'.html'``
showFig : bool
Shows the figure if ``True``.
**Default:** ``True``
Returns
-------
(str, dict)
A tuple consisting of a string containing the html for the figure and a dictionary containing the plot data.
See Also
--------
plotRaster :
iplotSpikeHist :
plotSpikeHist :
Examples
--------
>>> import netpyne, netpyne.examples.example
>>> out = netpyne.analysis.iplotRaster()
"""
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout
from bokeh.models import Legend
from bokeh.colors import RGB
from bokeh.models.annotations import Title
print('Plotting interactive raster ...')
TOOLS = 'hover,save,pan,box_zoom,reset,wheel_zoom',
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
popColorDict = None
if popColors is not None:
popColorDict = popColors.copy()
for pop, color in popColorDict.items():
popColorDict[pop] = RGB(*[round(f * 255) for f in color])
cells, cellGids, netStimLabels = getCellsInclude(include)
df = pd.DataFrame.from_records(cells)
df = pd.concat([df.drop('tags', axis=1), pd.DataFrame.from_records(df['tags'].tolist())], axis=1)
keep = ['pop', 'gid', 'conns']
if isinstance(orderBy, str) and orderBy not in cells[0]['tags']: # if orderBy property doesn't exist or is not numeric, use gid
orderBy = 'gid'
elif isinstance(orderBy, str) and not isinstance(cells[0]['tags'][orderBy], Number):
orderBy = 'gid'
if isinstance(orderBy, list):
keep = keep + list(set(orderBy) - set(keep))
elif orderBy not in keep:
keep.append(orderBy)
df = df[keep]
popLabels = [pop for pop in sim.net.allPops if pop in df['pop'].unique()] #preserves original ordering
if netStimLabels: popLabels.append('NetStims')
popColorsTmp = {popLabel: colors[ipop%len(colors)] for ipop,popLabel in enumerate(popLabels)} # dict with color for each pop
if popColorDict: popColorsTmp.update(popColorDict)
popColorDict = popColorsTmp
if len(cellGids) > 0:
gidColors = {cell['gid']: popColorDict[cell['tags']['pop']] for cell in cells} # dict with color for each gid
try:
sel, spkts, spkgids = getSpktSpkid(cellGids=cellGids, timeRange=timeRange, allCells=(include == ['allCells']))
except:
import sys
print((sys.exc_info()))
spkgids, spkts = [], []
sel = pd.DataFrame(columns=['spkt', 'spkid'])
sel['spkgidColor'] = sel['spkid'].map(gidColors)
sel['pop'] = sel['spkid'].map(df.set_index('gid')['pop'])
df['gidColor'] = df['pop'].map(popColorDict)
df.set_index('gid', inplace=True)
# Order by
if len(df) > 0:
ylabelText = 'Cells (ordered by %s)'%(orderBy)
df = df.sort_values(by=orderBy)
sel['spkind'] = sel['spkid'].apply(df.index.get_loc)
else:
sel = pd.DataFrame(columns=['spkt', 'spkid', 'spkind'])
ylabelText = ''
# Add NetStim spikes
numCellSpks = len(sel)
numNetStims = 0
for netStimLabel in netStimLabels:
netStimSpks = [spk for cell,stims in sim.allSimData['stims'].items() \
for stimLabel,stimSpks in stims.items() for spk in stimSpks if stimLabel == netStimLabel]
if len(netStimSpks) > 0:
# lastInd = max(spkinds) if len(spkinds)>0 else 0
lastInd = sel['spkind'].max() if len(sel['spkind']) > 0 else 0
spktsNew = netStimSpks
spkindsNew = [lastInd+1+i for i in range(len(netStimSpks))]
ns = pd.DataFrame(list(zip(spktsNew, spkindsNew)), columns=['spkt', 'spkind'])
ns['spkgidColor'] = popColorDict['netStims']
sel = pd.concat([sel, ns])
numNetStims += 1
else:
pass
#print netStimLabel+' produced no spikes'
if len(cellGids)>0 and numNetStims:
ylabelText = ylabelText + ' and NetStims (at the end)'
elif numNetStims:
ylabelText = ylabelText + 'NetStims'
if numCellSpks+numNetStims == 0:
print('No spikes available to plot raster')
return None
# Time Range
if timeRange == [0, sim.cfg.duration]:
pass
elif timeRange is None:
timeRange = [0, sim.cfg.duration]
else:
sel = sel.query('spkt >= @timeRange[0] and spkt <= @timeRange[1]')
# Limit to max spikes
if (len(sel)>maxSpikes):
print((' Showing only the first %i out of %i spikes' % (maxSpikes, len(sel)))) # Limit num of spikes
if numNetStims: # sort first if have netStims
sel = sel.sort_values(by='spkt')
sel = sel.iloc[:maxSpikes]
timeRange[1] = sel['spkt'].max()
# Plot stats
gidPops = df['pop'].tolist()
popNumCells = [float(gidPops.count(pop)) for pop in popLabels] if numCellSpks else [0] * len(popLabels)
totalSpikes = len(sel)
totalConnections = sum([len(conns) for conns in df['conns']])
numCells = len(cells)
firingRate = float(totalSpikes)/(numCells+numNetStims)/(timeRange[1]-timeRange[0])*1e3 if totalSpikes>0 else 0 # Calculate firing rate
connsPerCell = totalConnections/float(numCells) if numCells>0 else 0 # Calculate the number of connections per cell
if popRates:
avgRates = {}
tsecs = (timeRange[1]-timeRange[0])/1e3
for i, (pop, popNum) in enumerate(zip(popLabels, popNumCells)):
if numCells > 0 and pop != 'NetStims':
if numCellSpks == 0:
avgRates[pop] = 0
else:
avgRates[pop] = len([spkid for spkid in sel['spkind'].iloc[:numCellSpks-1] if df['pop'].iloc[int(spkid)]==pop])/popNum/tsecs
if numNetStims:
popNumCells[-1] = numNetStims
avgRates['NetStims'] = len([spkid for spkid in sel['spkind'].iloc[numCellSpks:]])/numNetStims/tsecs
if orderInverse:
y_range=(sel['spkind'].max(), sel['spkind'].min())
else:
y_range=(sel['spkind'].min(), sel['spkind'].max())
fig = figure(
title="Raster Plot",
tools=TOOLS,
active_drag = 'pan',
active_scroll = 'wheel_zoom',
tooltips=[('Cell GID', '@y'), ('Spike time', '@x')],
x_axis_label="Time (ms)",
y_axis_label=ylabelText,
x_range=(timeRange[0], timeRange[1]),
y_range=y_range,
toolbar_location='above')
t = Title()
if syncLines:
for spkt in sel['spkt'].tolist():
fig.line((spkt, spkt), (0, len(cells)+numNetStims), color='red', line_width=2)
print(syncMeasure())
t.text = 'cells=%i syns/cell=%0.1f rate=%0.1f Hz sync=%0.2f' % (numCells,connsPerCell,firingRate,syncMeasure())
else:
t.text = 'cells=%i syns/cell=%0.1f rate=%0.1f Hz' % (numCells,connsPerCell,firingRate)
fig.title = t
if spikeHist:
histo = np.histogram(sel['spkt'].tolist(), bins = np.arange(timeRange[0], timeRange[1], spikeHistBin))
histoT = histo[1][:-1]+spikeHistBin/2
histoCount = histo[0]
legendItems = []
grouped = sel.groupby('pop')
for name, group in grouped:
if popRates:
label = name + ' (%.3g Hz)'%(avgRates[name])
else:
label = name
fig.scatter(group['spkt'], group['spkind'], color=group['spkgidColor'], marker=marker, size=markerSize, legend_label=label)
#legendItems.append((label, [s]))
if spikeHist:
from bokeh.models import LinearAxis, Range1d
fig.extra_y_ranges={'spikeHist': Range1d(start=min(histoCount), end=max(histoCount))}
fig.add_layout(LinearAxis(y_range_name='spikeHist', axis_label='Spike count'), 'right')
fig.line (histoT, histoCount, line_width=2, y_range_name='spikeHist')
legend = Legend(items=legendItems, location=(10,0))
legend.click_policy='hide'
fig.add_layout(legend, 'right')
plot_layout = layout([fig], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="Raster Plot")
# save figure data
if saveData:
figData = {'spkTimes': sel['spkt'].tolist(), 'spkInds': sel['spkind'].tolist(), 'spkColors': sel['spkgidColor'].tolist(), 'cellGids': cellGids, 'sortedGids': df.index.tolist(), 'numNetStims': numNetStims,
'include': include, 'timeRange': timeRange, 'maxSpikes': maxSpikes, 'orderBy': orderBy, 'orderInverse': orderInverse, 'spikeHist': spikeHist,
'syncLines': syncLines}
_saveFigData(figData, saveData, 'raster')
# save figure
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename + '_iplot_raster.html'
file = open(filename, 'w')
file.write(html)
file.close()
if showFig: show(plot_layout)
return html, {'include': include, 'spkts': spkts, 'spkinds': sel['spkind'].tolist(), 'timeRange': timeRange}
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive dipole
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotDipole(expData={'label': 'Experiment', 'x':[], 'y':[]}, showFig=False):
'''
expData: experimental data; a dict with ['x'] and ['y'] 1-d vectors (either lists or np.arrays) of same length
showFig: show output figure in web browser (default: None)
'''
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
fig = figure(title="Dipole Plot", tools=TOOLS, toolbar_location='above', x_axis_label="Time (ms)", y_axis_label='Dipole (nAM x 3000)')
# renormalize the dipole and save
def baseline_renormalize():
# N_pyr cells in grid. This is PER LAYER
N_pyr = sim.cfg.N_pyr_x * sim.cfg.N_pyr_y
# dipole offset calculation: increasing number of pyr cells (L2 and L5, simultaneously)
# with no inputs resulted in an aggregate dipole over the interval [50., 1000.] ms that
# eventually plateaus at -48 fAm. The range over this interval is something like 3 fAm
# so the resultant correction is here, per dipole
# dpl_offset = N_pyr * 50.207
dpl_offset = {
# these values will be subtracted
'L2': N_pyr * 0.0443,
'L5': N_pyr * -49.0502
# 'L5': N_pyr * -48.3642,
# will be calculated next, this is a placeholder
# 'agg': None,
}
# L2 dipole offset can be roughly baseline shifted over the entire range of t
dpl = {'L2': np.array(sim.simData['dipole']['L2']),
'L5': np.array(sim.simData['dipole']['L5'])}
dpl['L2'] -= dpl_offset['L2']
# L5 dipole offset should be different for interval [50., 500.] and then it can be offset
# slope (m) and intercept (b) params for L5 dipole offset
# uncorrected for N_cells
# these values were fit over the range [37., 750.)
m = 3.4770508e-3
b = -51.231085
# these values were fit over the range [750., 5000]
t1 = 750.
m1 = 1.01e-4
b1 = -48.412078
t = sim.simData['t']
# piecewise normalization
dpl['L5'][t <= 37.] -= dpl_offset['L5']
dpl['L5'][(t > 37.) & (t < t1)] -= N_pyr * (m * t[(t > 37.) & (t < t1)] + b)
dpl['L5'][t >= t1] -= N_pyr * (m1 * t[t >= t1] + b1)
# recalculate the aggregate dipole based on the baseline normalized ones
dpl['agg'] = dpl['L2'] + dpl['L5']
return dpl
# convolve with a hamming window
def hammfilt(x, winsz):
from pylab import convolve
from numpy import hamming
win = hamming(winsz)
win /= sum(win)
return convolve(x,win,'same')
# baseline renormalize
dpl = baseline_renormalize()
# convert units from fAm to nAm, rescale and smooth
for key in dpl.keys():
dpl[key] *= 1e-6 * sim.cfg.dipole_scalefctr
if sim.cfg.dipole_smooth_win > 0:
dpl[key] = hammfilt(dpl[key], sim.cfg.dipole_smooth_win/sim.cfg.dt)
# Set index 0 to 0
dpl[key][0] = 0.0
# plot exp data
fig.line(expData['x'], expData['y'], color='black', legend=expData['label'])
# plot recorded dipole data
fig.line(sim.simData['t'], dpl['L2'], color='green', legend="L2Pyr", line_width=2.0)
fig.line(sim.simData['t'], dpl['L5'], color='red', legend="L5Pyr", line_width=2.0)
fig.line(sim.simData['t'], dpl['L2']+dpl['L5'], color='blue', legend="Aggregate", line_width=2.0)
fig.legend.location = "top_right"
fig.legend.click_policy = "hide"
plot_layout = layout(fig, sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="Dipole Plot")
if showFig:
show(fig)
return html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive dipole Spectrogram
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotDipoleSpectrogram(expData={'label': 'Experiment', 'x':[], 'y':[]}, minFreq = 1, maxFreq = 80, stepFreq = 1, norm = True, showFig=False):
'''
expData: experimental data; a dict with ['x'] and ['y'] 1-d vectors (either lists or np.arrays) of same length
showFig: show output figure in web browser (default: None)
'''
from .. import sim
from bokeh.plotting import figure, show
from bokeh.models import BasicTicker, ColorBar, ColumnDataSource, LinearColorMapper, PrintfTickFormatter
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout
# renormalize the dipole and save
def baseline_renormalize():
# N_pyr cells in grid. This is PER LAYER
N_pyr = sim.cfg.N_pyr_x * sim.cfg.N_pyr_y
# dipole offset calculation: increasing number of pyr cells (L2 and L5, simultaneously)
# with no inputs resulted in an aggregate dipole over the interval [50., 1000.] ms that
# eventually plateaus at -48 fAm. The range over this interval is something like 3 fAm
# so the resultant correction is here, per dipole
# dpl_offset = N_pyr * 50.207
dpl_offset = {
# these values will be subtracted
'L2': N_pyr * 0.0443,
'L5': N_pyr * -49.0502
# 'L5': N_pyr * -48.3642,
# will be calculated next, this is a placeholder
# 'agg': None,
}
# L2 dipole offset can be roughly baseline shifted over the entire range of t
dpl = {'L2': np.array(sim.simData['dipole']['L2']),
'L5': np.array(sim.simData['dipole']['L5'])}
dpl['L2'] -= dpl_offset['L2']
# L5 dipole offset should be different for interval [50., 500.] and then it can be offset
# slope (m) and intercept (b) params for L5 dipole offset
# uncorrected for N_cells
# these values were fit over the range [37., 750.)
m = 3.4770508e-3
b = -51.231085
# these values were fit over the range [750., 5000]
t1 = 750.
m1 = 1.01e-4
b1 = -48.412078
t = sim.simData['t']
# piecewise normalization
dpl['L5'][t <= 37.] -= dpl_offset['L5']
dpl['L5'][(t > 37.) & (t < t1)] -= N_pyr * (m * t[(t > 37.) & (t < t1)] + b)
dpl['L5'][t >= t1] -= N_pyr * (m1 * t[t >= t1] + b1)
# recalculate the aggregate dipole based on the baseline normalized ones
dpl['agg'] = dpl['L2'] + dpl['L5']
return dpl
# convolve with a hamming window
def hammfilt(x, winsz):
from pylab import convolve
from numpy import hamming
win = hamming(winsz)
win /= sum(win)
return convolve(x,win,'same')
# baseline renormalize
dpl = baseline_renormalize()
# convert units from fAm to nAm, rescale and smooth
for key in dpl.keys():
dpl[key] *= 1e-6 * sim.cfg.dipole_scalefctr
if sim.cfg.dipole_smooth_win > 0:
dpl[key] = hammfilt(dpl[key], sim.cfg.dipole_smooth_win/sim.cfg.dt)
# Set index 0 to 0
dpl[key][0] = 0.0
# plot recorded dipole data
#fig.line(sim.simData['t'], dpl['L2']+dpl['L5'], color='blue', legend="Aggregate", line_width=2.0)
# plot Morlet spectrogram
from ..support.morlet import MorletSpec, index2ms
spec = []
dplsum = dpl['L2'] + dpl['L5']
fs = int(1000.0 / sim.cfg.recordStep)
t_spec = np.linspace(0, index2ms(len(dplsum), fs), len(dplsum))
spec = MorletSpec(dplsum, fs, freqmin=minFreq, freqmax=maxFreq, freqstep=stepFreq)
f = np.array(range(minFreq, maxFreq+1, stepFreq)) # only used as output for user
vmin = np.array(spec.TFR).min()
vmax = np.array(spec.TFR).max()
T = [0, sim.cfg.duration]
F = spec.f
if norm:
spec.TFR = spec.TFR / vmax
S = spec.TFR
vc = [0, 1]
else:
S = spec.TFR
vc = [vmin, vmax]
# old code
fig = plt.imshow(S, extent=(np.amin(T), np.amax(T), np.amin(F), np.amax(F)), origin='lower', interpolation='None', aspect='auto', vmin=vc[0], vmax=vc[1],
cmap=plt.get_cmap('jet')) # viridis
plt.colorbar(label='Power')
plt.gca().invert_yaxis()
plt.ylabel('Hz')
plt.tight_layout()
if showFig:
plt.show()
# TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
# fig = figure(title="Dipole Spectrogram Plot", tools=TOOLS, toolbar_location='above',
# plot_width=S.shape[1], plot_height=S.shape[0], x_axis_label="Time (ms)", y_axis_label='Frequency (Hz)')
# fig.image(image=[S], x=[0], y=[0], dw=[S.shape[1]], dh=[S.shape[0]], palette='Spectral11')
# steps = list(range(0, S.shape[1], int(1./fs*1e6)))
# fig.xaxis.ticker = steps
# fig.xaxis.major_label_overrides = {t*40000/1e6: str(t) for t in T}
# fig.legend.location = "top_right"
# fig.legend.click_policy = "hide"
# plot_layout = layout(fig, sizing_mode='stretch_both')
# html = file_html(plot_layout, CDN, title="Dipole Spectrogram Plot")
# if showFig:
# show(fig)
return fig #html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive dipole Spectrogram
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotDipolePSD(expData={'label': 'Experiment', 'x':[], 'y':[]}, minFreq = 1, maxFreq = 80, stepFreq = 1, norm = True, showFig=False):
'''
expData: experimental data; a dict with ['x'] and ['y'] 1-d vectors (either lists or np.arrays) of same length
showFig: show output figure in web browser (default: None)
'''
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout
# renormalize the dipole and save
def baseline_renormalize():
# N_pyr cells in grid. This is PER LAYER
N_pyr = sim.cfg.N_pyr_x * sim.cfg.N_pyr_y
# dipole offset calculation: increasing number of pyr cells (L2 and L5, simultaneously)
# with no inputs resulted in an aggregate dipole over the interval [50., 1000.] ms that
# eventually plateaus at -48 fAm. The range over this interval is something like 3 fAm
# so the resultant correction is here, per dipole
# dpl_offset = N_pyr * 50.207
dpl_offset = {
# these values will be subtracted
'L2': N_pyr * 0.0443,
'L5': N_pyr * -49.0502
# 'L5': N_pyr * -48.3642,
# will be calculated next, this is a placeholder
# 'agg': None,
}
# L2 dipole offset can be roughly baseline shifted over the entire range of t
dpl = {'L2': np.array(sim.simData['dipole']['L2']),
'L5': np.array(sim.simData['dipole']['L5'])}
dpl['L2'] -= dpl_offset['L2']
# L5 dipole offset should be different for interval [50., 500.] and then it can be offset
# slope (m) and intercept (b) params for L5 dipole offset
# uncorrected for N_cells
# these values were fit over the range [37., 750.)
m = 3.4770508e-3
b = -51.231085
# these values were fit over the range [750., 5000]
t1 = 750.
m1 = 1.01e-4
b1 = -48.412078
t = sim.simData['t']
# piecewise normalization
dpl['L5'][t <= 37.] -= dpl_offset['L5']
dpl['L5'][(t > 37.) & (t < t1)] -= N_pyr * (m * t[(t > 37.) & (t < t1)] + b)
dpl['L5'][t >= t1] -= N_pyr * (m1 * t[t >= t1] + b1)
# recalculate the aggregate dipole based on the baseline normalized ones
dpl['agg'] = dpl['L2'] + dpl['L5']
return dpl
# convolve with a hamming window
def hammfilt(x, winsz):
from pylab import convolve
from numpy import hamming
win = hamming(winsz)
win /= sum(win)
return convolve(x,win,'same')
# baseline renormalize
dpl = baseline_renormalize()
# convert units from fAm to nAm, rescale and smooth
for key in dpl.keys():
dpl[key] *= 1e-6 * sim.cfg.dipole_scalefctr
if sim.cfg.dipole_smooth_win > 0:
dpl[key] = hammfilt(dpl[key], sim.cfg.dipole_smooth_win/sim.cfg.dt)
# Set index 0 to 0
dpl[key][0] = 0.0
# plot recorded dipole data
#fig.line(sim.simData['t'], dpl['L2']+dpl['L5'], color='blue', legend="Aggregate", line_width=2.0)
# plot Morlet spectrogram
from ..support.morlet import MorletSpec, index2ms
spec = []
dplsum = dpl['L2'] + dpl['L5']
fs = int(1000.0 / sim.cfg.recordStep)
t_spec = np.linspace(0, index2ms(len(dplsum), fs), len(dplsum))
spec = MorletSpec(dplsum, fs, freqmin=minFreq, freqmax=maxFreq, freqstep=stepFreq)
f = np.array(range(minFreq, maxFreq+1, stepFreq)) # only used as output for user
vmin = np.array(spec.TFR).min()
vmax = np.array(spec.TFR).max()
T = [0, sim.cfg.duration]
F = spec.f
S = spec.TFR
vc = [vmin, vmax]
Fs = int(1000.0/sim.cfg.recordStep)
signal = np.mean(S, 1)
# plot
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
fig = figure(title="Dipole PSD Plot", tools=TOOLS, toolbar_location='above', x_axis_label="Frequency (Hz)", y_axis_label='Power')
fig.line(F, signal, color='black', line_width=2.0)
fig.legend.location = "top_right"
fig.legend.click_policy = "hide"
plot_layout = layout(fig, sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="Dipole PSD Plot")
if showFig:
show(fig)
return html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive Spike Histogram
## ISSUES: Y scale, add colors to be effective
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotSpikeHist(include = ['allCells', 'eachPop'], legendLabels = [], timeRange = None, binSize = 5, overlay=True, yaxis = 'rate',
popColors=[], norm=False, smooth=None, filtFreq=False, filtOrder=3, saveData = None, saveFig = None, showFig = False):
'''
Plot spike histogram
- include (['all',|'allCells','allNetStims',|,120,|,'E1'|,('L2', 56)|,('L5',[4,5,6])]): List of data series to include.
Note: one line per item, not grouped (default: ['allCells', 'eachPop'])
- timeRange ([start:stop]): Time range of spikes shown; if None shows all (default: None)
- binSize (int): Size in ms of each bin (default: 5)
- overlay (True|False): Whether to overlay the data lines or plot in separate subplots (default: True)
- yaxis ('rate'|'count'): Units of y axis (firing rate in Hz, or spike count) (default: 'rate')
- popColors (dict): Dictionary with color (value) used for each population (key) (default: None)
- figSize ((width, height)): Size of figure (default: (10,8))
- saveData (None|True|'fileName'): File name where to save the final data used to generate the figure;
if set to True uses filename from simConfig (default: None)
- saveFig (None|True|'fileName'): File name where to save the figure;
if set to True uses filename from simConfig (default: None)
- showFig (True|False): Whether to show the figure or not (default: True)
- Returns figure handle
'''
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import gridplot
from bokeh.models import Legend
from bokeh.colors import RGB
print('Plotting interactive spike histogram...')
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
popColorDict=popColors.copy()
if popColorDict:
for pop, color in popColorDict.items():
if not isinstance(color, RGB):
popColorDict[pop] = RGB(*[round(f * 255) for f in color])
if timeRange is None:
timeRange = [0, sim.cfg.duration]
if 'eachPop' in include:
include.remove('eachPop')
for pop in sim.net.allPops: include.append(pop)
if yaxis == 'rate':
if norm:
yaxisLabel = 'Normalized firing rate'
else:
yaxisLabel = 'Avg cell firing rate (Hz)'
elif yaxis == 'count':
if norm:
yaxisLabel = 'Normalized spike count'
else:
yaxisLabel = 'Spike count'
figs=[]
if overlay:
figs.append(figure(title="Spike Historgram", tools=TOOLS, x_axis_label="Time (ms)", y_axis_label=yaxisLabel, toolbar_location='above'))
fig = figs[0]
legendItems = []
for iplot, subset in enumerate(include):
if not overlay:
figs.append(figure(title=str(subset), tools=TOOLS, x_axis_label="Time (ms)", y_axis_label=yaxisLabel))
fig = figs[iplot]
if isinstance(subset, list):
cells, cellGids, netStimLabels = getCellsInclude(subset)
else:
cells, cellGids, netStimLabels = getCellsInclude([subset])
numNetStims = 0
# Select cells to include
if len(cellGids) > 0:
try:
spkinds,spkts = list(zip(*[(spkgid,spkt) for spkgid,spkt in zip(sim.allSimData['spkid'],sim.allSimData['spkt']) if spkgid in cellGids]))
except:
spkinds,spkts = [],[]
else:
spkinds,spkts = [],[]
spkts, spkinds = list(spkts), list(spkinds)
numNetStims = 0
if 'stims' in sim.allSimData:
for netStimLabel in netStimLabels:
netStimSpks = [spk for cell,stims in sim.allSimData['stims'].items() \
for stimLabel,stimSpks in stims.items() for spk in stimSpks if stimLabel == netStimLabel]
if len(netStimSpks) > 0:
lastInd = max(spkinds) if len(spkinds)>0 else 0
spktsNew = netStimSpks
spkindsNew = [lastInd+1+i for i in range(len(netStimSpks))]
spkts.extend(spktsNew)
spkinds.extend(spkindsNew)
numNetStims += 1
histo = np.histogram(spkts, bins = np.arange(timeRange[0], timeRange[1], binSize))
histoT = histo[1][:-1]+binSize/2
histoCount = histo[0]
if yaxis=='rate':
histoCount = histoCount * (1000.0 / binSize) / (len(cellGids)+numNetStims) # convert to firing rate
if filtFreq:
from scipy import signal
fs = 1000.0/binSize
nyquist = fs/2.0
if isinstance(filtFreq, list): # bandpass
Wn = [filtFreq[0]/nyquist, filtFreq[1]/nyquist]
b, a = signal.butter(filtOrder, Wn, btype='bandpass')
elif isinstance(filtFreq, Number): # lowpass
Wn = filtFreq/nyquist
b, a = signal.butter(filtOrder, Wn)
histoCount = signal.filtfilt(b, a, histoCount)
if norm:
histoCount /= max(histoCount)
if smooth:
histoCount = _smooth1d(histoCount, smooth)[:len(histoT)]
if isinstance(subset, list):
color = colors[iplot%len(colors)]
else:
color = popColorDict[subset] if subset in popColorDict else colors[iplot%len(colors)]
label = legendLabels[iplot] if legendLabels else str(subset)
s = fig.line(histoT, histoCount, line_width=2.0, name=str(subset), color=color)
if overlay:
legendItems.append((str(label), [s]))
if overlay:
legend = Legend(items=legendItems, location=(10,0))
fig.add_layout(legend)
fig.legend.click_policy='hide'
fig.legend.location='top_right'
print(figs)
plot_layout = gridplot(figs, ncols=1, merge_tools=False, sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="Spike Historgram")
if showFig: show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename+'_'+'ispikeHist.html'
file = open(filename, 'w')
file.write(html)
file.close()
return html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive Rate PSD
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotRatePSD(include = ['allCells', 'eachPop'], timeRange = None, binSize = 5, maxFreq = 100, NFFT = 256, noverlap = 128, smooth = 0, overlay=True, ylim = None,
popColors = {}, figSize=(1000, 400), saveData = None, saveFig = None, showFig = False):
'''
Plot firing rate power spectral density (PSD)
- include (['all',|'allCells','allNetStims',|,120,|,'E1'|,('L2', 56)|,('L5',[4,5,6])]): List of data series to include.
Note: one line per item, not grouped (default: ['allCells', 'eachPop'])
- timeRange ([start:stop]): Time range of spikes shown; if None shows all (default: None)
- binSize (int): Size in ms of spike bins (default: 5)
- maxFreq (float): Maximum frequency to show in plot (default: 100)
- NFFT (float): The number of data points used in each block for the FFT (power of 2) (default: 256)
- smooth (int): Window size for smoothing; no smoothing if 0 (default: 0)
- overlay (True|False): Whether to overlay the data lines or plot in separate subplots (default: True)
- graphType ('line'|'bar'): Type of graph to use (line graph or bar plot) (default: 'line')
- yaxis ('rate'|'count'): Units of y axis (firing rate in Hz, or spike count) (default: 'rate')
- popColors (dict): Dictionary with color (value) used for each population (key) (default: None)
- figSize ((width, height)): Size of figure (default: (10,8))
- saveData (None|True|'fileName'): File name where to save the final data used to generate the figure;
if set to True uses filename from simConfig (default: None)
- saveFig (None|True|'fileName'): File name where to save the figure;
if set to True uses filename from simConfig (default: None)
- showFig (True|False): Whether to show the figure or not (default: True)
- Returns figure handle
'''
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout
from bokeh.colors import RGB
from bokeh.models import Legend
print('Plotting interactive firing rate power spectral density (PSD) ...')
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
popColorDict=popColors.copy()
if popColorDict:
for pop, color in popColorDict.items():
popColorDict[pop] = RGB(*[round(f * 255) for f in color])
# Replace 'eachPop' with list of pops
if 'eachPop' in include:
include.remove('eachPop')
for pop in sim.net.allPops: include.append(pop)
# time range
if timeRange is None:
timeRange = [0,sim.cfg.duration]
histData = []
allPower, allSignal, allFreqs = [], [], []
legendItems = []
figs=[]
if overlay:
figs.append(figure(title="PSD Rate Plot", tools=TOOLS, x_axis_label="Frequncy (Hz)", y_axis_label="Power Spectral Density (db/Hz)", toolbar_location='above'))
fig = figs[0]
legendItems = []
# Plot separate line for each entry in include
for iplot, subset in enumerate(include):
if not overlay:
figs.append(figure(title=str(subset), tools=TOOLS, x_axis_label="Frequency (HZ)", y_axis_label="Power Spectral Density (db/Hz)"))
fig = figs[iplot]
cells, cellGids, netStimLabels = getCellsInclude([subset])
numNetStims = 0
# Select cells to include
if len(cellGids) > 0:
try:
spkinds,spkts = list(zip(*[(spkgid,spkt) for spkgid,spkt in zip(sim.allSimData['spkid'],sim.allSimData['spkt']) if spkgid in cellGids]))
except:
spkinds,spkts = [],[]
else:
spkinds,spkts = [],[]
# Add NetStim spikes
spkts, spkinds = list(spkts), list(spkinds)
numNetStims = 0
if 'stims' in sim.allSimData:
for netStimLabel in netStimLabels:
netStimSpks = [spk for cell,stims in sim.allSimData['stims'].items() \
for stimLabel,stimSpks in stims.items() for spk in stimSpks if stimLabel == netStimLabel]
if len(netStimSpks) > 0:
lastInd = max(spkinds) if len(spkinds)>0 else 0
spktsNew = netStimSpks
spkindsNew = [lastInd+1+i for i in range(len(netStimSpks))]
spkts.extend(spktsNew)
spkinds.extend(spkindsNew)
numNetStims += 1
histo = np.histogram(spkts, bins = np.arange(timeRange[0], timeRange[1], binSize))
histoT = histo[1][:-1]+binSize/2
histoCount = histo[0]
histoCount = histoCount * (1000.0 / binSize) / (len(cellGids)+numNetStims) # convert to rates
histData.append(histoCount)
Fs = 1000.0 / binSize
power = mlab.psd(histoCount, Fs=Fs, NFFT=NFFT, detrend=mlab.detrend_none, window=mlab.window_hanning,
noverlap=noverlap, pad_to=None, sides='default', scale_by_freq=None)
if smooth:
signal = _smooth1d(10*np.log10(power[0]), smooth)
else:
signal = 10*np.log10(power[0])
freqs = power[1]
color = popColorDict[subset] if subset in popColorDict else colors[iplot%len(colors)]
allFreqs.append(freqs)
allPower.append(power)
allSignal.append(signal)
s = fig.line(freqs[freqs<maxFreq], signal[freqs<maxFreq], line_width = 2.0, color=color)
if overlay: legendItems.append((subset, [s]))
if overlay:
legend = Legend(items=legendItems)
legend.click_policy = 'hide'
fig.add_layout(legend)
plot_layout = layout(figs, ncols=1, plot_width=figSize[0], figHeight=figSize[1], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="PSD Rate Plot")
if showFig: show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename+'_'+'ipsd.html'
file = open(filename, 'w')
file.write(html)
file.close()
return html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive Traces
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotTraces(include=None, timeRange=None, overlay=False, oneFigPer='cell', rerun=False, colors=None, ylim=None, axis='on', fontSize=12, figSize=(10,8), saveData=None, saveFig=None, showFig=True, ylabel=None, linkAxes=False):
'''
Plot recorded traces
- include (['all',|'allCells','allNetStims',|,120,|,'E1'|,('L2', 56)|,('L5',[4,5,6])]): List of cells for which to plot
the recorded traces (default: [])
- timeRange ([start:stop]): Time range of spikes shown; if None shows all (default: None)
- overlay (True|False): Whether to overlay the data lines or plot in separate subplots (default: False)
- oneFigPer ('cell'|'trace'): Whether to plot one figure per cell (showing multiple traces)
or per trace (showing multiple cells) (default: 'cell')
- rerun (True|False): rerun simulation so new set of cells gets recorded (default: False)
- colors (list): List of normalized RGB colors to use for traces
- ylim (list): Y-axis limits
- axis ('on'|'off'): Whether to show axis or not; if not, then a scalebar is included (default: 'on')
- figSize ((width, height)): Size of figure (default: (10,8))
- saveData (None|True|'fileName'): File name where to save the final data used to generate the figure;
if set to True uses filename from simConfig (default: None)
- saveFig (None|True|'fileName'): File name where to save the figure;
if set to True uses filename from simConfig (default: None)
- showFig (True|False): Whether to show the figure or not (default: True)
- Returns figure handles
'''
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout, column
from bokeh.models import HoverTool
from bokeh.models import Legend
from bokeh.colors import RGB
print('Plotting interactive recorded cell traces per', oneFigPer)
TOOLS = 'save,pan,box_zoom,reset,wheel_zoom'
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
if include is None: # if none, record from whatever was recorded
if 'plotTraces' in sim.cfg.analysis and 'include' in sim.cfg.analysis['plotTraces']:
include = sim.cfg.analysis['plotTraces']['include'] + sim.cfg.recordCells
else:
include = sim.cfg.recordCells
if timeRange is None:
timeRange = [0,sim.cfg.duration]
tracesList = list(sim.cfg.recordTraces.keys())
tracesList.sort()
cells, cellGids, _ = getCellsInclude(include)
gidPops = {cell['gid']: cell['tags']['pop'] for cell in cells}
recordStep = sim.cfg.recordStep
figs = {}
tracesData = []
if ylabel is None:
y_axis_label = None
else:
y_axis_label = ylabel
if oneFigPer == 'cell':
for gid in cellGids:
if overlay:
figs['_gid_' + str(gid)] = figure(title = "Cell {}, Pop {}".format(gid, gidPops[gid]),
tools = TOOLS,
active_drag = 'pan',
active_scroll = 'wheel_zoom',
x_axis_label="Time (ms)",
y_axis_label=y_axis_label
)
fig = figs['_gid_' + str(gid)]
else:
figs['_gid_' + str(gid)] = []
for itrace, trace in enumerate(tracesList):
if 'cell_{}'.format(gid) in sim.allSimData[trace]:
fullTrace = sim.allSimData[trace]['cell_{}'.format(gid)]
if isinstance(fullTrace, dict):
data = [fullTrace[key][int(timeRange[0]/recordStep):int(timeRange[1]/recordStep)] for key in list(fullTrace.keys())]
lenData = len(data[0])
data = np.transpose(np.array(data))
else:
data = fullTrace[int(timeRange[0]/recordStep):int(timeRange[1]/recordStep)]
lenData = len(data)
t = np.arange(timeRange[0], timeRange[1]+recordStep, recordStep)
tracesData.append({'t': t, 'cell_'+str(gid)+'_'+trace: data})
if overlay:
fig.line(t[:lenData], data, line_width=2, line_color=colors[itrace], legend_label=trace)
hover = HoverTool(tooltips=[('Time', '@x'), ('Value', '@y')], mode='vline')
fig.add_tools(hover)
fig.legend.click_policy="hide"
else:
subfig = figure(title = "Cell {}, Pop {}".format(gid, gidPops[gid]),
tools = TOOLS,
active_drag = 'pan',
active_scroll = 'wheel_zoom',
x_axis_label="Time (ms)",
y_axis_label=y_axis_label
)
subfig.line(t[:lenData], data, line_width=2, line_color=colors[itrace], legend_label=trace)
if linkAxes:
if itrace > 0:
subfig.x_range = figs['_gid_' + str(gid)][0].x_range
subfig.y_range = figs['_gid_' + str(gid)][0].y_range
hover = HoverTool(tooltips=[('Time', '@x'), ('Value', '@y')], mode='vline')
subfig.add_tools(hover)
subfig.legend.click_policy="hide"
figs['_gid_' + str(gid)].append(subfig)
elif oneFigPer == 'trace':
for itrace, trace in enumerate(tracesList):
if overlay:
figs['_trace_' + str(trace)] = figure(title = str(trace),
tools = TOOLS,
active_drag = 'pan',
active_scroll = 'wheel_zoom',
x_axis_label="Time (ms)",
y_axis_label=y_axis_label
)
fig = figs['_trace_' + str(trace)]
else:
figs['_trace_' + str(trace)] = []
for igid, gid in enumerate(cellGids):
if 'cell_'+str(gid) in sim.allSimData[trace]:
fullTrace = sim.allSimData[trace]['cell_'+str(gid)]
if isinstance(fullTrace, dict):
data = [fullTrace[key][int(timeRange[0]/recordStep):int(timeRange[1]/recordStep)] for key in list(fullTrace.keys())]
lenData = len(data[0])
data = np.transpose(np.array(data))
else:
data = np.array(fullTrace[int(timeRange[0]/recordStep):int(timeRange[1]/recordStep)])
lenData = len(data)
t = np.arange(timeRange[0], timeRange[1]+recordStep, recordStep)
tracesData.append({'t': t, 'cell_'+str(gid)+'_'+trace: data})
if overlay:
fig.line(t[:lenData], data, line_width=2, line_color=colors[igid], legend_label='Cell %d, Pop %s '%(int(gid), gidPops[gid]))
hover = HoverTool(tooltips=[('Time', '@x'), ('Value', '@y')], mode='vline')
fig.add_tools(hover)
fig.legend.click_policy="hide"
else:
subfig = figure(title = str(trace),
tools = TOOLS,
active_drag = 'pan',
active_scroll = 'wheel_zoom',
x_axis_label="Time (ms)",
y_axis_label=y_axis_label
)
if isinstance(data, list):
for tl,dl in zip(t,data):
subfig.line(tl[:len(dl)], dl, line_width=1.5, line_color=colors[igid], legend_label='Cell %d, Pop %s '%(int(gid), gidPops[gid]))
else:
subfig.line(t[:len(data)], data, line_width=1.5, line_color=colors[igid], legend_label='Cell %d, Pop %s '%(int(gid), gidPops[gid]))
if linkAxes:
if igid > 0:
subfig.x_range = figs['_trace_' + str(trace)][0].x_range
subfig.y_range = figs['_trace_' + str(trace)][0].y_range
hover = HoverTool(tooltips=[('Time', '@x'), ('Value', '@y')], mode='vline')
subfig.add_tools(hover)
subfig.legend.click_policy="hide"
figs['_trace_' + str(trace)].append(subfig)
for figLabel, figObj in figs.items():
if overlay:
plot_layout = layout(figObj, sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title=figLabel)
overlay_text = '_overlay'
else:
plot_layout = column(*figObj, sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title=figLabel)
overlay_text = ''
if showFig: show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename + '_iplot_traces' + figLabel + overlay_text + '.html'
file = open(filename, 'w')
file.write(html)
file.close()
return html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive LFP
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotLFP(electrodes = ['avg', 'all'], plots = ['timeSeries', 'PSD', 'spectrogram', 'locations'], timeRange = None, NFFT = 256, noverlap = 128,
nperseg = 256, maxFreq = 100, smooth = 0, separation = 1.0, includeAxon=True, logx=False, logy=False, norm=False, dpi = 200, overlay=False, filtFreq = False, filtOrder=3, detrend=False, colors = None, saveData = None, saveFig = None, showFig = False):
from .. import sim
from bokeh.plotting import figure, show
from bokeh.resources import CDN
from bokeh.embed import file_html
from bokeh.layouts import layout, column, row
from bokeh.colors import RGB
print('Plotting interactive LFP ...')
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
if timeRange is None:
timeRange = [0,sim.cfg.duration]
lfp = np.array(sim.allSimData['LFP'])[int(timeRange[0]/sim.cfg.recordStep):int(timeRange[1]/sim.cfg.recordStep),:]
if not colors:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
# electrode selection
if 'all' in electrodes:
electrodes.remove('all')
electrodes.extend(list(range(int(sim.net.recXElectrode.nsites))))
figs = {}
data = {'lfp': lfp}
# time series plot
# TODO add scalebar
if 'timeSeries' in plots:
figs['timeSeries'] = figure(title="LFP Time Series Plot", tools=TOOLS, x_axis_label="Time (ms)", y_axis_label="LFP electrode")
figs['timeSeries'].yaxis.major_tick_line_color = None
figs['timeSeries'].yaxis.minor_tick_line_color = None
figs['timeSeries'].yaxis.major_label_text_font_size = '0pt'
figs['timeSeries'].yaxis.axis_line_color = None
figs['timeSeries'].ygrid.visible = False
ydisp = np.absolute(lfp).max() * separation
offset = 1.0*ydisp
t = np.arange(timeRange[0], timeRange[1], sim.cfg.recordStep)
for i,elec in enumerate(electrodes[::-1]):
if elec == 'avg':
dplSum = np.mean(lfp, axis=1)
color = 'black'
lw=1.0
# figs['timeSeries'].line(t, -lfpPlot+(i*ydisp), line_color=color)
elif isinstance(elec, Number) and elec <= sim.net.recXElectrode.nsites:
pass
lfpPlot = lfp[:, elec]
color = colors[i%len(colors)]
lw=1.0
legend=str(elec)
# if len(electrodes) > 1:
# legend=str(elec)
# else:
# legend=None
figs['timeSeries'].line(t, lfpPlot+(i*ydisp), color=color, name=str(elec), legend=legend)
data['lfpPlot'] = lfpPlot
data['ydisp'] = ydisp
data['t'] = t
figs['timeSeries'].legend.location = "top_right"
figs['timeSeries'].legend.click_policy = "hide"
plot_layout = layout(figs['timeSeries'], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="Time Series LFP Plot")
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename+'_'+'lfp_timeseries.png'
file = open(filename, 'w')
file.write(html)
file.close()
# PSD ----------------------------------
# TODO fix final layout of plots
if 'PSD' in plots:
figs['psd'] = []
allFreqs = []
allSignal = []
data['allFreqs'] = allFreqs
data['allSignal'] = allSignal
for i,elec in enumerate(electrodes):
p = figure(title="Electrode {}".format(str(elec)), tools=TOOLS, x_axis_label="Frequency (Hz)", y_axis_label="db/Hz")
if elec == 'avg':
lfpPlot = np.mean(lfp, axis=1)
color = 'black'
lw=1.5
elif isinstance(elec, Number) and elec <= sim.net.recXElectrode.nsites:
lfpPlot = lfp[:, elec]
color = colors[i%len(colors)]
lw=1.5
Fs = int(1000.0/sim.cfg.recordStep)
power = mlab.psd(lfpPlot, Fs=Fs, NFFT=NFFT, detrend=mlab.detrend_none, window=mlab.window_hanning,
noverlap=noverlap, pad_to=None, sides='default', scale_by_freq=None)
if smooth:
signal = _smooth1d(10*np.log10(power[0]), smooth)
else:
signal = 10*np.log10(power[0])
freqs = power[1]
allFreqs.append(freqs)
allSignal.append(signal)
p.line(freqs[freqs<maxFreq], signal[freqs<maxFreq], color=color)
figs['psd'].append(p)
# format plot
plot_layout = column(figs['psd'] )
html = file_html(plot_layout, CDN, title="LFP Power Spectral Density")
show(plot_layout)
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename+'_'+'lfp_psd.png'
file = open(filename, 'w')
file.write(html)
file.close()
# Spectrogram ------------------------------
# TODO impprove the color mapper to be more detailed
if 'spectrogram' in plots:
import matplotlib.cm as cm
from bokeh.transform import linear_cmap
from bokeh.models import ColorBar
from scipy import signal as spsig
# numCols = np.round(len(electrodes) / maxPlots) + 1
figs['spectro'] = []
#t = np.arange(timeRange[0], timeRange[1], sim.cfg.recordStep)
logx_spec = []
for i,elec in enumerate(electrodes):
if elec == 'avg':
lfpPlot = np.mean(lfp, axis=1)
elif isinstance(elec, Number) and elec <= sim.net.recXElectrode.nsites:
lfpPlot = lfp[:, elec]
# creates spectrogram over a range of data
# from: http://joelyancey.com/lfp-python-practice/
fs = int(1000.0/sim.cfg.recordStep)
f, t_spec, x_spec = spsig.spectrogram(lfpPlot, fs=fs, window='hanning',
detrend=mlab.detrend_none, nperseg=nperseg, noverlap=noverlap, nfft=NFFT, mode='psd')
x_mesh, y_mesh = np.meshgrid(t_spec*1000.0, f[f<maxFreq])
logx_spec.append(10*np.log10(x_spec[f<maxFreq]))
vmin = np.array(logx_spec).min()
vmax = np.array(logx_spec).max()
for i,elec in enumerate(electrodes):
p = figure(title="Electrode {}".format(str(elec)), tools=TOOLS, x_range=(0, timeRange[1]), y_range=(0, maxFreq),
x_axis_label = "Time (ms)", y_axis_label = "Frequency(Hz)")
mapper = linear_cmap (field_name='dB/Hz', palette='Spectral11', low=vmin, high=vmax)
color_bar = ColorBar(color_mapper=mapper['transform'], width=8, location=(0,0), label_standoff=7, major_tick_line_color=None)
p.image(image=[x_mesh, y_mesh, logx_spec[i]], x=0, y=0, color_mapper=mapper['transform'], dw=timeRange[1], dh=100)
p.add_layout(color_bar, 'right')
figs['spectro'].append(p)
plot_layout = column(figs['spectro'] )
html = file_html(plot_layout, CDN, title="LFP Power Spectral Density")
show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename+'_'+'lfp_spectrogram.png'
file = open(filename, 'w')
file.write(html)
file.close()
return html
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive connectivity
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotConn(includePre=['all'], includePost=['all'], feature='strength', orderBy='gid', figSize=(10,10), groupBy='pop', groupByIntervalPre=None, groupByIntervalPost=None, removeWeightNorm=False, graphType='matrix', synOrConn='syn', synMech=None, connsFile=None, tagsFile=None, clim=None, fontSize=12, saveData=None, saveFig=None, showFig=False):
'''
Plot network connectivity
- includePre (['all',|'allCells','allNetStims',|,120,|,'E1'|,('L2', 56)|,('L5',[4,5,6])]): Cells to show (default: ['all'])
- includePost (['all',|'allCells','allNetStims',|,120,|,'E1'|,('L2', 56)|,('L5',[4,5,6])]): Cells to show (default: ['all'])
- feature ('weight'|'delay'|'numConns'|'probability'|'strength'|'convergence'|'divergence'): Feature to show in connectivity matrix;
the only features applicable to groupBy='cell' are 'weight', 'delay' and 'numConns'; 'strength' = weight * probability (default: 'strength')
- groupBy ('pop'|'cell'|'y'|: Show matrix for individual cells, populations, or by other numeric tag such as 'y' (default: 'pop')
- groupByInterval (int or float): Interval of groupBy feature to group cells by in conn matrix, e.g. 100 to group by cortical depth in steps of 100 um (default: None)
- orderBy ('gid'|'y'|'ynorm'|...): Unique numeric cell property to order x and y axes by, e.g. 'gid', 'ynorm', 'y' (requires groupBy='cells') (default: 'gid')
- graphType ('matrix','bar','pie'): Type of graph to represent data (default: 'matrix')
- synOrConn ('syn'|'conn'): Use synapses or connections; note 1 connection can have multiple synapses (default: 'syn')
- figSize ((width, height)): Size of figure (default: (10,10))
- synMech (['AMPA', 'GABAA',...]): Show results only for these syn mechs (default: None)
- saveData (None|True|'fileName'): File name where to save the final data used to generate the figure;
if set to True uses filename from simConfig (default: None)
- saveFig (None|True|'fileName'): File name where to save the figure;
if set to True uses filename from simConfig (default: None)
- showFig (True|False): Whether to show the figure or not (default: True)
- Returns figure handles
'''
from netpyne import sim
from netpyne.analysis import network
from bokeh.plotting import figure, show
from bokeh.transform import linear_cmap
from bokeh.palettes import Viridis256
from bokeh.palettes import viridis
from bokeh.models import ColorBar
from bokeh.embed import file_html
from bokeh.resources import CDN
from bokeh.layouts import layout
from bokeh.colors import RGB
print('Plotting interactive connectivity matrix...')
if connsFile and tagsFile:
connMatrix, pre, post = network._plotConnCalculateFromFile(includePre, includePost, feature, orderBy, groupBy, groupByIntervalPre, groupByIntervalPost, synOrConn, synMech, connsFile, tagsFile, removeWeightNorm)
else:
connMatrix, pre, post = network._plotConnCalculateFromSim(includePre, includePost, feature, orderBy, groupBy, groupByIntervalPre, groupByIntervalPost, synOrConn, synMech, removeWeightNorm)
if connMatrix is None:
print(" Error calculating connMatrix in iplotConn()")
return None
# TODO: set plot font size in Bokeh
# for groupBy = 'cell' (needed to properly format data for Bokeh)
if groupBy == 'cell':
pre = [str(cell['gid']) for cell in pre]
post = [str(cell['gid']) for cell in post]
# matrix plot
if graphType == 'matrix':
pandas_data = pd.DataFrame(data=connMatrix, index=pre, columns=post)
pandas_data.index.name = 'pre'
pandas_data.columns.name = 'post'
bokeh_data = pandas_data.reset_index().melt(id_vars=['pre'], value_name=feature)
# the colormapper doesn't work if both high and low values are the same
low = np.nanmin(connMatrix)
high = np.nanmax(connMatrix)
if low == high:
low = low - 0.1 * low
high = high + 0.1 * high
conn_colormapper = linear_cmap(
field_name=feature,
palette=Viridis256,
low=low,
high=high,
nan_color='white'
)
conn_colorbar = ColorBar(color_mapper=conn_colormapper['transform'])
conn_colorbar.location = (0, 0)
conn_colorbar.title = feature
fig = figure(
x_range=pandas_data.columns.values,
y_range=np.flip(pandas_data.index.values),
tools = 'hover,save,pan,box_zoom,reset,wheel_zoom',
active_drag = 'pan',
active_scroll = 'wheel_zoom',
tooltips=[('Pre', '@pre'), ('Post', '@post'), (feature, '@' + feature)],
title='Connection ' + feature + ' matrix',
toolbar_location='below',
x_axis_location='above'
)
fig.xaxis.axis_label = pandas_data.columns.name
fig.yaxis.axis_label = pandas_data.index.name
fig.yaxis.major_label_orientation = 'horizontal'
fig.xaxis.major_label_orientation = 'vertical'
fig.rect(
source=bokeh_data,
x='post',
y='pre',
width=1,
height=1,
color=conn_colormapper
)
fig.add_layout(conn_colorbar, 'right')
# TODO: add grid lines?
# stacked bar graph
elif graphType == 'bar':
if groupBy == 'pop':
popsPre, popsPost = pre, post
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
bar_colors = colors[0:len(popsPre)]
data = {'post' : popsPost}
for popIndex, pop in enumerate(popsPre):
data[pop] = connMatrix[popIndex, :]
fig = figure(
x_range=popsPost,
title='Connection ' + feature + ' stacked bar graph',
toolbar_location=None,
tools='hover,save,pan,box_zoom,reset,wheel_zoom',
active_drag='pan',
active_scroll = 'wheel_zoom',
tooltips=[('Pre', '$name'), ('Post', '@post'), (feature, '@$name')],
)
fig.vbar_stack(
popsPre,
x='post',
width=0.9,
color=bar_colors,
source=data,
legend_label=popsPre,
)
fig.xgrid.grid_line_color = None
fig.legend.title = 'Pre'
fig.xaxis.axis_label = 'Post'
fig.yaxis.axis_label = feature
elif groupBy == 'cell':
print(' Error: plotConn graphType="bar" with groupBy="cell" not yet implemented')
return None
elif graphType == 'pie':
print(' Error: plotConn graphType="pie" not yet implemented')
return None
plot_layout = layout([fig], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title='Connection ' + feature + ' matrix')
#save figure data
if saveData:
figData = {'connMatrix': connMatrix, 'feature': feature, 'groupBy': groupBy,
'includePre': includePre, 'includePost': includePost, 'saveData': saveData, 'saveFig': saveFig, 'showFig': showFig}
_saveFigData(figData, saveData, 'conn')
# save figure
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename+'_iplot_conn_'+groupBy+'_'+feature+'_'+graphType+'.html'
file = open(filename, 'w')
file.write(html)
file.close()
# show fig
if showFig: show(fig)
return html