"""Class to handle the dipoles."""
# Authors: Mainak Jas <mjas@mgh.harvard.edu>
# Sam Neymotin <samnemo@gmail.com>
import warnings
import numpy as np
from copy import deepcopy
from .externals.mne import _check_option
from .viz import plot_dipole, plot_psd, plot_tfr_morlet
[docs]def simulate_dipole(net, tstop, dt=0.025, n_trials=None, record_vsec=False,
record_isec=False, postproc=False):
"""Simulate a dipole given the experiment parameters.
Parameters
----------
net : Network object
The Network object specifying how cells are
connected.
tstop : float
The simulation stop time (ms).
dt : float
The integration time step of h.CVode (ms)
n_trials : int | None
The number of trials to simulate. If None, the 'N_trials' value
of the ``params`` used to create ``net`` is used (must be >0)
record_vsec : 'all' | 'soma' | False
Option to record voltages from all sections ('all'), or just
the soma ('soma'). Default: False.
record_isec : 'all' | 'soma' | False
Option to record voltages from all sections ('all'), or just
the soma ('soma'). Default: False.
postproc : bool
If True, smoothing (``dipole_smooth_win``) and scaling
(``dipole_scalefctr``) values are read from the parameter file, and
applied to the dipole objects before returning. Note that this setting
only affects the dipole waveforms, and not somatic voltages, possible
extracellular recordings etc. The preferred way is to use the
:meth:`~hnn_core.dipole.Dipole.smooth` and
:meth:`~hnn_core.dipole.Dipole.scale` methods instead. Default: False.
Returns
-------
dpls: list
List of dipole objects for each trials
"""
from .parallel_backends import _BACKEND, JoblibBackend
if _BACKEND is None:
_BACKEND = JoblibBackend(n_jobs=1)
if n_trials is None:
n_trials = net._params['N_trials']
if n_trials < 1:
raise ValueError("Invalid number of simulations: %d" % n_trials)
if not net.connectivity:
warnings.warn('No connections instantiated in network. Consider using '
'net = jones_2009_model() or net = law_2021_model() to '
'create a predefined network from published models.',
UserWarning)
# ADD DRIVE WARNINGS HERE
if not net.external_drives and not net.external_biases:
warnings.warn('No external drives or biases loaded', UserWarning)
for drive_name, drive in net.external_drives.items():
if 'tstop' in drive['dynamics']:
if drive['dynamics']['tstop'] is None:
drive['dynamics']['tstop'] = tstop
for bias_name, bias in net.external_biases.items():
for cell_type, bias_cell_type in bias.items():
if bias_cell_type['tstop'] is None:
bias_cell_type['tstop'] = tstop
if bias_cell_type['tstop'] < 0.:
raise ValueError('End time of tonic input cannot be negative')
duration = bias_cell_type['tstop'] - bias_cell_type['t0']
if duration < 0.:
raise ValueError('Duration of tonic input cannot be negative')
net._instantiate_drives(n_trials=n_trials, tstop=tstop)
net._reset_rec_arrays()
_check_option('record_vsec', record_vsec, ['all', 'soma', False])
net._params['record_vsec'] = record_vsec
_check_option('record_isec', record_isec, ['all', 'soma', False])
net._params['record_isec'] = record_isec
if postproc:
warnings.warn('The postproc-argument is deprecated and will be removed'
' in a future release of hnn-core. Please define '
'smoothing and scaling explicitly using Dipole methods.',
DeprecationWarning)
dpls = _BACKEND.simulate(net, tstop, dt, n_trials, postproc)
return dpls
[docs]def read_dipole(fname):
"""Read dipole values from a file and create a Dipole instance.
Parameters
----------
fname : str
Full path to the input file (.txt)
Returns
-------
dpl : Dipole
The instance of Dipole class
"""
dpl_data = np.loadtxt(fname, dtype=float)
ncols = dpl_data.shape[1]
if ncols not in (2, 4):
raise ValueError(
f'Data are supposed to have 2 or 4 columns while we have {ncols}.')
dpl = Dipole(dpl_data[:, 0], dpl_data[:, 1:])
return dpl
[docs]def average_dipoles(dpls):
"""Compute dipole averages over a list of Dipole objects.
Parameters
----------
dpls : list of Dipole objects
Contains list of dipole objects, each with a `data` member containing
'L2', 'L5' and 'agg' components
Returns
-------
dpl : instance of Dipole
A new dipole object with each component of `dpl.data` representing the
average over the same components in the input list
"""
scale_applied = dpls[0].scale_applied
for dpl_idx, dpl in enumerate(dpls):
if dpl.scale_applied != scale_applied:
raise RuntimeError('All dipoles must be scaled equally!')
if not isinstance(dpl, Dipole):
raise ValueError(
f"All elements in the list should be instances of "
f"Dipole. Got {type(dpl)}")
if dpl.nave > 1:
raise ValueError("Dipole at index %d was already an average of %d"
" trials. Cannot reaverage" %
(dpl_idx, dpl.nave))
avg_data = list()
layers = dpl.data.keys()
for layer in layers:
avg_data.append(
np.mean(np.array([dpl.data[layer] for dpl in dpls]), axis=0)
)
avg_data = np.c_[avg_data].T
avg_dpl = Dipole(dpls[0].times, avg_data)
# The averaged scale should equal all scals in the input dpl list.
avg_dpl.scale_applied = scale_applied
# set nave to the number of trials averaged in this dipole
avg_dpl.nave = len(dpls)
return avg_dpl
def _rmse(dpl, exp_dpl, tstart=0.0, tstop=0.0, weights=None):
""" Calculates RMSE between data in dpl and exp_dpl
Parameters
----------
dpl : instance of Dipole
A dipole object with simulated data
exp_dpl : instance of Dipole
A dipole object with experimental data
tstart : None | float
Time at beginning of range over which to calculate RMSE
tstop : None | float
Time at end of range over which to calculate RMSE
weights : None | array
An array of weights to be applied to each point in
simulated dpl. Must have length >= dpl.data
If None, weights will be replaced with 1's for typical RMSE
calculation.
Returns
-------
err : float
Weighted RMSE between data in dpl and exp_dpl
"""
from scipy import signal
exp_times = exp_dpl.times
sim_times = dpl.times
# do tstart and tstop fall within both datasets?
# if not, use the closest data point as the new tstop/tstart
for tseries in [exp_times, sim_times]:
if tstart < tseries[0]:
tstart = tseries[0]
if tstop > tseries[-1]:
tstop = tseries[-1]
# make sure start and end times are valid for both dipoles
exp_start_index = (np.abs(exp_times - tstart)).argmin()
exp_end_index = (np.abs(exp_times - tstop)).argmin()
exp_length = exp_end_index - exp_start_index
sim_start_index = (np.abs(sim_times - tstart)).argmin()
sim_end_index = (np.abs(sim_times - tstop)).argmin()
sim_length = sim_end_index - sim_start_index
if weights is None:
# weighted RMSE with weights of all 1's is equivalent to
# normal RMSE
weights = np.ones(len(sim_times[0:sim_end_index]))
weights = weights[sim_start_index:sim_end_index]
dpl1 = dpl.data['agg'][sim_start_index:sim_end_index]
dpl2 = exp_dpl.data['agg'][exp_start_index:exp_end_index]
if (sim_length > exp_length):
# downsample simulation timeseries to match exp data
dpl1 = signal.resample(dpl1, exp_length)
weights = signal.resample(weights, exp_length)
indices = np.where(weights < 1e-4)
weights[indices] = 0
elif (sim_length < exp_length):
# downsample exp timeseries to match simulation data
dpl2 = signal.resample(dpl2, sim_length)
return np.sqrt((weights * ((dpl1 - dpl2) ** 2)).sum() / weights.sum())
[docs]class Dipole(object):
"""Dipole class.
An instance of the ``Dipole``-class contains the simulated dipole moment
timecourses for L2 and L5 pyramidal cells, as well as their aggregate
(``'agg'``). The units of the dipole moment are in ``nAm``
(1e-9 Ampere-meters).
Parameters
----------
times : array (n_times,)
The time vector (in ms)
data : array, shape (n_times x n_layers)
The data. The first column represents 'agg' (the total diple),
the second 'L2' layer and the last one 'L5' layer. For experimental
data, it can contain only one column.
nave : int
Number of trials that were averaged to produce this Dipole. Defaults
to 1
Attributes
----------
times : array-like
The time vector (in ms)
sfreq : float
The sampling frequency (in Hz)
data : dict of array
Dipole moment timecourse arrays with keys 'agg', 'L2' and 'L5'
nave : int
Number of trials that were averaged to produce this Dipole
scale_applied : int or float
The total factor by which the dipole has been scaled (using
:meth:`~hnn_core.dipole.Dipole.scale`).
"""
def __init__(self, times, data, nave=1): # noqa: D102
self.times = np.array(times)
if data.ndim == 1:
data = data[:, None]
if data.shape[1] == 3:
self.data = {'agg': data[:, 0], 'L2': data[:, 1], 'L5': data[:, 2]}
elif data.shape[1] == 1:
self.data = {'agg': data[:, 0]}
self.nave = nave
self.sfreq = 1000. / (times[1] - times[0]) # NB assumes len > 1
self.scale_applied = 1 # for visualisation
[docs] def copy(self):
"""Return a copy of the Dipole instance
Returns
-------
dpl_copy : instance of Dipole
A copy of the Dipole instance.
"""
return deepcopy(self)
def _post_proc(self, window_len, fctr):
"""Apply scaling and smoothing from param-files (DEPRECATE)
Parameters
----------
window_len : int
Smoothing window in ms
fctr : int
Scaling factor
"""
self.scale(fctr)
if window_len > 0: # this is to allow param-files with len==0
self.smooth(window_len)
def _convert_fAm_to_nAm(self):
"""The NEURON simulator output is in fAm, convert to nAm
NB! Must be run `after` :meth:`Dipole.baseline_renormalization`
"""
for key in self.data.keys():
self.data[key] *= 1e-6
[docs] def scale(self, factor):
"""Scale (multiply) the dipole moment by a fixed factor
The attribute ``Dipole.scale_applied`` is updated to reflect factors
applied and displayed in plots.
Parameters
----------
factor : int
Scaling factor, applied to the data in-place.
"""
for key in self.data.keys():
self.data[key] *= factor
self.scale_applied *= factor
return self
[docs] def smooth(self, window_len):
"""Smooth the dipole waveform using Hamming-windowed convolution
Note that this method operates in-place, i.e., it will alter the data.
If you prefer a filtered copy, consider using the
:meth:`~hnn_core.dipole.Dipole.copy`-method.
Parameters
----------
window_len : float
The length (in ms) of a `~numpy.hamming` window to convolve the
data with.
Returns
-------
dpl_copy : instance of Dipole
A copy of the modified Dipole instance.
"""
from .utils import smooth_waveform
for key in self.data.keys():
self.data[key] = smooth_waveform(self.data[key], window_len,
self.sfreq)
return self
[docs] def savgol_filter(self, h_freq):
"""Smooth the dipole waveform using Savitzky-Golay filtering
Note that this method operates in-place, i.e., it will alter the data.
If you prefer a filtered copy, consider using the
:meth:`~hnn_core.dipole.Dipole.copy`-method. The high-frequency cutoff
value of a Savitzky-Golay filter is approximate; see the SciPy
reference: :func:`~scipy.signal.savgol_filter`.
Parameters
----------
h_freq : float or None
Approximate high cutoff frequency in Hz. Note that this
is not an exact cutoff, since Savitzky-Golay filtering
is done using polynomial fits
instead of FIR/IIR filtering. This parameter is thus used to
determine the length of the window over which a 5th-order
polynomial smoothing is applied.
Returns
-------
dpl_copy : instance of Dipole
A copy of the modified Dipole instance.
"""
from .utils import _savgol_filter
if h_freq < 0:
raise ValueError('h_freq cannot be negative')
elif h_freq > 0.5 * self.sfreq:
raise ValueError(
'h_freq must be less than half the sample rate')
for key in self.data.keys():
self.data[key] = _savgol_filter(self.data[key],
h_freq,
self.sfreq)
return self
[docs] def plot(self, tmin=None, tmax=None, layer='agg', decim=None, ax=None,
color='k', show=True):
"""Simple layer-specific plot function.
Parameters
----------
tmin : float or None
Start time of plot (in ms). If None, plot entire simulation.
tmax : float or None
End time of plot (in ms). If None, plot entire simulation.
layer : str
The layer to plot. Can be one of 'agg', 'L2', and 'L5'
decimate : int
Factor by which to decimate the raw dipole traces (optional)
ax : instance of matplotlib figure | None
The matplotlib axis
color : tuple of float
RGBA value to use for plotting. By default, 'k' (black)
show : bool
If True, show the figure
Returns
-------
fig : instance of plt.fig
The matplotlib figure handle.
"""
return plot_dipole(self, tmin=tmin, tmax=tmax, ax=ax, layer=layer,
decim=decim, color=color, show=show)
[docs] def plot_psd(self, fmin=0, fmax=None, tmin=None, tmax=None, layer='agg',
color=None, label=None, ax=None, show=True):
"""Plot power spectral density (PSD) of dipole time course
Applies `~scipy.signal.periodogram` from SciPy with
``window='hamming'``.
Note that no spectral averaging is applied across time, as most
``hnn_core`` simulations are short-duration. However, passing a list of
`Dipole` instances will plot their average (Hamming-windowed) power,
which resembles the `Welch`-method applied over time.
Parameters
----------
dpl : instance of Dipole | list of Dipole instances
The Dipole object.
fmin : float
Minimum frequency to plot (in Hz). Default: 0 Hz
fmax : float
Maximum frequency to plot (in Hz). Default: None (plot up to
Nyquist)
tmin : float or None
Start time of data to include (in ms). If None, use entire
simulation.
tmax : float or None
End time of data to include (in ms). If None, use entire
simulation.
layer : str, default 'agg'
The layer to plot. Can be one of 'agg', 'L2', and 'L5'
color : str | tuple | None
The line color of PSD
label : str | None
Line label for PSD
ax : instance of matplotlib figure | None
The matplotlib axis.
show : bool
If True, show the figure
Returns
-------
fig : instance of matplotlib Figure
The matplotlib figure handle.
"""
return plot_psd(self, fmin=fmin, fmax=fmax, tmin=tmin, tmax=tmax,
layer=layer, color=color, label=label, ax=ax,
show=show)
[docs] def plot_tfr_morlet(self, freqs, n_cycles=7., tmin=None, tmax=None,
layer='agg', decim=None, padding='zeros', ax=None,
colormap='inferno', colorbar=True,
colorbar_inside=False, show=True):
"""Plot Morlet time-frequency representation of dipole time course
NB: Calls `~mne.time_frequency.tfr_array_morlet`, so ``mne`` must be
installed.
Parameters
----------
dpl : instance of Dipole | list of Dipole instances
The Dipole object. If a list of dipoles is given, the power is
calculated separately for each trial, then averaged.
freqs : array
Frequency range of interest.
n_cycles : float or array of float, default 7.0
Number of cycles. Fixed number or one per frequency.
tmin : float or None
Start time of plot in milliseconds. If None, plot entire
simulation.
tmax : float or None
End time of plot in milliseconds. If None, plot entire simulation.
layer : str, default 'agg'
The layer to plot. Can be one of 'agg', 'L2', and 'L5'
decim : int or list of int or None (default)
Optional (integer) factor by which to decimate the raw dipole
traces. The SciPy function :func:`~scipy.signal.decimate` is used,
which recommends values <13. To achieve higher decimation factors,
a list of ints can be provided. These are applied successively.
padding : str or None
Optional padding of the dipole time course beyond the plotting
limits. Possible values are: 'zeros' for padding with 0's
(default), 'mirror' for mirror-image padding.
ax : instance of matplotlib figure | None
The matplotlib axis
colormap : str
The name of a matplotlib colormap, e.g., 'viridis'. Default:
'inferno'
colorbar : bool
If True (default), adjust figure to include colorbar.
colorbar_inside: bool, default False
Put the color inside the heatmap if True.
show : bool
If True, show the figure
Returns
-------
fig : instance of matplotlib Figure
The matplotlib figure handle.
"""
return plot_tfr_morlet(
self, freqs, n_cycles=n_cycles, tmin=tmin, tmax=tmax,
layer=layer, decim=decim, padding=padding, ax=ax,
colormap=colormap, colorbar=colorbar,
colorbar_inside=colorbar_inside, show=show)
def _baseline_renormalize(self, N_pyr_x, N_pyr_y):
"""Only baseline renormalize if the units are fAm.
Parameters
----------
N_pyr_x : int
Nr of cells (x)
N_pyr_y : int
Nr of cells (y)
"""
# N_pyr cells in grid. This is PER LAYER
N_pyr = N_pyr_x * 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
self.data['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
# piecewise normalization
self.data['L5'][self.times <= 37.] -= dpl_offset['L5']
self.data['L5'][(self.times > 37.) & (self.times < t1)] -= N_pyr * \
(m * self.times[(self.times > 37.) & (self.times < t1)] + b)
self.data['L5'][self.times >= t1] -= N_pyr * \
(m1 * self.times[self.times >= t1] + b1)
# recalculate the aggregate dipole based on the baseline
# normalized ones
self.data['agg'] = self.data['L2'] + self.data['L5']
[docs] def write(self, fname):
"""Write dipole values to a file.
Parameters
----------
fname : str
Full path to the output file (.txt)
Outputs
-------
A tab separatd txt file where rows correspond
to samples and columns correspond to
1) time (s),
2) aggregate current dipole (scaled nAm),
3) L2/3 current dipole (scaled nAm), and
4) L5 current dipole (scaled nAm)
"""
if self.nave > 1:
warnings.warn("Saving Dipole to file that is an average of %d"
" trials" % self.nave)
X = [self.times]
fmt = ['%3.3f']
for data in self.data.values():
X.append(data)
fmt.append('%5.4f')
X = np.r_[X].T
np.savetxt(fname, X, fmt=fmt, delimiter='\t')