"""
Module for analyzing and plotting LFP-related results
"""
try:
basestring
except NameError:
basestring = str
from netpyne import __gui__
if __gui__:
import matplotlib.pyplot as plt
from matplotlib import mlab
import numpy as np
from numbers import Number
from .utils import colorList, exception, _saveFigData, _showFigure, _smooth1d
# -------------------------------------------------------------------------------------------------------------------
## Plot LFP (time-resolved, power spectral density, time-frequency and 3D locations)
# -------------------------------------------------------------------------------------------------------------------
[docs]
@exception
def plotLFP(
timeRange=None,
electrodes=['avg', 'all'],
plots=['timeSeries', 'PSD', 'spectrogram', 'locations'],
inputLFP=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,
filtFreq=False,
filtOrder=3,
detrend=False,
transformMethod='morlet',
maxPlots=8,
overlay=False,
colors=None,
figSize=(8, 8),
fontSize=14,
lineWidth=1.5,
dpi=200,
saveData=None,
saveFig=None,
showFig=True,
):
"""
Function for plotting local field potentials (LFP)
Parameters
----------
timeRange : list [start, stop]
Time range to plot.
**Default:**
``None`` plots entire time range
electrodes : list
List of electrodes to include; ``'avg'`` is the average of all electrodes; ``'all'`` is each electrode separately.
**Default:** ``['avg', 'all']``
plots : list
List of plot types to show.
**Default:** ``['timeSeries', 'PSD', 'spectrogram', 'locations']``
NFFT : int (power of 2)
Number of data points used in each block for the PSD and time-freq FFT.
**Default:** ``256``
noverlap : int (<nperseg)
Number of points of overlap between segments for PSD and time-freq.
**Default:** ``128``
nperseg : int
Length of each segment for time-freq.
**Default:** ``256``
minFreq : float
Minimum frequency shown in plot for PSD and time-freq.
**Default:** ``1``
maxFreq : float
Maximum frequency shown in plot for PSD and time-freq.
**Default:** ``100``
stepFreq : float
Step frequency.
**Default:** ``1``
smooth : int
Window size for smoothing LFP; no smoothing if ``0``
**Default:** ``0``
separation : float
Separation factor between time-resolved LFP plots; multiplied by max LFP value.
**Default:** ``1.0``
includeAxon : bool
Whether to show the axon in the location plot.
**Default:** ``True``
logx : bool
Whether to make x-axis logarithmic
**Default:** ``False``
**Options:** ``<option>`` <description of option>
logy : bool
Whether to make y-axis logarithmic
**Default:** ``False``
normSignal : bool
Whether to normalize the signal.
**Default:** ``False``
normPSD : bool
Whether to normalize the power spectral density.
**Default:** ``False``
normSpec : bool
Needs documentation.
**Default:** ``False``
filtFreq : int or list
Frequency for low-pass filter (int) or frequencies for bandpass filter in a list: [low, high]
**Default:** ``False`` does not filter the data
filtOrder : int
Order of the filter defined by `filtFreq`.
**Default:** ``3``
detrend : bool
Whether to detrend.
**Default:** ``False``
transformMethod : str
Transform method.
**Default:** ``'morlet'``
**Options:** ``'fft'``
maxPlots : int
Maximum number of subplots. Currently unused.
**Default:** ``8``
overlay : bool
Whether to overlay plots or use subplots.
**Default:** ``False`` overlays plots.
colors : list
List of normalized RGB colors to use.
**Default:** ``None`` uses standard colors
figSize : list [width, height]
Size of figure in inches.
**Default:** ``(10, 8)``
fontSize : int
Font size on figure.
**Default:** ``14``
lineWidth : int
Line width.
**Default:** ``1.5``
dpi : int
Resolution of figure in dots per inch.
**Default:** ``100``
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
-------
(figs, dict)
A tuple consisting of the Matplotlib figure handles and a dictionary containing the plot data
See Also
--------
iplotLFP :
"""
# Note: should probably split funcs for signal, psd, spectrogram and locs
from .. import sim
from ..support.scalebar import add_scalebar
print('Plotting LFP ...')
if not colors:
colors = colorList
# set font size
plt.rcParams.update({'font.size': fontSize})
# time range
if timeRange is None:
timeRange = [0, sim.cfg.duration]
if inputLFP is not None:
lfp = inputLFP[int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), :]
else:
lfp = np.array(sim.allSimData['LFP'])[
int(timeRange[0] / sim.cfg.recordStep) : int(timeRange[1] / sim.cfg.recordStep), :
]
if filtFreq:
from scipy import signal
fs = 1000.0 / sim.cfg.recordStep
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)
for i in range(lfp.shape[1]):
lfp[:, i] = signal.filtfilt(b, a, lfp[:, i])
if detrend:
from scipy import signal
for i in range(lfp.shape[1]):
lfp[:, i] = signal.detrend(lfp[:, i])
if normSignal:
for i in range(lfp.shape[1]):
offset = min(lfp[:, i])
if offset <= 0:
lfp[:, i] += abs(offset)
lfp[:, i] /= max(lfp[:, i])
# electrode selection
if 'all' in electrodes:
electrodes.remove('all')
electrodes.extend(list(range(int(sim.net.recXElectrode.nsites))))
# plotting
figs = []
# maxPlots = 8.0
data = {'lfp': lfp} # returned data
# time series -----------------------------------------
if 'timeSeries' in plots:
ydisp = np.absolute(lfp).max() * separation
offset = 1.0 * ydisp
t = np.arange(timeRange[0], timeRange[1], sim.cfg.recordStep)
if figSize:
figs.append(plt.figure(figsize=figSize))
for i, elec in enumerate(electrodes):
if elec == 'avg':
lfpPlot = np.mean(lfp, axis=1)
color = 'k'
lw = 1.0
elif isinstance(elec, Number) and (inputLFP is not None or elec <= sim.net.recXElectrode.nsites):
lfpPlot = lfp[:, elec]
color = colors[i % len(colors)]
lw = 1.0
if len(t) < len(lfpPlot):
lfpPlot = lfpPlot[: len(t)]
plt.plot(t[0 : len(lfpPlot)], -lfpPlot + (i * ydisp), color=color, linewidth=lw)
if len(electrodes) > 1:
plt.text(
timeRange[0] - 0.07 * (timeRange[1] - timeRange[0]),
(i * ydisp),
elec,
color=color,
ha='center',
va='top',
fontsize=fontSize,
fontweight='bold',
)
ax = plt.gca()
data['lfpPlot'] = lfpPlot
data['ydisp'] = ydisp
data['t'] = t
# format plot
if len(electrodes) > 1:
plt.text(
timeRange[0] - 0.14 * (timeRange[1] - timeRange[0]),
(len(electrodes) * ydisp) / 2.0,
'LFP electrode',
color='k',
ha='left',
va='bottom',
fontsize=fontSize,
rotation=90,
)
plt.ylim(-offset, (len(electrodes)) * ydisp)
else:
plt.suptitle('LFP Signal', fontSize=fontSize, fontweight='bold')
ax.invert_yaxis()
plt.xlabel('time (ms)', fontsize=fontSize)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
plt.subplots_adjust(bottom=0.1, top=1.0, right=1.0)
# calculate scalebar size and add scalebar
round_to_n = lambda x, n, m: int(np.ceil(round(x, -int(np.floor(np.log10(abs(x)))) + (n - 1)) / m)) * m
scaley = 1000.0 # values in mV but want to convert to uV
m = 10.0
sizey = 100 / scaley
while sizey > 0.25 * ydisp:
try:
sizey = round_to_n(0.2 * ydisp * scaley, 1, m) / scaley
except:
sizey /= 10.0
m /= 10.0
labely = '%.3g $\mu$V' % (sizey * scaley) # )[1:]
if len(electrodes) > 1:
add_scalebar(
ax,
hidey=True,
matchy=False,
hidex=False,
matchx=False,
sizex=0,
sizey=-sizey,
labely=labely,
unitsy='$\mu$V',
scaley=scaley,
loc=3,
pad=0.5,
borderpad=0.5,
sep=3,
prop=None,
barcolor="black",
barwidth=2,
)
else:
add_scalebar(
ax,
hidey=True,
matchy=False,
hidex=True,
matchx=True,
sizex=None,
sizey=-sizey,
labely=labely,
unitsy='$\mu$V',
scaley=scaley,
unitsx='ms',
loc=3,
pad=0.5,
borderpad=0.5,
sep=3,
prop=None,
barcolor="black",
barwidth=2,
)
# save figure
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename + '_LFP_timeseries.png'
plt.savefig(filename, dpi=dpi)
# PSD ----------------------------------
if 'PSD' in plots:
if overlay:
figs.append(plt.figure(figsize=figSize))
else:
numCols = 1 # np.round(len(electrodes) / maxPlots) + 1
figs.append(plt.figure(figsize=(figSize[0] * numCols, figSize[1])))
# import seaborn as sb
allFreqs = []
allSignal = []
data['allFreqs'] = allFreqs
data['allSignal'] = allSignal
for i, elec in enumerate(electrodes):
if elec == 'avg':
lfpPlot = np.mean(lfp, axis=1)
elif isinstance(elec, Number) and (inputLFP is not None or elec <= sim.net.recXElectrode.nsites):
lfpPlot = lfp[:, elec]
# Morlet wavelet transform method
if transformMethod == 'morlet':
from ..support.morlet import MorletSpec, index2ms
Fs = int(1000.0 / sim.cfg.recordStep)
# t_spec = np.linspace(0, index2ms(len(lfpPlot), Fs), len(lfpPlot))
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)
# ALTERNATIVE PSD CALCULATION USING WELCH
# from http://joelyancey.com/lfp-python-practice/
# from scipy import signal as spsig
# Fs = int(1000.0/sim.cfg.recordStep)
# maxFreq=100
# f, psd = spsig.welch(lfpPlot, Fs, nperseg=100)
# plt.semilogy(f,psd,'k')
# sb.despine()
# plt.xlim((0,maxFreq))
# plt.yticks(size=fontsiz)
# plt.xticks(size=fontsiz)
# plt.ylabel('$uV^{2}/Hz$',size=fontsiz)
if normPSD:
vmax = np.max(allSignal)
for i, s in enumerate(allSignal):
allSignal[i] = allSignal[i] / vmax
for i, elec in enumerate(electrodes):
if not overlay:
plt.subplot(int(np.ceil(len(electrodes) / numCols)), numCols, i + 1)
if elec == 'avg':
color = 'k'
elif isinstance(elec, Number) and (inputLFP is not None or elec <= sim.net.recXElectrode.nsites):
color = colors[i % len(colors)]
freqs = allFreqs[i]
signal = allSignal[i]
plt.plot(
freqs[freqs < maxFreq],
signal[freqs < maxFreq],
linewidth=lineWidth,
color=color,
label='Electrode %s' % (str(elec)),
)
plt.xlim([0, maxFreq])
if len(electrodes) > 1 and not overlay:
plt.title('Electrode %s' % (str(elec)), fontsize=fontSize)
plt.ylabel(ylabel, fontsize=fontSize)
# format plot
plt.xlabel('Frequency (Hz)', fontsize=fontSize)
if overlay:
plt.legend(fontsize=fontSize)
plt.tight_layout()
plt.suptitle('LFP Power Spectral Density', fontsize=fontSize, fontweight='bold') # add yaxis in opposite side
plt.subplots_adjust(bottom=0.08, top=0.92)
if logx:
pass
# from IPython import embed; embed()
# save figure
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename + '_LFP_psd.png'
plt.savefig(filename, dpi=dpi)
# Spectrogram ------------------------------
if 'spectrogram' in plots:
import matplotlib.cm as cm
numCols = 1 # np.round(len(electrodes) / maxPlots) + 1
figs.append(plt.figure(figsize=(figSize[0] * numCols, figSize[1])))
# Morlet wavelet transform method
if transformMethod == 'morlet':
from ..support.morlet import MorletSpec, index2ms
spec = []
freqList = None
if logy:
freqList = np.logspace(np.log10(minFreq), np.log10(maxFreq), int((maxFreq - minFreq) / stepFreq))
for i, elec in enumerate(electrodes):
if elec == 'avg':
lfpPlot = np.mean(lfp, axis=1)
elif isinstance(elec, Number) and (inputLFP is not None or 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, lfreq=freqList)
)
f = (
freqList if freqList is not None else np.array(range(minFreq, maxFreq + 1, stepFreq))
) # only used as output for user
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):
plt.subplot(int(np.ceil(len(electrodes) / numCols)), numCols, i + 1)
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]
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('viridis'),
)
if normSpec:
plt.colorbar(label='Normalized power')
else:
plt.colorbar(label='Power')
plt.ylabel('Hz')
plt.tight_layout()
if len(electrodes) > 1:
plt.title('Electrode %s' % (str(elec)), fontsize=fontSize - 2)
# FFT transform method
elif transformMethod == 'fft':
from scipy import signal as spsig
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])
spec.append(10 * np.log10(x_spec[f < maxFreq]))
vmin = np.array(spec).min()
vmax = np.array(spec).max()
for i, elec in enumerate(electrodes):
plt.subplot(np.ceil(len(electrodes) / numCols), numCols, i + 1)
plt.pcolormesh(x_mesh, y_mesh, spec[i], cmap=cm.viridis, vmin=vmin, vmax=vmax)
plt.colorbar(label='dB/Hz', ticks=[np.ceil(vmin), np.floor(vmax)])
if logy:
plt.yscale('log')
plt.ylabel('Log-frequency (Hz)')
if isinstance(logy, list):
yticks = tuple(logy)
plt.yticks(yticks, yticks)
else:
plt.ylabel('(Hz)')
if len(electrodes) > 1:
plt.title('Electrode %s' % (str(elec)), fontsize=fontSize - 2)
plt.xlabel('time (ms)', fontsize=fontSize)
plt.tight_layout()
plt.suptitle('LFP spectrogram', size=fontSize, fontweight='bold')
plt.subplots_adjust(bottom=0.08, top=0.90)
# save figure
if saveFig:
if isinstance(saveFig, basestring):
filename = saveFig
else:
filename = sim.cfg.filename + '_LFP_timefreq.png'
plt.savefig(filename, dpi=dpi)
# locations ------------------------------
if 'locations' in plots:
try:
cvals = [] # used to store total transfer resistance
for cell in sim.net.compartCells:
trSegs = list(
np.sum(sim.net.recXElectrode.getTransferResistance(cell.gid) * 1e3, axis=0)
) # convert from Mohm to kilohm
if not includeAxon:
i = 0
for secName, sec in cell.secs.items():
nseg = sec['hObj'].nseg # .geom.nseg
if 'axon' in secName:
for j in range(i, i + nseg):
del trSegs[j]
i += nseg
cvals.extend(trSegs)
includePost = [c.gid for c in sim.net.compartCells]
fig = sim.analysis.plotShape(
includePost=includePost,
showElectrodes=electrodes,
cvals=cvals,
includeAxon=includeAxon,
dpi=dpi,
fontSize=fontSize,
saveFig=saveFig,
showFig=showFig,
figSize=figSize,
)[0]
figs.append(fig)
except:
print(' Failed to plot LFP locations...')
outputData = {
'LFP': lfp,
'electrodes': electrodes,
'timeRange': timeRange,
'saveData': saveData,
'saveFig': saveFig,
'showFig': showFig,
}
if 'timeSeries' in plots:
outputData.update({'t': t})
if 'PSD' in plots:
outputData.update({'allFreqs': allFreqs, 'allSignal': allSignal})
if 'spectrogram' in plots:
outputData.update({'spec': spec, 't': t_spec * 1000.0, 'freqs': f[f <= maxFreq]})
# save figure data
if saveData:
figData = outputData
_saveFigData(figData, saveData, 'lfp')
# show fig
if showFig:
_showFigure()
return figs, outputData