"""
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__
try:
basestring
except NameError:
basestring = str
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, _guiTheme
import numpy as np
import pandas as pd
from bokeh.themes import built_in_themes
from bokeh.io import curdoc
from bokeh.palettes import Viridis256
from bokeh.models import HoverTool
bokeh_theme = curdoc().theme
# -------------------------------------------------------------------------------------------------------------------
## 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, **kwargs):
"""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 ...')
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
TOOLS = 'hover,save,pan,box_zoom,reset,wheel_zoom',
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
else:
colors = kwargs['palette']
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=[] if include == ['allCells'] else cellGids, timeRange=timeRange)
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 = None,
active_scroll = None,
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, **kwargs):
"""
iplotDipole
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
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
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, **kwargs):
"""
iplotDipoleSpectrogram
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
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
# 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, **kwargs):
"""
iplotDipolePSD
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
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
# 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
# -------------------------------------------------------------------------------------------------------------------
[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, **kwargs):
"""
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...')
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
else:
colors = kwargs['palette']
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 Histogram", 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 Histogram")
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={}, saveData=None, saveFig=None, showFig=False, **kwargs):
"""
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) ...')
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
else:
colors = kwargs['palette']
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="Frequency (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])
# 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, 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, **kwargs):
"""
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)
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
TOOLS = 'save,pan,box_zoom,reset,wheel_zoom'
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
else:
colors = kwargs['palette']
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 = None,
active_scroll = None,
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 = None,
active_scroll = None,
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 = None,
active_scroll = None,
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 = None,
active_scroll = None,
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'], timeRange=None, NFFT=256, noverlap=128, nperseg=256, minFreq=1, maxFreq=100, stepFreq=1, smooth=0, separation=1.0, includeAxon=True, logx=False, logy=False, normSignal=False, normPSD=False, normSpec=False, overlay=False, filtFreq=False, filtOrder=3, detrend=False, transformMethod='morlet', colors=None, saveData=None, saveFig=None, showFig=False, **kwargs):
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 ...')
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList]
else:
colors = kwargs['palette']
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),:]
# 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,
active_drag = None,
active_scroll = None,
x_axis_label="Time (ms)",
y_axis_label="LFP electrode",
toolbar_location="above")
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
hover = HoverTool(tooltips=[('Time', '@x'), ('Value', '@y')], mode='vline')
figs['timeSeries'].add_tools(hover)
figs['timeSeries'].legend.click_policy="hide"
ydisp = np.absolute(lfp).max() * separation
offset = 1.0*ydisp
t = np.arange(timeRange[0], timeRange[1], sim.cfg.recordStep)
avg_color = 'black'
if 'theme' in kwargs:
if kwargs['theme'] == 'gui' or 'dark' in kwargs['theme']:
avg_color = 'white'
for i,elec in enumerate(electrodes[::-1]):
if elec == 'avg':
dplSum = np.mean(lfp, axis=1)
color = avg_color
elif isinstance(elec, Number) and elec <= sim.net.recXElectrode.nsites:
pass
lfpPlot = lfp[:, elec]
color = colors[i%len(colors)]
legend=str(elec)
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 showFig: show(plot_layout)
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename + '_iLFP_timeseries.html'
file = open(filename, 'w')
file.write(html)
file.close()
# PSD ----------------------------------
if 'PSD' in plots:
figs['psd'] = []
allFreqs = []
allSignal = []
data['allFreqs'] = allFreqs
data['allSignal'] = allSignal
avg_color = 'black'
if 'theme' in kwargs:
if kwargs['theme'] == 'gui' or 'dark' in kwargs['theme']:
avg_color = 'white'
for i,elec in enumerate(electrodes):
p = figure(title="Electrode {}".format(str(elec)),
tools=TOOLS,
active_drag = None,
active_scroll = None,
x_axis_label="Frequency (Hz)",
y_axis_label="db/Hz",
toolbar_location="above")
hover = HoverTool(tooltips=[('Frequency', '@x'), ('Power', '@y')], mode='vline')
p.add_tools(hover)
p.legend.click_policy="hide"
if elec == 'avg':
lfpPlot = np.mean(lfp, axis=1)
color = avg_color
elif isinstance(elec, Number) and elec <= sim.net.recXElectrode.nsites:
lfpPlot = lfp[:, elec]
color = colors[i%len(colors)]
# Morlet wavelet transform method
if transformMethod == 'morlet':
from ..support.morlet import MorletSpec, index2ms
Fs = int(1000.0/sim.cfg.recordStep)
morletSpec = MorletSpec(lfpPlot, Fs, freqmin=minFreq, freqmax=maxFreq, freqstep=stepFreq)
freqs = F = morletSpec.f
spec = morletSpec.TFR
signal = np.mean(spec, 1)
ylabel = 'Power'
# FFT transform method
elif transformMethod == 'fft':
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]
ylabel = 'Power (dB/Hz)'
allFreqs.append(freqs)
allSignal.append(signal)
p.line(freqs[freqs<maxFreq], signal[freqs<maxFreq], color=color, line_width=3)
figs['psd'].append(p)
plot_layout = column(figs['psd'], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="LFP Power Spectral Density")
if showFig: show(plot_layout)
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename+'_iLFP_psd.html'
file = open(filename, 'w')
file.write(html)
file.close()
# Spectrogram ------------------------------
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
figs['spectro'] = []
logx_spec = []
# Morlet wavelet transform method
if transformMethod == 'morlet':
from ..support.morlet import MorletSpec, index2ms
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]
fs = int(1000.0/sim.cfg.recordStep)
t_spec = np.linspace(0, index2ms(len(lfpPlot), fs), len(lfpPlot))
spec.append(MorletSpec(lfpPlot, fs, freqmin=minFreq, freqmax=maxFreq, freqstep=stepFreq))
f = np.array(range(minFreq, maxFreq+1, stepFreq))
x_mesh, y_mesh = np.meshgrid(t_spec*1000.0, f[f<maxFreq])
vmin = np.array([s.TFR for s in spec]).min()
vmax = np.array([s.TFR for s in spec]).max()
for i,elec in enumerate(electrodes):
T = timeRange
F = spec[i].f
if normSpec:
spec[i].TFR = spec[i].TFR / vmax
S = spec[i].TFR
vc = [0, 1]
else:
S = spec[i].TFR
vc = [vmin, vmax]
p = figure(
title="Electrode {}".format(str(elec)),
tools=TOOLS,
active_drag = None,
active_scroll = None,
x_range=(0, timeRange[1]),
y_range=(0, maxFreq),
x_axis_label = "Time (ms)",
y_axis_label = "Frequency(Hz)",
toolbar_location="above",
#tooltips = [("Time", "$x"), ("Frequency", "$y"), ("Power", "@image")],
)
mapper = linear_cmap(field_name='dB/Hz', palette=Viridis256, low=vmin, high=vmax)
color_bar = ColorBar(color_mapper=mapper['transform'], location=(0,0), label_standoff=15, title='Power')
p.image(image=[x_mesh, y_mesh, S], x=0, y=0, color_mapper=mapper['transform'], dw=timeRange[1], dh=100)
p.add_layout(color_bar, 'right')
figs['spectro'].append(p)
# FFT transform method
elif transformMethod == 'fft':
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,
active_drag = None,
active_scroll = None,
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=Viridis256, 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'], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="LFP Power Spectral Density")
if showFig: show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename + '_iLFP_timefreq.html'
file = open(filename, 'w')
file.write(html)
file.close()
# -------------------------------------------------------------------------------------------------------------------
## 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, **kwargs):
"""
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 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
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 = None,
active_scroll = None,
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
colors = colorList
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= None,
active_scroll = None,
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
# -------------------------------------------------------------------------------------------------------------------
## Plot 2D representation of network cell positions and connections
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplot2Dnet(include=['allCells'], view='xy', showConns=True, popColors=None, tagsFile=None, figSize=(12,12), fontSize=12, saveData=None, saveFig=None, showFig=True, **kwargs):
"""Plots 2D representation of network cell positions and connections.
Parameters
----------
include : list
List of presynaptic cells to include.
**Default:** ``['allCells']``
**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
view : str
Perspective of view.
**Default:** ``'xy'`` front view,
**Options:** ``'xz'`` top-down view
showConns : bool
Whether to show connections or not.
**Default:** ``True``
popColors : dict
Dictionary with custom color (value) used for each population (key).
**Default:** ``None`` uses standard colors
tagsFile : str
Path to a saved tags file to use in connectivity plot.
**Default:** ``None``
figSize : list [width, height]
Size of figure in inches.
**Default:** ``(12, 12)``
fontSize : int
Font size on figure.
**Default:** ``12``
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.ext'`` saves to a custom path and filename, valid file extensions are ``'.png'``, ``'.jpg'``, ``'.eps'``, and ``'.tiff'``
showFig : bool
Shows the figure if ``True``.
**Default:** ``True``
Returns
-------
(fig, dict)
A tuple consisting of the matplotlib figure handle and a dictionary containing the plot data.
See Also
--------
iplot2Dnet :
Examples
--------
>>> import netpyne, netpyne.examples.example
>>> out = netpyne.analysis.plot2Dnet()
"""
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 2D representation of network cell locations and connections...')
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
TOOLS = 'hover,save,pan,box_zoom,reset,wheel_zoom',
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList] # bokeh only handles integer rgb values from 0-255
else:
colors = kwargs['palette']
# front view
if view == 'xy':
ycoord = 'y'
elif view == 'xz':
ycoord = 'z'
if tagsFile:
print('Loading tags file...')
import json
with open(tagsFile, 'r') as fileObj: tagsTmp = json.load(fileObj)['tags']
tagsFormat = tagsTmp.pop('format', [])
tags = {int(k): v for k,v in tagsTmp.items()} # find method to load json with int keys?
del tagsTmp
# set indices of fields to read compact format (no keys)
missing = []
popIndex = tagsFormat.index('pop') if 'pop' in tagsFormat else missing.append('pop')
xIndex = tagsFormat.index('x') if 'x' in tagsFormat else missing.append('x')
yIndex = tagsFormat.index('y') if 'y' in tagsFormat else missing.append('y')
zIndex = tagsFormat.index('z') if 'z' in tagsFormat else missing.append('z')
if len(missing) > 0:
print("Missing:")
print(missing)
return None, None, None
# find pre and post cells
if tags:
cellGids = getCellsIncludeTags(include, tags, tagsFormat)
popLabels = list(set([tags[gid][popIndex] for gid in cellGids]))
# pop and cell colors
popColorsTmp = {popLabel: colorList[ipop%len(colorList)] for ipop,popLabel in enumerate(popLabels)} # dict with color for each pop
if popColors: popColorsTmp.update(popColors)
popColors = popColorsTmp
cellColors = [popColors[tags[gid][popIndex]] for gid in cellGids]
# cell locations
posX = [tags[gid][xIndex] for gid in cellGids] # get all x positions
if ycoord == 'y':
posY = [tags[gid][yIndex] for gid in cellGids] # get all y positions
elif ycoord == 'z':
posY = [tags[gid][zIndex] for gid in cellGids] # get all y positions
else:
print('Error loading tags from file')
return None
else:
cells, cellGids, _ = getCellsInclude(include)
selectedPops = [cell['tags']['pop'] for cell in cells]
popLabels = [pop for pop in sim.net.allPops if pop in selectedPops] # preserves original ordering
# pop and cell colors
popColorsTmp = {popLabel: colorList[ipop%len(colorList)] for ipop,popLabel in enumerate(popLabels)} # dict with color for each pop
if popColors: popColorsTmp.update(popColors)
popColors = popColorsTmp
cellColors = ['#%02x%02x%02x' % tuple([int(col * 255) for col in popColors[cell['tags']['pop']]]) for cell in cells]
# cell locations
posX = [cell['tags']['x'] for cell in cells] # get all x positions
posY = [cell['tags'][ycoord] for cell in cells] # get all y positions
fig = figure(
title="2D Network representation",
tools=TOOLS,
active_drag = None,
active_scroll = None,
tooltips=[('y location', '@y'), ('x location', '@x')],
x_axis_label="x (um)",
y_axis_label='y (um)',
x_range=[min(posX)-0.05*max(posX),1.05*max(posX)],
y_range=[1.05*max(posY), min(posY)-0.05*max(posY)],
toolbar_location='above')
if 'radius' in kwargs:
radius = kwargs['radius']
else:
radius = 1.0
fig.scatter(posX, posY, radius=radius, fill_color=cellColors, line_color=None) # plot cell soma positions
posXpre, posYpre = [], []
posXpost, posYpost = [], []
if showConns and not tagsFile:
for postCell in cells:
for con in postCell['conns']: # plot connections between cells
if not isinstance(con['preGid'], basestring) and con['preGid'] in cellGids:
posXpre,posYpre = next(((cell['tags']['x'],cell['tags'][ycoord]) for cell in cells if cell['gid']==con['preGid']), None)
posXpost, posYpost = postCell['tags']['x'], postCell['tags'][ycoord]
if kwargs.get('theme', '') == 'gui':
color = 'yellow'
else:
color='red'
if con['synMech'] in ['inh', 'GABA', 'GABAA', 'GABAB']:
if kwargs.get('theme', '') == 'gui':
color = 'lightcyan'
else:
color = 'blue'
width = 0.1 #50*con['weight']
fig.line([posXpre, posXpost], [posYpre, posYpost], color=color, line_width=width) # plot line from pre to post
fontsiz = fontSize
for popLabel in popLabels:
fig.line(0,0,color=tuple([int(x*255) for x in popColors[popLabel]]), legend_label=popLabel)
legend = Legend(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 = {'posX': posX, 'posY': posY, 'posX': cellColors, 'posXpre': posXpre, 'posXpost': posXpost, 'posYpre': posYpre, 'posYpost': posYpost,
'include': include, 'saveData': saveData, 'saveFig': saveFig, 'showFig': showFig}
_saveFigData(figData, saveData, '2Dnet')
# save figure
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename+'_iplot_2Dnet_'+groupBy+'_'+feature+'_'+graphType+'.html'
file = open(filename, 'w')
file.write(html)
file.close()
# show fig
if showFig: show(fig)
return html, {'include': include, 'posX': posX, 'posY': posY, 'posXpre': posXpre, 'posXpost': posXpost, 'posYpre': posYpre, 'posYpost': posYpost}
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive RxD concentration
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotRxDConcentration(speciesLabel, regionLabel, plane='xy', saveFig=None, showFig=True, **kwargs):
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
from bokeh.transform import linear_cmap
from bokeh.models import ColorBar
print('Plotting interactive RxD concentration ...')
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
if not 'palette' in kwargs:
from bokeh.palettes import Viridis256
colors = Viridis256
#colors = [RGB(*[round(f * 255) for f in color]) for color in colorList]
else:
colors = kwargs['palette']
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
species = sim.net.rxd['species'][speciesLabel]['hObj']
region = sim.net.rxd['regions'][regionLabel]['hObj']
plane2mean = {'xz': 1, 'xy': 2}
data = species[region].states3d[:].mean(plane2mean[plane]).T
data = np.flipud(data) # This does not flip the axis values, should improve
dh, dw = np.shape(data)
low = np.nanmin(data)
high = np.nanmax(data)
if low == high:
low = low - 0.1 * low
high = high + 0.1 * high
conc_colormapper = linear_cmap(
field_name = '[' + species.name + '] (mM)',
palette = colors,
low = low,
high = high,
nan_color = 'white',
)
conc_colorbar = ColorBar(
color_mapper = conc_colormapper['transform'],
label_standoff = 12,
title_standoff = 12,
)
conc_colorbar.title = '[' + species.name + '] (mM)'
fig = figure(
title = 'RxD: ' + species.name + ' concentration',
toolbar_location = 'above',
tools = 'hover,save,pan,box_zoom,reset,wheel_zoom',
active_drag = None,
active_scroll = None,
tooltips = [("x", "$x"), ("y", "$y"), ("value", "@image")],
match_aspect = True,
x_axis_label = plane[0] + " location (um)",
y_axis_label = plane[1] + " location (um)",
)
fig.image(image=[data], x=0, y=0, dw=dw, dh=dh, color_mapper=conc_colormapper['transform'], level="image")
fig.add_layout(conc_colorbar, 'right')
plot_layout = layout([fig], sizing_mode='scale_height')
html = file_html(plot_layout, CDN, title="RxD Concentration")
if showFig:
show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename + '_RxD_concentration.html'
outfile = open(filename, 'w')
outfile.write(html)
outfile.close()
return html, data
# -------------------------------------------------------------------------------------------------------------------
## Plot interactive spike statistics
# -------------------------------------------------------------------------------------------------------------------
[docs]@exception
def iplotSpikeStats(include=['eachPop', 'allCells'], statDataIn={}, timeRange=None, graphType='boxplot', stats=['rate', 'isicv'], bins=50, histlogy=False, histlogx=False, histmin=0.0, density=False, includeRate0=False, legendLabels=None, normfit=False, histShading=True, xlim=None, popColors={}, saveData=None, saveFig=None, showFig=True, **kwargs):
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
from bokeh.transform import linear_cmap
from bokeh.transform import factor_cmap
from bokeh.models import ColorBar
from bokeh.palettes import Spectral6
from bokeh.models.mappers import CategoricalColorMapper
print('Plotting interactive spike statistics ...')
TOOLS = "pan,wheel_zoom,box_zoom,reset,save,box_select"
xlabels = {'rate': 'Rate (Hz)', 'isicv': 'Irregularity (ISI CV)', 'sync': 'Synchrony', 'pairsync': 'Pairwise synchrony'}
# 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]
for stat in stats:
if 'theme' in kwargs:
if kwargs['theme'] != 'default':
if kwargs['theme'] == 'gui':
from bokeh.themes import Theme
theme = Theme(json=_guiTheme)
else:
theme = kwargs['theme']
curdoc().theme = theme
else:
curdoc().theme = bokeh_theme
else:
curdoc().theme = bokeh_theme
if not 'palette' in kwargs:
colors = [RGB(*[round(f * 255) for f in color]) for color in colorList]
else:
colors = kwargs['palette']
xlabel = xlabels[stat]
statData = []
gidsData = []
ynormsData = []
# Calculate data for each entry in include
for iplot,subset in enumerate(include):
if stat in statDataIn:
statData = statDataIn[stat]['statData']
gidsData = statDataIn[stat].get('gidsData', [])
ynormsData = statDataIn[stat].get('ynormsData', [])
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 = [],[]
# 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
try:
spkts,spkinds = list(zip(*[(spkt, spkind) for spkt, spkind in zip(spkts, spkinds)
if timeRange[0] <= spkt <= timeRange[1]]))
except:
pass
# rate stats
if stat == 'rate':
toRate = 1e3/(timeRange[1]-timeRange[0])
if includeRate0:
rates = [spkinds.count(gid)*toRate for gid in cellGids] \
if len(spkinds)>0 else [0]*len(cellGids) #cellGids] #set(spkinds)]
else:
rates = [spkinds.count(gid)*toRate for gid in set(spkinds)] \
if len(spkinds)>0 else [0] #cellGids] #set(spkinds)]
statData.append(rates)
# Inter-spike interval (ISI) coefficient of variation (CV) stats
elif stat == 'isicv':
import numpy as np
spkmat = [[spkt for spkind,spkt in zip(spkinds,spkts) if spkind==gid]
for gid in set(spkinds)]
isimat = [[t - s for s, t in zip(spks, spks[1:])] for spks in spkmat if len(spks)>10]
isicv = [np.std(x) / np.mean(x) if len(x)>0 else 0 for x in isimat] # if len(x)>0]
statData.append(isicv)
# synchrony
elif stat in ['sync', 'pairsync']:
try:
import pyspike
except:
print("Error: plotSpikeStats() requires the PySpike python package \
to calculate synchrony (try: pip install pyspike)")
return 0
spkmat = [pyspike.SpikeTrain([spkt for spkind,spkt in zip(spkinds,spkts)
if spkind==gid], timeRange) for gid in set(spkinds)]
if stat == 'sync':
# (SPIKE-Sync measure)' # see http://www.scholarpedia.org/article/Measures_of_spike_train_synchrony
syncMat = [pyspike.spike_sync(spkmat)]
#graphType = 'bar'
elif stat == 'pairsync':
# (SPIKE-Sync measure)' # see http://www.scholarpedia.org/article/Measures_of_spike_train_synchrony
syncMat = np.mean(pyspike.spike_sync_matrix(spkmat), 0)
statData.append(syncMat)
# boxplot
if graphType == 'boxplot':
line_width = 2
line_color = 'black'
if 'theme' in kwargs:
if kwargs['theme'] == 'gui' or 'dark' in kwargs['theme']:
line_color = 'lightgray'
labels = legendLabels if legendLabels else include
data_lists = statData[::-1]
data_lists.reverse()
box_colors = colors[0:len(labels)]
if include[0] == 'allCells':
del box_colors[-1]
box_colors.insert(0, ('darkslategray'))
data = {}
for label, data_list in zip(labels, data_lists):
data[label] = data_list
fig = figure(
title = 'Spike statistics: ' + xlabels[stat],
toolbar_location = 'above',
tools = 'hover,save,pan,box_zoom,reset,wheel_zoom',
active_drag = None,
active_scroll = None,
tooltips = [(stat, "$y")],
x_axis_label = 'Population',
y_axis_label = xlabel,
x_range = labels
)
fig.xgrid.visible = False
df = pd.DataFrame(dict([ (k,pd.Series(v)) for k,v in data.items() ]))
q1 = df.quantile(q=0.25)
q2 = df.quantile(q=0.5)
q3 = df.quantile(q=0.75)
iqr = q3 - q1
upper = q3 + 1.5 * iqr
lower = q1 - 1.5 * iqr
qmin = df.quantile(q=0.00)
qmax = df.quantile(q=1.00)
qmean = df.mean()
out_highs = df[df > upper]
out_lows = df[df < lower]
# outliers
for index, label in enumerate(labels):
out_x = []
out_y = []
out_high = out_highs[label].dropna()
if not out_high.empty:
upper[label] = df[df < upper][label].max()
for val in out_high:
out_x.append(label)
out_y.append(val)
else:
upper[label] = qmax[label]
out_low = out_lows[label].dropna()
if not out_low.empty:
lower[label] = df[df > lower][label].min()
for val in out_low:
out_x.append(label)
out_y.append(val)
else:
lower[label] = qmin[label]
if out_x:
fig.circle_x(out_x, out_y, size=10, fill_color=box_colors[index], fill_alpha=0.8, line_color=line_color)
# stems
fig.segment(labels, upper, labels, q3, line_color=line_color, line_width=line_width)
fig.segment(labels, lower, labels, q1, line_color=line_color, line_width=line_width)
# boxes
fig.vbar(labels, 0.7, q2, q3, line_color=line_color, line_width=line_width, fill_color=box_colors)
fig.vbar(labels, 0.7, q1, q2, line_color=line_color, line_width=line_width, fill_color=box_colors)
# whiskers (almost-0 height rects simpler than segments)
fig.rect(labels, lower, 20, 1, width_units='screen', height_units='screen', line_color=line_color, line_width=line_width)
fig.rect(labels, upper, 20, 1, width_units='screen', height_units='screen', line_color=line_color, line_width=line_width)
# means
fig.circle_cross(labels, qmean, size=10, fill_color='white', fill_alpha=0.5, line_color='black')
else:
raise Exception('Only boxplot is currently supported in iplotSpikeStats.')
plot_layout = layout([fig], sizing_mode='stretch_both')
html = file_html(plot_layout, CDN, title="Spike Statistics")
if showFig:
show(plot_layout)
if saveFig:
if isinstance(saveFig, str):
filename = saveFig
else:
filename = sim.cfg.filename + '_spike_' + stat + '.html'
outfile = open(filename, 'w')
outfile.write(html)
outfile.close()
return html