diff --git a/PointNICE/__init__.py b/PointNICE/__init__.py index c637f95..165fd6f 100644 --- a/PointNICE/__init__.py +++ b/PointNICE/__init__.py @@ -1,17 +1,16 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-06 13:36:00 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-29 17:38:00 +# @Last Modified time: 2018-04-03 20:21:05 ''' Import the core classes, generic utilities and algorithmic constants. ''' from .bls import BilayerSonophore from .solvers import * from .neurons import * from .utils import * from .constants import * from .plt import * -from .matlab_utils import * diff --git a/PointNICE/constants.py b/PointNICE/constants.py index 11f0233..3bf797d 100644 --- a/PointNICE/constants.py +++ b/PointNICE/constants.py @@ -1,54 +1,54 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-04 13:23:31 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-27 19:17:22 +# @Last Modified time: 2018-04-03 19:45:11 ''' Algorithmic constants used in the package. ''' # Fitting and pre-processing LJFIT_PM_MAX = 1e8 # intermolecular pressure at the deflection lower bound for LJ fitting (Pa) PNET_EQ_MAX = 1e-1 # error threshold for net pressure at computed equilibrium position (Pa) PMAVG_STD_ERR_MAX = 1500 # error threshold in nonlinear fit of molecular pressure (Pa) # Mechanical simulations Z_ERR_MAX = 1e-11 # periodic convergence threshold for deflection (m) NG_ERR_MAX = 1e-24 # periodic convergence threshold for gas content (mol) NCYCLES_MAX = 50 # max number of acoustic cycles in mechanical simulations CHARGE_RANGE = (-120e-5, 70e-5) # physiological charge range constraining the membrane (C/m2) # E-STIM simulations DT_ESTIM = 1e-4 # A-STIM simulations SOLVER_NSTEPS = 1000 # maximum number of steps allowed during one call to the LSODA/DOP853 solvers CLASSIC_TARGET_DT = 1e-8 # target temporal resolution for output arrays of classic simulations NPC_FULL = 1000 # nb of samples per acoustic period in full system NPC_HH = 40 # nb of samples per acoustic period in HH system DQ_UPDATE = 1e-5 # charge evolution threshold between two hybrid integrations (C/m2) DT_UPDATE = 5e-4 # time interval between two hybrid integrations (s) DT_EFF = 5e-5 # time step for effective integration (s) # Spike detection SPIKE_MIN_QAMP = 10e-5 # threshold amplitude for spike detection on charge signal (C/m2) -SPIKE_MIN_QPROM = 2e-5 # threshold prominence for spike detection on charge signal (C/m2) +SPIKE_MIN_QPROM = 5e-5 # threshold prominence for spike detection on charge signal (C/m2) SPIKE_MIN_VAMP = 10.0 # threshold amplitude for spike detection on potential signal (mV) SPIKE_MIN_DT = 1e-3 # minimal time interval for spike detection on charge signal (s) # Titrations TITRATION_T_OFFSET = 50e-3 # offset period for titration procedures (s) TITRATION_ASTIM_A_MAX = 3e5 # initial acoustic pressure upper bound for titration (Pa) TITRATION_ASTIM_DA_MAX = 1e3 # acoustic pressure search range threshold for titration (Pa) TITRATION_ESTIM_A_MAX = 50.0 # initial current density upper bound for titration (mA/m2) TITRATION_ESTIM_DA_MAX = 0.1 # current density search range threshold for titration (mA/m2) TITRATION_T_MAX = 2e-1 # initial stimulus duration upper bound for titration (s) TITRATION_DT_THR = 1e-3 # stimulus duration search range threshold for titration (s) TITRATION_DDC_THR = 0.01 # stimulus duty cycle search range threshold for titration (-) TITRATION_DC_MAX = 1.0 # initial stimulus duty cycle upper bound for titration (-) diff --git a/PointNICE/matlab_utils.py b/PointNICE/matlab_utils.py deleted file mode 100644 index 4da595b..0000000 --- a/PointNICE/matlab_utils.py +++ /dev/null @@ -1,69 +0,0 @@ -# -*- coding: utf-8 -*- -# @Author: Theo -# @Date: 2018-03-21 18:46:09 -# @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-29 18:57:42 - -''' Module of utilitary functions calling MATLAB routines internally. ''' - -import os -import logging -import pickle -import numpy as np -import matlab.engine -from .utils import InputError - -# Get package logger -logger = logging.getLogger('PointNICE') - - -def findPeaksMATLAB(y, mpp, mpd, mph): - ''' Find peaks in a signal by calling internally the matlab "finpeaks" function. - - :param y: input signal - :param mpp: minimal peak prominence (signal units) - :param mpd: minimal distance between consecutive peaks (number of indexes) - :param mph: minimal peak height (signal units) - :return: 4-tuple with indexes, values, witdths and prominences of peaks - ''' - - # Start MATLAB matlab_engine if not existing - if 'matlab_eng' not in globals(): - global matlab_eng - logger.info('starting matlab engine') - matlab_eng = matlab.engine.start_matlab() - logger.info('matlab engine started') - - _, ipeaks, widths, proms = np.array(matlab_eng.findpeaks(matlab.double(y.tolist()), - 'MinPeakProminence', mpp, - 'MinPeakDistance', mpd, - 'MinPeakHeight', mph, - nargout=4)) - return (np.array(ipeaks[0], dtype=int) - 1, np.array(widths[0]), np.array(proms[0])) - - -def detectSpikesMATLAB(t, y, mph, mpp, mtd): - ''' Detect spikes on a signal using the MATLAB findpeaks function, - and return their indexes, width and prominence. - - :param t: regular time vector (s) - :param y: signal vector (signal units) - :param mph: minimal peak height (signal units) - :param mpp: minimal peak prominence (signal units) - :param mtd: minimal distance between consecutive peaks (time difference) - :return: 3-tuple with indexes, widths and prominences of spikes - ''' - dt = t[1] - t[0] - ipeaks, widthpeaks, prompeaks = findPeaksMATLAB(y, mpp, int(np.ceil(mtd / dt)), mph) - return (ipeaks, widthpeaks * dt, prompeaks) - - -def detectSpikesInFile(filename, mph, mpp, mtd): - if not os.path.isfile(filename): - raise InputError('File does not exist') - with open(filename, 'rb') as fh: - frame = pickle.load(fh) - df = frame['data'] - t = df['t'].values - Qm = df['Qm'].values - return detectSpikesMATLAB(t, Qm, mph, mpp, mtd) diff --git a/PointNICE/plt/pltutils.py b/PointNICE/plt/pltutils.py index 7504e80..0ec49e2 100644 --- a/PointNICE/plt/pltutils.py +++ b/PointNICE/plt/pltutils.py @@ -1,963 +1,955 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-23 14:55:37 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-30 17:59:34 +# @Last Modified time: 2018-04-02 15:20:09 ''' Plotting utilities ''' import sys import os import pickle import ntpath import re import logging import tkinter as tk from tkinter import filedialog import numpy as np from scipy.interpolate import interp2d # from scipy.optimize import brentq import matplotlib import matplotlib.pyplot as plt from matplotlib.patches import Rectangle import matplotlib.cm as cm from matplotlib.ticker import FormatStrFormatter import pandas as pd from .. import neurons from ..utils import getNeuronsDict, getLookupDir, rescale, InputError from ..bls import BilayerSonophore from .pltvars import pltvars matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 # Get package logger logger = logging.getLogger('PointNICE') # Define global variables neuron = None bls = None timeunits = {'ASTIM': 't_ms', 'ESTIM': 't_ms', 'MECH': 't_us'} # Regular expression for input files rgxp = re.compile('(ESTIM|ASTIM)_([A-Za-z]*)_(.*).pkl') rgxp_mech = re.compile('(MECH)_(.*).pkl') # nb = '[0-9]*[.]?[0-9]+' # rgxp_ASTIM = re.compile('(ASTIM)_(\w+)_(PW|CW)_({0})nm_({0})kHz_({0})kPa_({0})ms(.*)_(\w+).pkl'.format(nb)) # rgxp_ESTIM = re.compile('(ESTIM)_(\w+)_(PW|CW)_({0})mA_per_m2_({0})ms(.*).pkl'.format(nb)) # rgxp_PW = re.compile('_PRF({0})kHz_DC({0})_(PW|CW)_(\d+)kHz_(\d+)kPa_(\d+)ms_(.*).pkl'.format(nb)) # Figure naming conventions ESTIM_CW_title = '{} neuron: CW E-STIM {:.2f}mA/m2, {:.0f}ms' ESTIM_PW_title = '{} neuron: PW E-STIM {:.2f}mA/m2, {:.0f}ms, {:.2f}kHz PRF, {:.0f}% DC' ASTIM_CW_title = '{} neuron: CW A-STIM {:.0f}kHz, {:.0f}kPa, {:.0f}ms' ASTIM_PW_title = '{} neuron: PW A-STIM {:.0f}kHz, {:.0f}kPa, {:.0f}ms, {:.2f}kHz PRF, {:.2f}% DC' MECH_title = '{:.0f}nm BLS structure: MECH-STIM {:.0f}kHz, {:.0f}kPa' class InteractiveLegend(object): """ Class defining an interactive matplotlib legend, where lines visibility can be toggled by simply clicking on the corresponding legend label. Other graphic objects can also be associated to the toggle of a specific line Adapted from: http://stackoverflow.com/questions/31410043/hiding-lines-after-showing-a-pyplot-figure """ def __init__(self, legend, aliases): self.legend = legend self.fig = legend.axes.figure self.lookup_artist, self.lookup_handle = self._build_lookups(legend) self._setup_connections() self.handles_aliases = aliases self.update() def _setup_connections(self): for artist in self.legend.texts + self.legend.legendHandles: artist.set_picker(10) # 10 points tolerance self.fig.canvas.mpl_connect('pick_event', self.on_pick) def _build_lookups(self, legend): ''' Method of the InteractiveLegend class building the legend lookups. ''' labels = [t.get_text() for t in legend.texts] handles = legend.legendHandles label2handle = dict(zip(labels, handles)) handle2text = dict(zip(handles, legend.texts)) lookup_artist = {} lookup_handle = {} for artist in legend.axes.get_children(): if artist.get_label() in labels: handle = label2handle[artist.get_label()] lookup_handle[artist] = handle lookup_artist[handle] = artist lookup_artist[handle2text[handle]] = artist lookup_handle.update(zip(handles, handles)) lookup_handle.update(zip(legend.texts, handles)) return lookup_artist, lookup_handle def on_pick(self, event): handle = event.artist if handle in self.lookup_artist: artist = self.lookup_artist[handle] artist.set_visible(not artist.get_visible()) self.update() def update(self): for artist in self.lookup_artist.values(): handle = self.lookup_handle[artist] if artist.get_visible(): handle.set_visible(True) if artist in self.handles_aliases: for al in self.handles_aliases[artist]: al.set_visible(True) else: handle.set_visible(False) if artist in self.handles_aliases: for al in self.handles_aliases[artist]: al.set_visible(False) self.fig.canvas.draw() def show(self): ''' showing the interactive legend ''' plt.show() def rescaleColorset(cset, lb=0, ub=255, lb_new=0., ub_new=1.): return list(tuple((x - lb) / (ub - lb) * (ub_new - lb_new) + lb_new for x in c) for c in cset) def getPatchesLoc(t, states): ''' Determine the location of stimulus patches. :param t: simulation time vector (s). :param states: a vector of stimulation state (ON/OFF) at each instant in time. :return: 3-tuple with number of patches, timing of STIM-ON an STIM-OFF instants. ''' # Compute states derivatives and identify bounds indexes of pulses dstates = np.diff(states) ipatch_on = np.insert(np.where(dstates > 0.0)[0] + 1, 0, 0) ipatch_off = np.where(dstates < 0.0)[0] # Get time instants for pulses ON and OFF npatches = ipatch_on.size tpatch_on = t[ipatch_on] tpatch_off = t[ipatch_off] # return 3-tuple with #patches, pulse ON and pulse OFF instants return (npatches, tpatch_on, tpatch_off) def SaveFigDialog(dirname, filename): """ Open a FileSaveDialogBox to set the directory and name of the figure to be saved. The default directory and filename are given, and the default extension is ".pdf" :param dirname: default directory :param filename: default filename :return: full path to the chosen filename """ root = tk.Tk() root.withdraw() filename_out = filedialog.asksaveasfilename(defaultextension=".pdf", initialdir=dirname, initialfile=filename) return filename_out def plotComp(varname, filepaths, labels=None, fs=15, lw=2, colors=None, lines=None, patches='one', xticks=None, yticks=None, blacklegend=False, straightlegend=False, showfig=True, inset=None): ''' Compare profiles of several specific output variables of NICE simulations. :param varname: name of variable to extract and compare :param filepaths: list of full paths to output data files to be compared :param labels: list of labels to use in the legend :param fs: labels fontsize :param patches: string indicating whether to indicate periods of stimulation with colored rectangular patches ''' # Input check 1: variable name if varname not in pltvars: raise InputError('Unknown plot variable: "{}"'.format(varname)) pltvar = pltvars[varname] # Input check 2: labels if labels: if len(labels) != len(filepaths): raise InputError('Invalid labels ({}): not matching number of compared files ({})' .format(len(labels), len(filepaths))) if not all(isinstance(x, str) for x in labels): raise InputError('Invalid labels: must be string typed') # Input check 3: line styles and colors if colors is None: colors = ['C{}'.format(j) for j in range(len(filepaths))] if lines is None: lines = ['-'] * len(filepaths) # Input check 4: STIM-ON patches greypatch = False if patches == 'none': patches = [False] * len(filepaths) elif patches == 'all': patches = [True] * len(filepaths) elif patches == 'one': patches = [True] + [False] * (len(filepaths) - 1) greypatch = True elif isinstance(patches, list): if len(patches) != len(filepaths): raise InputError('Invalid patches ({}): not matching number of compared files ({})' .format(len(patches), len(filepaths))) if not all(isinstance(p, bool) for p in patches): raise InputError('Invalid patch sequence: all list items must be boolean typed') else: raise InputError('Invalid patches: must be either "none", all", "one", or a boolean list') # Initialize figure and axis fig, ax = plt.subplots(figsize=(11, 4)) ax.set_zorder(0) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) if 'min' in pltvar and 'max' in pltvar: # optional min and max on y-axis ax.set_ylim(pltvar['min'], pltvar['max']) if pltvar['unit']: # y-label with optional unit ax.set_ylabel('$\\rm {}\ ({})$'.format(pltvar['label'], pltvar['unit']), fontsize=fs) else: ax.set_ylabel('$\\rm {}$'.format(pltvar['label']), fontsize=fs) if xticks is not None: # optional x-ticks ax.set_xticks(xticks) if yticks is not None: # optional y-ticks ax.set_yticks(yticks) else: ax.locator_params(axis='y', nbins=2) if any(ax.get_yticks() < 0): ax.yaxis.set_major_formatter(FormatStrFormatter('%+.0f')) for tick in ax.xaxis.get_major_ticks() + ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Optional inset axis if inset is not None: inset_ax = fig.add_axes(ax.get_position()) inset_ax.set_xlim(inset['xlims'][0], inset['xlims'][1]) inset_ax.set_ylim(inset['ylims'][0], inset['ylims'][1]) inset_ax.set_xticks([]) inset_ax.set_yticks([]) # inset_ax.patch.set_alpha(1.0) inset_ax.set_zorder(1) inset_ax.add_patch(Rectangle((inset['xlims'][0], inset['ylims'][0]), inset['xlims'][1] - inset['xlims'][0], inset['ylims'][1] - inset['ylims'][0], color='w')) # Retrieve neurons dictionary neurons_dict = getNeuronsDict() # Loop through data files aliases = {} for j, filepath in enumerate(filepaths): # Retrieve sim type pkl_filename = ntpath.basename(filepath) mo1 = rgxp.fullmatch(pkl_filename) mo2 = rgxp_mech.fullmatch(pkl_filename) if mo1: mo = mo1 elif mo2: mo = mo2 else: logger.error('Error: "%s" file does not match regexp pattern', pkl_filename) sys.exit(1) sim_type = mo.group(1) if sim_type not in ('MECH', 'ASTIM', 'ESTIM'): raise InputError('Invalid simulation type: {}'.format(sim_type)) if j == 0: sim_type_ref = sim_type t_plt = pltvars[timeunits[sim_type]] elif sim_type != sim_type_ref: raise InputError('Invalid comparison: different simulation types') # Load data logger.info('Loading data from "%s"', pkl_filename) with open(filepath, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] meta = frame['meta'] # Extract variables t = df['t'].values states = df['states'].values nsamples = t.size # Initialize neuron object if ESTIM or ASTIM sim type if sim_type in ['ASTIM', 'ESTIM']: neuron_name = mo.group(2) global neuron neuron = neurons_dict[neuron_name]() Cm0 = neuron.Cm0 Qm0 = Cm0 * neuron.Vm0 * 1e-3 # Extract neuron states if needed if 'alias' in pltvar and 'neuron_states' in pltvar['alias']: neuron_states = [df[sn].values for sn in neuron.states_names] else: Cm0 = meta['Cm0'] Qm0 = meta['Qm0'] # Initialize BLS if needed if sim_type in ['MECH', 'ASTIM'] and 'alias' in pltvar and 'bls' in pltvar['alias']: global bls bls = BilayerSonophore(meta['a'], meta['Fdrive'], Cm0, Qm0) # Determine patches location npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) # Add onset to time and states vectors if t_plt['onset'] > 0.0: t = np.insert(t, 0, -t_plt['onset']) states = np.insert(states, 0, 0) # Set x-axis label ax.set_xlabel('$\\rm {}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) # Extract variable to plot if 'alias' in pltvar: var = eval(pltvar['alias']) elif 'key' in pltvar: var = df[pltvar['key']].values elif 'constant' in pltvar: var = eval(pltvar['constant']) * np.ones(nsamples) else: var = df[varname].values if var.size == t.size - 1: var = np.insert(var, 0, var[0]) # Determine legend label if labels: label = labels[j] else: if sim_type == 'ESTIM': if meta['DC'] == 1.0: label = ESTIM_CW_title.format(neuron_name, meta['Astim'], meta['tstim'] * 1e3) else: label = ESTIM_PW_title.format(neuron_name, meta['Astim'], meta['tstim'] * 1e3, meta['PRF'] * 1e-3, meta['DC'] * 1e2) elif sim_type == 'ASTIM': if meta['DC'] == 1.0: label = ASTIM_CW_title.format(neuron_name, meta['Fdrive'] * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3) else: label = ASTIM_PW_title.format(neuron_name, meta['Fdrive'] * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3, meta['PRF'] * 1e-3, meta['DC'] * 1e2) elif sim_type == 'MECH': label = MECH_title.format(meta['a'] * 1e9, meta['Fdrive'] * 1e-3, meta['Adrive'] * 1e-3) # Plot trace handle = ax.plot(t * t_plt['factor'], var * pltvar['factor'], linewidth=lw, linestyle=lines[j], color=colors[j], label=label) if inset is not None: inset_window = np.logical_and(t > (inset['xlims'][0] / t_plt['factor']), t < (inset['xlims'][1] / t_plt['factor'])) inset_ax.plot(t[inset_window] * t_plt['factor'], var[inset_window] * pltvar['factor'], linewidth=lw, linestyle=lines[j], color=colors[j]) # Add optional STIM-ON patches if patches[j]: (ybottom, ytop) = ax.get_ylim() la = [] color = '#8A8A8A' if greypatch else handle[0].get_color() for i in range(npatches): la.append(ax.add_patch(Rectangle((tpatch_on[i] * t_plt['factor'], ybottom), (tpatch_off[i] - tpatch_on[i]) * t_plt['factor'], ytop - ybottom, color=color, alpha=0.1))) aliases[handle[0]] = la if inset is not None: cond_on = np.logical_and(tpatch_on > (inset['xlims'][0] / t_plt['factor']), tpatch_on < (inset['xlims'][1] / t_plt['factor'])) cond_off = np.logical_and(tpatch_off > (inset['xlims'][0] / t_plt['factor']), tpatch_off < (inset['xlims'][1] / t_plt['factor'])) cond_glob = np.logical_and(tpatch_on < (inset['xlims'][0] / t_plt['factor']), tpatch_off > (inset['xlims'][1] / t_plt['factor'])) cond_onoff = np.logical_or(cond_on, cond_off) cond = np.logical_or(cond_onoff, cond_glob) npatches_inset = np.sum(cond) for i in range(npatches_inset): inset_ax.add_patch(Rectangle((tpatch_on[cond][i] * t_plt['factor'], ybottom), (tpatch_off[cond][i] - tpatch_on[cond][i]) * t_plt['factor'], ytop - ybottom, color=color, alpha=0.1)) fig.tight_layout() # Optional operations on inset: if inset is not None: - # Re-position inset axis and materialize it with a dashed frame + # Re-position inset axis axpos = ax.get_position() left, right, = rescale(inset['xcoords'], ax.get_xlim()[0], ax.get_xlim()[1], axpos.x0, axpos.x0 + axpos.width) bottom, top, = rescale(inset['ycoords'], ax.get_ylim()[0], ax.get_ylim()[1], axpos.y0, axpos.y0 + axpos.height) inset_ax.set_position([left, bottom, right - left, top - bottom]) - inset_ax.axis('off') - inner_factor = 0.02 - dx = inset['xlims'][1] - inset['xlims'][0] - innerx = [inset['xlims'][0] + inner_factor * dx, inset['xlims'][1] - inner_factor * dx] - dy = inset['ylims'][1] - inset['ylims'][0] - innery = [inset['ylims'][0] + inner_factor * dy, inset['ylims'][1] - inner_factor * dy] - inset_ax.plot(innerx, [innery[0]] * 2, linestyle='--', color='k') - inset_ax.plot(innerx, [innery[1]] * 2, linestyle='--', color='k') - inset_ax.plot([innerx[0]] * 2, innery, linestyle='--', color='k') - inset_ax.plot([innerx[1]] * 2, innery, linestyle='--', color='k') - - # Materialize inset target region with dashed frame - ax.plot(inset['xlims'], [inset['ylims'][0]] * 2, linestyle='--', color='k') - ax.plot(inset['xlims'], [inset['ylims'][1]] * 2, linestyle='--', color='k') - ax.plot([inset['xlims'][0]] * 2, inset['ylims'], linestyle='--', color='k') - ax.plot([inset['xlims'][1]] * 2, inset['ylims'], linestyle='--', color='k') + for i in inset_ax.spines.values(): + i.set_linewidth(2) + + # Materialize inset target region with contour frame + ax.plot(inset['xlims'], [inset['ylims'][0]] * 2, linestyle='-', color='k') + ax.plot(inset['xlims'], [inset['ylims'][1]] * 2, linestyle='-', color='k') + ax.plot([inset['xlims'][0]] * 2, inset['ylims'], linestyle='-', color='k') + ax.plot([inset['xlims'][1]] * 2, inset['ylims'], linestyle='-', color='k') # Link target and inset with dashed lines if possible if inset['xcoords'][1] < inset['xlims'][0]: ax.plot([inset['xcoords'][1], inset['xlims'][0]], [inset['ycoords'][0], inset['ylims'][0]], linestyle='--', color='k') ax.plot([inset['xcoords'][1], inset['xlims'][0]], [inset['ycoords'][1], inset['ylims'][1]], linestyle='--', color='k') elif inset['xcoords'][0] > inset['xlims'][1]: ax.plot([inset['xcoords'][0], inset['xlims'][1]], [inset['ycoords'][0], inset['ylims'][0]], linestyle='--', color='k') ax.plot([inset['xcoords'][0], inset['xlims'][1]], [inset['ycoords'][1], inset['ylims'][1]], linestyle='--', color='k') else: logger.warning('Inset x-coordinates intersect with those of target region') # Create interactive legend leg = ax.legend(loc=1, fontsize=fs, frameon=False) if blacklegend: for l in leg.get_lines(): l.set_color('k') if straightlegend: for l in leg.get_lines(): l.set_linestyle('-') interactive_legend = InteractiveLegend(ax.legend_, aliases) if showfig: plt.show() return fig def plotBatch(directory, filepaths, vars_dict=None, plt_show=True, plt_save=False, ask_before_save=True, fig_ext='png', tag='fig', fs=15, lw=2, title=True, show_patches=True): ''' Plot a figure with profiles of several specific NICE output variables, for several NICE simulations. :param positions: subplot indexes of each variable :param filepaths: list of full paths to output data files to be compared :param vars_dict: dict of lists of variables names to extract and plot together :param plt_show: boolean stating whether to show the created figures :param plt_save: boolean stating whether to save the created figures :param ask_before_save: boolean stating whether to show the created figures :param fig_ext: file extension for the saved figures :param tag: suffix added to the end of the figures name :param fs: labels font size :param lw: curves line width :param title: boolean stating whether to display a general title on the figures :param show_patches: boolean indicating whether to indicate periods of stimulation with colored rectangular patches ''' # Check validity of plot variables if vars_dict: yvars = list(sum(list(vars_dict.values()), [])) for key in yvars: if key not in pltvars: raise InputError('Unknown plot variable: "{}"'.format(key)) # Dictionary of neurons neurons_dict = getNeuronsDict() # Loop through data files for filepath in filepaths: # Get code from file name pkl_filename = ntpath.basename(filepath) filecode = pkl_filename[0:-4] # Retrieve sim type mo1 = rgxp.fullmatch(pkl_filename) mo2 = rgxp_mech.fullmatch(pkl_filename) if mo1: mo = mo1 elif mo2: mo = mo2 else: logger.error('Error: "%s" file does not match regexp pattern', pkl_filename) sys.exit(1) sim_type = mo.group(1) if sim_type not in ('MECH', 'ASTIM', 'ESTIM'): raise InputError('Invalid simulation type: {}'.format(sim_type)) # Load data logger.info('Loading data from "%s"', pkl_filename) with open(filepath, 'rb') as fh: frame = pickle.load(fh) df = frame['data'] meta = frame['meta'] # Extract variables logger.info('Extracting variables') t = df['t'].values states = df['states'].values nsamples = t.size # Initialize channel mechanism if sim_type in ['ASTIM', 'ESTIM']: neuron_name = mo.group(2) global neuron neuron = neurons_dict[neuron_name]() neuron_states = [df[sn].values for sn in neuron.states_names] Cm0 = neuron.Cm0 Qm0 = Cm0 * neuron.Vm0 * 1e-3 t_plt = pltvars['t_ms'] else: Cm0 = meta['Cm0'] Qm0 = meta['Qm0'] t_plt = pltvars['t_us'] # Initialize BLS if sim_type in ['MECH', 'ASTIM']: global bls Fdrive = meta['Fdrive'] a = meta['a'] bls = BilayerSonophore(a, Fdrive, Cm0, Qm0) # Determine patches location npatches, tpatch_on, tpatch_off = getPatchesLoc(t, states) # Adding onset to time and states vectors if t_plt['onset'] > 0.0: t = np.insert(t, 0, -t_plt['onset']) states = np.insert(states, 0, 0) # Determine variables to plot if not provided if not vars_dict: if sim_type == 'ASTIM': vars_dict = {'Z': ['Z'], 'Q_m': ['Qm']} elif sim_type == 'ESTIM': vars_dict = {'V_m': ['Vm']} elif sim_type == 'MECH': vars_dict = {'P_{AC}': ['Pac'], 'Z': ['Z'], 'n_g': ['ng']} if sim_type in ['ASTIM', 'ESTIM'] and hasattr(neuron, 'pltvars_scheme'): vars_dict.update(neuron.pltvars_scheme) labels = list(vars_dict.keys()) naxes = len(vars_dict) # Plotting if naxes == 1: _, ax = plt.subplots(figsize=(11, 4)) axes = [ax] else: _, axes = plt.subplots(naxes, 1, figsize=(11, min(3 * naxes, 9))) for i in range(naxes): ax = axes[i] ax_pltvars = [pltvars[j] for j in vars_dict[labels[i]]] nvars = len(ax_pltvars) # X-axis if i < naxes - 1: ax.get_xaxis().set_ticklabels([]) else: ax.set_xlabel('${}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Y-axis if ax_pltvars[0]['unit']: ax.set_ylabel('${}\ ({})$'.format(labels[i], ax_pltvars[0]['unit']), fontsize=fs) else: ax.set_ylabel('${}$'.format(labels[i]), fontsize=fs) if 'min' in ax_pltvars[0] and 'max' in ax_pltvars[0]: ax_min = min([ap['min'] for ap in ax_pltvars]) ax_max = max([ap['max'] for ap in ax_pltvars]) ax.set_ylim(ax_min, ax_max) ax.locator_params(axis='y', nbins=2) for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Time series icolor = 0 for j in range(nvars): # Extract variable pltvar = ax_pltvars[j] if 'alias' in pltvar: var = eval(pltvar['alias']) elif 'key' in pltvar: var = df[pltvar['key']].values elif 'constant' in pltvar: var = eval(pltvar['constant']) * np.ones(nsamples) else: var = df[vars_dict[labels[i]][j]].values if var.size == t.size - 1: var = np.insert(var, 0, var[0]) # Plot variable if 'constant' in pltvar: ax.plot(t * t_plt['factor'], var * pltvar['factor'], '--', c='black', lw=lw, label='${}$'.format(pltvar['label'])) else: ax.plot(t * t_plt['factor'], var * pltvar['factor'], c='C{}'.format(icolor), lw=lw, label='${}$'.format(pltvar['label'])) icolor += 1 # Patches if show_patches == 1: (ybottom, ytop) = ax.get_ylim() for j in range(npatches): ax.add_patch(Rectangle((tpatch_on[j] * t_plt['factor'], ybottom), (tpatch_off[j] - tpatch_on[j]) * t_plt['factor'], ytop - ybottom, color='#8A8A8A', alpha=0.1)) # Legend if nvars > 1: ax.legend(fontsize=fs, loc=7, ncol=nvars // 4 + 1) # Title if title: if sim_type == 'ESTIM': if meta['DC'] == 1.0: fig_title = ESTIM_CW_title.format(neuron.name, meta['Astim'], meta['tstim'] * 1e3) else: fig_title = ESTIM_PW_title.format(neuron.name, meta['Astim'], meta['tstim'] * 1e3, meta['PRF'] * 1e-3, meta['DC'] * 1e2) elif sim_type == 'ASTIM': if meta['DC'] == 1.0: fig_title = ASTIM_CW_title.format(neuron.name, Fdrive * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3) else: fig_title = ASTIM_PW_title.format(neuron.name, Fdrive * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3, meta['PRF'] * 1e-3, meta['DC'] * 1e2) elif sim_type == 'MECH': fig_title = MECH_title.format(a * 1e9, Fdrive * 1e-3, meta['Adrive'] * 1e-3) axes[0].set_title(fig_title, fontsize=fs) plt.tight_layout() # Save figure if needed (automatic or checked) if plt_save: if ask_before_save: plt_filename = SaveFigDialog(directory, '{}_{}.{}'.format(filecode, tag, fig_ext)) else: plt_filename = '{}/{}_{}.{}'.format(directory, filecode, tag, fig_ext) if plt_filename: plt.savefig(plt_filename) logger.info('Saving figure as "{}"'.format(plt_filename)) plt.close() # Show all plots if needed if plt_show: plt.show() def plotGatingKinetics(neuron, fs=15): ''' Plot the voltage-dependent steady-states and time constants of activation and inactivation gates of the different ionic currents involved in a specific neuron's membrane. :param neuron: specific channel mechanism object :param fs: labels and title font size ''' # Input membrane potential vector Vm = np.linspace(-100, 50, 300) xinf_dict = {} taux_dict = {} logger.info('Computing %s neuron gating kinetics', neuron.name) names = neuron.states_names print(names) for xname in names: Vm_state = True # Names of functions of interest xinf_func_str = xname.lower() + 'inf' taux_func_str = 'tau' + xname.lower() alphax_func_str = 'alpha' + xname.lower() betax_func_str = 'beta' + xname.lower() # derx_func_str = 'der' + xname.upper() # 1st choice: use xinf and taux function if hasattr(neuron, xinf_func_str) and hasattr(neuron, taux_func_str): xinf_func = getattr(neuron, xinf_func_str) taux_func = getattr(neuron, taux_func_str) xinf = np.array([xinf_func(v) for v in Vm]) if isinstance(taux_func, float): taux = taux_func * np.ones(len(Vm)) else: taux = np.array([taux_func(v) for v in Vm]) # 2nd choice: use alphax and betax functions elif hasattr(neuron, alphax_func_str) and hasattr(neuron, betax_func_str): alphax_func = getattr(neuron, alphax_func_str) betax_func = getattr(neuron, betax_func_str) alphax = np.array([alphax_func(v) for v in Vm]) if isinstance(betax_func, float): betax = betax_func * np.ones(len(Vm)) else: betax = np.array([betax_func(v) for v in Vm]) taux = 1.0 / (alphax + betax) xinf = taux * alphax # # 3rd choice: use derX choice # elif hasattr(neuron, derx_func_str): # derx_func = getattr(neuron, derx_func_str) # xinf = brentq(lambda x: derx_func(neuron.Vm, x), 0, 1) else: Vm_state = False if not Vm_state: logger.error('no function to compute %s-state gating kinetics', xname) else: xinf_dict[xname] = xinf taux_dict[xname] = taux fig, axes = plt.subplots(2) fig.suptitle('{} neuron: gating dynamics'.format(neuron.name)) ax = axes[0] ax.get_xaxis().set_ticklabels([]) ax.set_ylabel('$X_{\infty}$', fontsize=fs) for xname in names: if xname in xinf_dict: ax.plot(Vm, xinf_dict[xname], lw=2, label='$' + xname + '_{\infty}$') ax.legend(fontsize=fs, loc=7) ax = axes[1] ax.set_xlabel('$V_m\ (mV)$', fontsize=fs) ax.set_ylabel('$\\tau_X\ (ms)$', fontsize=fs) for xname in names: if xname in taux_dict: ax.plot(Vm, taux_dict[xname] * 1e3, lw=2, label='$\\tau_{' + xname + '}$') ax.legend(fontsize=fs, loc=7) plt.show() def plotRateConstants(neuron, fs=15): ''' Plot the voltage-dependent activation and inactivation rate constants for each gate of all ionic currents involved in a specific neuron's membrane. :param neuron: specific channel mechanism object :param fs: labels and title font size ''' # Input membrane potential vector Vm = np.linspace(-250, 250, 300) alphax_dict = {} betax_dict = {} logger.info('Computing %s neuron gating kinetics', neuron.name) names = neuron.states_names for xname in names: Vm_state = True # Names of functions of interest xinf_func_str = xname.lower() + 'inf' taux_func_str = 'tau' + xname.lower() alphax_func_str = 'alpha' + xname.lower() betax_func_str = 'beta' + xname.lower() # 1st choice: use alphax and betax functions if hasattr(neuron, alphax_func_str) and hasattr(neuron, betax_func_str): alphax_func = getattr(neuron, alphax_func_str) betax_func = getattr(neuron, betax_func_str) alphax = np.array([alphax_func(v) for v in Vm]) betax = np.array([betax_func(v) for v in Vm]) # 2nd choice: use xinf and taux function elif hasattr(neuron, xinf_func_str) and hasattr(neuron, taux_func_str): xinf_func = getattr(neuron, xinf_func_str) taux_func = getattr(neuron, taux_func_str) xinf = np.array([xinf_func(v) for v in Vm]) taux = np.array([taux_func(v) for v in Vm]) alphax = xinf / taux betax = 1.0 / taux - alphax else: Vm_state = False if not Vm_state: logger.error('no function to compute %s-state gating kinetics', xname) else: alphax_dict[xname] = alphax betax_dict[xname] = betax naxes = len(alphax_dict) _, axes = plt.subplots(naxes, figsize=(11, min(3 * naxes, 9))) for i, xname in enumerate(alphax_dict.keys()): ax1 = axes[i] if i == 0: ax1.set_title('{} neuron: rate constants'.format(neuron.name)) if i == naxes - 1: ax1.set_xlabel('$V_m\ (mV)$', fontsize=fs) else: ax1.get_xaxis().set_ticklabels([]) ax1.set_ylabel('$\\alpha_{' + xname + '}\ (ms^{-1})$', fontsize=fs, color='C0') for label in ax1.get_yticklabels(): label.set_color('C0') ax1.plot(Vm, alphax_dict[xname] * 1e-3, lw=2) ax2 = ax1.twinx() ax2.set_ylabel('$\\beta_{' + xname + '}\ (ms^{-1})$', fontsize=fs, color='C1') for label in ax2.get_yticklabels(): label.set_color('C1') ax2.plot(Vm, betax_dict[xname] * 1e-3, lw=2, color='C1') plt.tight_layout() plt.show() def setGrid(n): ''' Determine number of rows and columns in figure grid, based on number of variables to plot. ''' if n <= 3: return (1, n) else: return ((n - 1) // 3 + 1, 3) def plotEffCoeffs(neuron, Fdrive, a=32e-9, fs=12): ''' Plot the profiles of all effective coefficients of a specific neuron for a given frequency. For each coefficient X, one line chart per lookup amplitude is plotted, using charge as the input variable on the abscissa and a linear color code for the amplitude value. :param neuron: channel mechanism object :param Fdrive: acoustic drive frequency (Hz) :param a: sonophore diameter (m) :param fs: figure fontsize ''' # Check lookup file existence lookup_file = '{}_lookups_a{:.1f}nm.pkl'.format(neuron.name, a * 1e9) lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) if not os.path.isfile(lookup_path): raise InputError('Missing lookup file: "{}"'.format(lookup_file)) # Load coefficients with open(lookup_path, 'rb') as fh: lookup_dict = pickle.load(fh) # Retrieve 1D inputs from lookup dictionary freqs = lookup_dict['f'] amps = lookup_dict['A'] charges = lookup_dict['Q'] # Check that frequency is within lookup range margin = 1e-9 # adding margin to compensate for eventual round error frange = (freqs.min() - margin, freqs.max() + margin) if Fdrive < frange[0] or Fdrive > frange[1]: raise InputError(('Invalid frequency: {:.2f} kHz (must be within ' + '{:.1f} kHz - {:.1f} MHz lookup interval)') .format(Fdrive * 1e-3, frange[0] * 1e-3, frange[1] * 1e-6)) # Define coefficients list coeffs_list = ['V', 'ng', *neuron.coeff_names] AQ_coeffs = {} # If Fdrive in lookup frequencies, simply project (A, Q) dataset at that frequency if Fdrive in freqs: iFdrive = np.searchsorted(freqs, Fdrive) logger.info('Using lookups directly at %.2f kHz', freqs[iFdrive] * 1e-3) for cn in coeffs_list: AQ_coeffs[cn] = np.squeeze(lookup_dict[cn][iFdrive, :, :]) # Otherwise, project 2 (A, Q) interpolation datasets at Fdrive bounding values # indexes in lookup frequencies onto two 1D charge-based interpolation datasets, and # interpolate between them afterwards else: ilb = np.searchsorted(freqs, Fdrive) - 1 logger.info('Interpolating lookups between %.2f kHz and %.2f kHz', freqs[ilb] * 1e-3, freqs[ilb + 1] * 1e-3) for cn in coeffs_list: AQ_slice = [] for iAdrive in range(len(amps)): fQ_slice = np.squeeze(lookup_dict[cn][:, iAdrive, :]) itrp = interp2d(freqs, charges, fQ_slice.T) Q_vect = itrp(Fdrive, charges) AQ_slice.append(Q_vect) AQ_coeffs[cn] = np.squeeze(np.array([AQ_slice])) # Replace dict key AQ_coeffs['Veff'] = AQ_coeffs.pop('V') # Plotting logger.info('plotting') Amin, Amax = amps.min(), amps.max() ncoeffs = len(coeffs_list) nrows, ncols = setGrid(ncoeffs) xvar = pltvars['Qm'] mymap = cm.get_cmap('viridis') sm_amp = plt.cm.ScalarMappable(cmap=mymap, norm=plt.Normalize(Amin * 1e-3, Amax * 1e-3)) sm_amp._A = [] fig, _ = plt.subplots(figsize=(5 * ncols, 1.8 * nrows), squeeze=False) for j, cn in enumerate(AQ_coeffs.keys()): ax = plt.subplot2grid((nrows, ncols), (j // 3, j % 3)) # ax = axes[j // 3, j % 3] yvar = pltvars[cn] ax.set_xlabel('${}\ ({})$'.format(xvar['label'], xvar['unit']), fontsize=fs) ax.set_ylabel('${}\ ({})$'.format(yvar['label'], yvar['unit']), fontsize=fs) for i, Adrive in enumerate(amps): ax.plot(charges * xvar['factor'], AQ_coeffs[cn][i, :] * yvar['factor'], c=mymap(rescale(Adrive, Amin, Amax))) plt.tight_layout() fig.suptitle('{} neuron: effective coefficients @ {:.0f} kHz'.format( neuron.name, Fdrive * 1e-3)) fig.subplots_adjust(top=0.9, right=0.85) cbar_ax = fig.add_axes([0.87, 0.05, 0.02, 0.9]) fig.add_axes() fig.colorbar(sm_amp, cax=cbar_ax) cbar_ax.set_ylabel('$A_{drive} \ (kPa)$', fontsize=fs) plt.show() diff --git a/PointNICE/solvers/simutils.py b/PointNICE/solvers/simutils.py index 9313861..569ff66 100644 --- a/PointNICE/solvers/simutils.py +++ b/PointNICE/solvers/simutils.py @@ -1,1249 +1,1403 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-22 14:33:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-29 19:20:45 +# @Last Modified time: 2018-04-03 20:11:06 """ Utility functions used in simulations """ import os import time import logging import pickle import shutil import tkinter as tk from tkinter import filedialog import numpy as np from openpyxl import load_workbook import lockfile import pandas as pd from ..bls import BilayerSonophore from .SolverUS import SolverUS from .SolverElec import SolverElec from ..constants import * from ..utils import getNeuronsDict, InputError # Get package logger logger = logging.getLogger('PointNICE') # Naming nomenclature for output files MECH_code = 'MECH_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.1f}nCcm2' ESTIM_CW_code = 'ESTIM_{}_CW_{:.1f}mA_per_m2_{:.0f}ms' ESTIM_PW_code = 'ESTIM_{}_PW_{:.1f}mA_per_m2_{:.0f}ms_PRF{:.2f}kHz_DC{:.2f}' ASTIM_CW_code = 'ASTIM_{}_CW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_{}' ASTIM_PW_code = 'ASTIM_{}_PW_{:.0f}nm_{:.0f}kHz_{:.0f}kPa_{:.0f}ms_PRF{:.2f}kHz_DC{:.2f}_{}' # Parameters units ASTIM_params = { 'f': {'index': 0, 'factor': 1e-3, 'unit': 'kHz'}, 'A': {'index': 1, 'factor': 1e-3, 'unit': 'kPa'}, 't': {'index': 2, 'factor': 1e3, 'unit': 'ms'}, 'PRF': {'index': 4, 'factor': 1e-3, 'unit': 'kHz'}, 'DC': {'index': 5, 'factor': 1e2, 'unit': '%'} } ESTIM_params = { 'A': {'index': 0, 'factor': 1e0, 'unit': 'mA/m2'}, 't': {'index': 1, 'factor': 1e3, 'unit': 'ms'}, 'PRF': {'index': 3, 'factor': 1e-3, 'unit': 'kHz'}, 'DC': {'index': 4, 'factor': 1e2, 'unit': '%'} } # Default geometry default_diam = 32e-9 default_embedding = 0.0e-6 def setBatchDir(): ''' Select batch directory for output files.α :return: full path to batch directory ''' root = tk.Tk() root.withdraw() batch_dir = filedialog.askdirectory() if not batch_dir: raise InputError('No output directory chosen') return batch_dir def checkBatchLog(batch_dir, batch_type): ''' Check for appropriate log file in batch directory, and create one if it is absent. :param batch_dir: full path to batch directory :param batch_type: type of simulation batch :return: 2 tuple with full path to log file and boolean stating if log file was created ''' # Check for directory existence if not os.path.isdir(batch_dir): raise InputError('"{}" output directory does not exist'.format(batch_dir)) # Determine log template from batch type if batch_type == 'MECH': logfile = 'log_MECH.xlsx' elif batch_type == 'A-STIM': logfile = 'log_ASTIM.xlsx' elif batch_type == 'E-STIM': logfile = 'log_ESTIM.xlsx' else: raise InputError('Unknown batch type', batch_type) # Get template in package subdirectory this_dir, _ = os.path.split(__file__) parent_dir = os.path.abspath(os.path.join(this_dir, os.pardir)) logsrc = parent_dir + '/templates/' + logfile assert os.path.isfile(logsrc), 'template log file "{}" not found'.format(logsrc) # Copy template in batch directory if no appropriate log file logdst = batch_dir + '/' + logfile is_log = os.path.isfile(logdst) if not is_log: shutil.copy2(logsrc, logdst) return (logdst, not is_log) def createSimQueue(amps, durations, offsets, PRFs, DCs): ''' Create a serialized 2D array of all parameter combinations for a series of individual parameter sweeps, while avoiding repetition of CW protocols for a given PRF sweep. :param amps: list (or 1D-array) of acoustic amplitudes :param durations: list (or 1D-array) of stimulus durations :param offsets: list (or 1D-array) of stimulus offsets (paired with durations array) :param PRFs: list (or 1D-array) of pulse-repetition frequencies :param DCs: list (or 1D-array) of duty cycle values :return: 2D-array with (amplitude, duration, offset, PRF, DC) for each stimulation protocol ''' # Convert input to 1D-arrays amps = np.array(amps) durations = np.array(durations) offsets = np.array(offsets) PRFs = np.array(PRFs) DCs = np.array(DCs) # Create index arrays iamps = range(len(amps)) idurs = range(len(durations)) # Create empty output matrix queue = np.empty((1, 5)) # Continuous protocols if 1.0 in DCs: nCW = len(amps) * len(durations) arr1 = np.ones(nCW) iCW_queue = np.array(np.meshgrid(iamps, idurs)).T.reshape(nCW, 2) CW_queue = np.vstack((amps[iCW_queue[:, 0]], durations[iCW_queue[:, 1]], offsets[iCW_queue[:, 1]], PRFs.min() * arr1, arr1)).T queue = np.vstack((queue, CW_queue)) # Pulsed protocols if np.any(DCs != 1.0): pulsed_DCs = DCs[DCs != 1.0] iPRFs = range(len(PRFs)) ipulsed_DCs = range(len(pulsed_DCs)) nPW = len(amps) * len(durations) * len(PRFs) * len(pulsed_DCs) iPW_queue = np.array(np.meshgrid(iamps, idurs, iPRFs, ipulsed_DCs)).T.reshape(nPW, 4) PW_queue = np.vstack((amps[iPW_queue[:, 0]], durations[iPW_queue[:, 1]], offsets[iPW_queue[:, 1]], PRFs[iPW_queue[:, 2]], pulsed_DCs[iPW_queue[:, 3]])).T queue = np.vstack((queue, PW_queue)) # Return return queue[1:, :] def xlslog(filename, sheetname, data): """ Append log data on a new row to specific sheet of excel workbook, using a lockfile to avoid read/write errors between concurrent processes. :param filename: absolute or relative path to the Excel workbook :param sheetname: name of the Excel spreadsheet to which data is appended :param data: data structure to be added to specific columns on a new row :return: boolean indicating success (1) or failure (0) of operation """ try: lock = lockfile.FileLock(filename) lock.acquire() wb = load_workbook(filename) ws = wb[sheetname] keys = data.keys() i = 1 row_data = {} for k in keys: row_data[k] = data[k] i += 1 ws.append(row_data) wb.save(filename) lock.release() return 1 except PermissionError: # If file cannot be accessed for writing because already opened logger.warning('Cannot write to "%s". Close the file and type "Y"', filename) user_str = input() if user_str in ['y', 'Y']: return xlslog(filename, sheetname, data) else: return 0 def detectPeaks(x, mph=None, mpd=1, threshold=0, edge='rising', kpsh=False, valley=False, ax=None): ''' Detect peaks in data based on their amplitude and other features. - From Marco Duarte: + Adapted from Marco Duarte: http://nbviewer.jupyter.org/github/demotu/BMC/blob/master/notebooks/DetectPeaks.ipynb :param x: 1D array_like data. :param mph: minimum peak height (default = None). :param mpd: minimum peak distance in indexes (default = 1) :param threshold : minimum peak prominence (default = 0) :param edge : for a flat peak, keep only the rising edge ('rising'), only the falling edge ('falling'), both edges ('both'), or don't detect a flat peak (None). (default = 'rising') :param kpsh: keep peaks with same height even if they are closer than `mpd` (default = False). :param valley: detect valleys (local minima) instead of peaks (default = False). :param show: plot data in matplotlib figure (default = False). :param ax: a matplotlib.axes.Axes instance, optional (default = None). :return: 1D array with the indices of the peaks ''' + print('min peak height:', mph, ', min peak distance:', mpd, + ', min peak prominence:', threshold) + # Convert input to numpy array x = np.atleast_1d(x).astype('float64') - if x.size < 3: - return np.array([], dtype=int) + + # Revert signal sign for valley detection if valley: x = -x - # find indices of all peaks - dx = x[1:] - x[:-1] - # handle NaN's - indnan = np.where(np.isnan(x))[0] - if indnan.size: - x[indnan] = np.inf - dx[np.where(np.isnan(dx))[0]] = np.inf + + # Differentiate signal + dx = np.diff(x) + + # Find indices of all peaks with edge criterion ine, ire, ife = np.array([[], [], []], dtype=int) if not edge: ine = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) > 0))[0] else: if edge.lower() in ['rising', 'both']: ire = np.where((np.hstack((dx, 0)) <= 0) & (np.hstack((0, dx)) > 0))[0] if edge.lower() in ['falling', 'both']: ife = np.where((np.hstack((dx, 0)) < 0) & (np.hstack((0, dx)) >= 0))[0] ind = np.unique(np.hstack((ine, ire, ife))) - # handle NaN's - if ind.size and indnan.size: - # NaN's and values close to NaN's cannot be peaks - ind = ind[np.in1d(ind, np.unique(np.hstack((indnan, indnan - 1, indnan + 1))), invert=True)] - # first and last values of x cannot be peaks + + # Remove first and last values of x if they are detected as peaks if ind.size and ind[0] == 0: ind = ind[1:] if ind.size and ind[-1] == x.size - 1: ind = ind[:-1] - # remove peaks < minimum peak height + + print('{} raw peaks'.format(ind.size)) + + # Remove peaks < minimum peak height if ind.size and mph is not None: ind = ind[x[ind] >= mph] - # remove peaks - neighbors < threshold + print('{} height-filtered peaks'.format(ind.size)) + + # Remove peaks - neighbors < threshold if ind.size and threshold > 0: dx = np.min(np.vstack([x[ind] - x[ind - 1], x[ind] - x[ind + 1]]), axis=0) ind = np.delete(ind, np.where(dx < threshold)[0]) - # detect small peaks closer than minimum peak distance + print('{} prominence-filtered peaks'.format(ind.size)) + + # Detect small peaks closer than minimum peak distance if ind.size and mpd > 1: ind = ind[np.argsort(x[ind])][::-1] # sort ind by peak height idel = np.zeros(ind.size, dtype=bool) for i in range(ind.size): if not idel[i]: # keep peaks with the same height if kpsh is True idel = idel | (ind >= ind[i] - mpd) & (ind <= ind[i] + mpd) \ & (x[ind[i]] > x[ind] if kpsh else True) idel[i] = 0 # Keep current peak # remove the small peaks and sort back the indices by their occurrence ind = np.sort(ind[~idel]) + print('{} distance-filtered peaks'.format(ind.size)) return ind def detectPeaksTime(t, y, mph, mtd, mpp=0): """ Extension of the detectPeaks function to detect peaks in data based on their amplitude and time difference, with a non-uniform time vector. :param t: time vector (not necessarily uniform) :param y: signal :param mph: minimal peak height :param mtd: minimal time difference :mpp: minmal peak prominence :return: array of peak indexes """ # Determine whether time vector is uniform (threshold in time step variation) dt = np.diff(t) if (dt.max() - dt.min()) / dt.min() < 1e-2: isuniform = True else: isuniform = False if isuniform: - ipeaks = detectPeaks(y, mph, mpd=np.ceil(mtd / t[1] - t[0]), threshold=mpp) + print('uniform time vector') + dt = t[1] - t[0] + mpd = int(np.ceil(mtd / dt)) + ipeaks = detectPeaks(y, mph, mpd=mpd, threshold=mpp) else: + print('non-uniform time vector') # Detect peaks on signal with no restriction on inter-peak distance irawpeaks = detectPeaks(y, mph, mpd=1, threshold=mpp) npeaks = irawpeaks.size if npeaks > 0: # Filter relevant peaks with temporal distance ipeaks = [irawpeaks[0]] for i in range(1, npeaks): i1 = ipeaks[-1] i2 = irawpeaks[i] if t[i2] - t[i1] < mtd: if y[i2] > y[i1]: ipeaks[-1] = i2 else: ipeaks.append(i2) else: ipeaks = [] ipeaks = np.array(ipeaks) - if ipeaks.size > 0: - return ipeaks - else: - return None + return ipeaks def detectSpikes(t, Qm, min_amp, min_dt): ''' Detect spikes on a charge density signal, and return their number, latency and rate. :param t: time vector (s) :param Qm: charge density vector (C/m2) :param min_amp: minimal charge amplitude to detect spikes (C/m2) :param min_dt: minimal time interval between 2 spikes (s) :return: 3-tuple with number of spikes, latency (s) and spike rate (sp/s) ''' i_spikes = detectPeaksTime(t, Qm, min_amp, min_dt) - if i_spikes is not None: + if len(i_spikes) > 0: latency = t[i_spikes[0]] # s n_spikes = i_spikes.size if n_spikes > 1: first_to_last_spike = t[i_spikes[-1]] - t[i_spikes[0]] # s spike_rate = (n_spikes - 1) / first_to_last_spike # spikes/s else: spike_rate = 'N/A' else: latency = 'N/A' spike_rate = 'N/A' n_spikes = 0 return (n_spikes, latency, spike_rate) +def findPeaks(y, mph=None, mpd=None, mpp=None): + ''' Detect peaks in a signal based on their height, prominence and/or separating distance. + + :param y: signal vector + :param mph: minimum peak height (in signal units, default = None). + :param mpd: minimum inter-peak distance (in indexes, default = None) + :param mpp: minimum peak prominence (in signal units, default = None) + :return: 4-tuple of arrays with the indexes of peaks occurence, peaks prominence, + peaks width at half-prominence and peaks half-prominence bounds (left and right) + + Adapted from: + - Marco Duarte's detect_peaks function + (http://nbviewer.jupyter.org/github/demotu/BMC/blob/master/notebooks/DetectPeaks.ipynb) + - MATLAB findpeaks function (https://ch.mathworks.com/help/signal/ref/findpeaks.html) + ''' + + # Find all peaks and valleys + dy = np.diff(y) + s = np.sign(dy) + ipeaks = np.where(np.diff(s) < 0)[0] + 1 + ivalleys = np.where(np.diff(s) > 0)[0] + 1 + + if ipeaks.size == 0: + return (None,) * 4 + + # Ensure that peaks and valleys are coherently distributed + assert abs(ipeaks.size - ivalleys.size) <= 1, 'Number of peaks and valleys not matching' + njoints = min(ipeaks.size, ivalleys.size) + if ipeaks[0] < ivalleys[0]: + assert np.all(ivalleys[:njoints] > ipeaks[:njoints]), 'Peaks and valleys not interpersed' + else: + assert np.all(ivalleys[:njoints] < ipeaks[:njoints]), 'Peaks and valleys not interpersed' + + # Remove peaks < minimum peak height + if mph is not None: + ipeaks = ipeaks[y[ipeaks] >= mph] + if ipeaks.size == 0: + return (None,) * 4 + + # Detect small peaks closer than minimum peak distance + if mpd is not None: + ipeaks = ipeaks[np.argsort(y[ipeaks])][::-1] # sort ipeaks by descending peak height + idel = np.zeros(ipeaks.size, dtype=bool) # initialize boolean deletion array (all false) + for i in range(ipeaks.size): # for each peak + if not idel[i]: # if not marked for deletion + closepeaks = (ipeaks >= ipeaks[i] - mpd) & (ipeaks <= ipeaks[i] + mpd) # close peaks + idel = idel | closepeaks # mark for deletion along with previously marked peaks + # idel = idel | (ipeaks >= ipeaks[i] - mpd) & (ipeaks <= ipeaks[i] + mpd) + idel[i] = 0 # keep current peak + # remove the small peaks and sort back the indices by their occurrence + ipeaks = np.sort(ipeaks[~idel]) + + # Detect smallest valleys between consecutive relevant peaks + ibottomvalleys = [] + if ipeaks[0] > ivalleys[0]: + itrappedvalleys = ivalleys[ivalleys < ipeaks[0]] + ibottomvalleys.append(itrappedvalleys[np.argmin(y[itrappedvalleys])]) + for i, j in zip(ipeaks[:-1], ipeaks[1:]): + itrappedvalleys = ivalleys[np.logical_and(ivalleys > i, ivalleys < j)] + ibottomvalleys.append(itrappedvalleys[np.argmin(y[itrappedvalleys])]) + if ipeaks[-1] < ivalleys[-1]: + itrappedvalleys = ivalleys[ivalleys > ipeaks[-1]] + ibottomvalleys.append(itrappedvalleys[np.argmin(y[itrappedvalleys])]) + ipeaks = ipeaks + ivalleys = np.array(ibottomvalleys, dtype=int) + + # Ensure each peak is bounded by two valleys, adding signal boundaries as valleys if necessary + if ipeaks[0] < ivalleys[0]: + ivalleys = np.insert(ivalleys, 0, 0) + if ipeaks[-1] > ivalleys[-1]: + ivalleys = np.insert(ivalleys, ivalleys.size, y.size - 1) + + # Remove peaks < minimum peak prominence + if mpp is not None: + + # Compute peaks prominences as difference between peaks and their closest valley + prominences = y[ipeaks] - np.amax((y[ivalleys[:-1]], y[ivalleys[1:]]), axis=0) + + # initialize peaks and valleys deletion tables + idelp = np.zeros(ipeaks.size, dtype=bool) + idelv = np.zeros(ivalleys.size, dtype=bool) + + # for each peak (sorted by ascending prominence order) + for ind in np.argsort(prominences): + ipeak = ipeaks[ind] # get peak index + + # get peak bases as first valleys on either side not marked for deletion + indleftbase = ind + indrightbase = ind + 1 + while idelv[indleftbase]: + indleftbase -= 1 + while idelv[indrightbase]: + indrightbase += 1 + ileftbase = ivalleys[indleftbase] + irightbase = ivalleys[indrightbase] + + # Compute peak prominence and mark for deletion if < mpp + indmaxbase = indleftbase if y[ileftbase] > y[irightbase] else indrightbase + if y[ipeak] - y[ivalleys[indmaxbase]] < mpp: + idelp[ind] = True # mark peak for deletion + idelv[indmaxbase] = True # mark highest surrouding valley for deletion + + # remove irrelevant peaks and valleys, and sort back the indices by their occurrence + ipeaks = np.sort(ipeaks[~idelp]) + ivalleys = np.sort(ivalleys[~idelv]) + + if ipeaks.size == 0: + return (None,) * 4 + + # Compute peaks prominences and reference half-prominence levels + prominences = y[ipeaks] - np.amax((y[ivalleys[:-1]], y[ivalleys[1:]]), axis=0) + refheights = y[ipeaks] - prominences / 2 + + # Compute half-prominence bounds + ibounds = np.empty((ipeaks.size, 2)) + for i in range(ipeaks.size): + + # compute the index of the left-intercept at half max + ileft = ipeaks[i] + while ileft >= ivalleys[i] and y[ileft] > refheights[i]: + ileft -= 1 + # ileft = findLeftIntercept(y, ipeaks[i], ivalleys[i], refheights[i]) + if ileft < ivalleys[i]: # intercept exactly on valley + ibounds[i, 0] = ivalleys[i] + else: # interpolate intercept linearly between signal boundary points + a = (y[ileft + 1] - y[ileft]) / 1 + b = y[ileft] - a * ileft + ibounds[i, 0] = (refheights[i] - b) / a + + # compute the index of the right-intercept at half max + iright = ipeaks[i] + while iright <= ivalleys[i + 1] and y[iright] > refheights[i]: + iright += 1 + # iright = findRightIntercept(y, ipeaks[i], ivalleys[i + 1], refheights[i]) + if iright > ivalleys[i + 1]: # intercept exactly on valley + ibounds[i, 1] = ivalleys[i + 1] + else: # interpolate intercept linearly between signal boundary points + a = (y[iright + 1] - y[iright]) / 1 + b = y[iright] - a * iright + ibounds[i, 1] = (refheights[i] - b) / a + + # Compute peaks widths at half-prominence + widths = np.diff(ibounds, axis=1) + + return (ipeaks, prominences, widths, ibounds) + + def runMech(batch_dir, log_filepath, bls, Fdrive, Adrive, Qm): ''' Run a single simulation of the mechanical system with specific parameters and an imposed value of charge density, and save the results in a PKL file. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param bls: BilayerSonophore instance :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param Qm: applided membrane charge density (C/m2) :return: full path to the output file ''' simcode = MECH_code.format(bls.a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, Qm * 1e5) # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Run simulation tstart = time.time() (t, y, states) = bls.run(Fdrive, Adrive, Qm) (Z, ng) = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f seconds', tcomp) # Store dataframe and metadata df = pd.DataFrame({'t': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng}) meta = {'a': bls.a, 'd': bls.d, 'Cm0': bls.Cm0, 'Qm0': bls.Qm0, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'Qm': Qm, 'tcomp': tcomp} # Export into to PKL file output_filepath = '{}/{}.pkl'.format(batch_dir, simcode) with open(output_filepath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.debug('simulation data exported to "%s"', output_filepath) # Compute key output metrics Zmax = np.amax(Z) Zmin = np.amin(Z) Zabs_max = np.amax(np.abs([Zmin, Zmax])) eAmax = bls.arealstrain(Zabs_max) Tmax = bls.TEtot(Zabs_max) Pmmax = bls.PMavgpred(Zmin) ngmax = np.amax(ng) dUdtmax = np.amax(np.abs(np.diff(U) / np.diff(t)**2)) # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': bls.a * 1e9, 'D': bls.d * 1e6, 'E': Fdrive * 1e-3, 'F': Adrive * 1e-3, 'G': Qm * 1e5, 'H': t.size, 'I': tcomp, 'J': bls.kA + bls.kA_tissue, 'K': Zmax * 1e9, 'L': eAmax, 'M': Tmax * 1e3, 'N': (ngmax - bls.ng0) / bls.ng0, 'O': Pmmax * 1e-3, 'P': dUdtmax } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) return output_filepath def runMechBatch(batch_dir, log_filepath, Cm0, Qm0, stim_params, a=default_diam, d=default_embedding): ''' Run batch simulations of the mechanical system with imposed values of charge density, for various sonophore spans and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param Cm0: membrane resting capacitance (F/m2) :param Qm0: membrane resting charge density (C/m2) :param stim_params: dictionary containing sweeps for all stimulation parameters :param a: BLS in-plane diameter (m) :param d: depth of embedding tissue around plasma membrane (m) ''' # Checking validity of stimulation parameters mandatory_params = ['freqs', 'amps', 'charges'] for mp in mandatory_params: if mp not in stim_params: raise InputError('Missing stimulation parameter field: "{}"'.format(mp)) # Define logging format MECH_log = ('Mechanical simulation %u/%u (a = %.1f nm, d = %.1f um, f = %.2f kHz, ' 'A = %.2f kPa, Q = %.1f nC/cm2)') logger.info("Starting mechanical simulation batch") # Unpack stimulation parameters amps = np.array(stim_params['amps']) charges = np.array(stim_params['charges']) # Generate simulations queue nA = len(amps) nQ = len(charges) sim_queue = np.array(np.meshgrid(amps, charges)).T.reshape(nA * nQ, 2) nqueue = sim_queue.shape[0] # Run simulations nsims = len(stim_params['freqs']) * nqueue simcount = 0 filepaths = [] for Fdrive in stim_params['freqs']: # Create BilayerSonophore instance (modulus of embedding tissue depends on frequency) bls = BilayerSonophore(a, Fdrive, Cm0, Qm0, d) for i in range(nqueue): simcount += 1 Adrive, Qm = sim_queue[i, :] # Log logger.info(MECH_log, simcount, nsims, a * 1e9, d * 1e6, Fdrive * 1e-3, Adrive * 1e-3, Qm * 1e5) # Run simulation try: output_filepath = runMech(batch_dir, log_filepath, bls, Fdrive, Adrive, Qm) filepaths.append(output_filepath) except (Warning, AssertionError) as inst: logger.warning('Integration error: %s. Continue batch? (y/n)', extra={inst}) user_str = input() if user_str not in ['y', 'Y']: return filepaths return filepaths def runEStim(batch_dir, log_filepath, solver, neuron, Astim, tstim, toffset, PRF, DC): ''' Run a single E-STIM simulation a given neuron for specific stimulation parameters, and save the results in a PKL file. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param solver: SolverElec instance :param Astim: pulse amplitude (mA/m2) :param tstim: pulse duration (s) :param toffset: offset duration (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :return: full path to the output file ''' if DC == 1.0: simcode = ESTIM_CW_code.format(neuron.name, Astim, tstim * 1e3) else: simcode = ESTIM_PW_code.format(neuron.name, Astim, tstim * 1e3, PRF * 1e-3, DC) # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Run simulation tstart = time.time() (t, y, states) = solver.run(neuron, Astim, tstim, toffset, PRF, DC) Vm, *channels = y tcomp = time.time() - tstart logger.debug('completed in %.2f seconds', tcomp) # Store dataframe and metadata df = pd.DataFrame({'t': t, 'states': states, 'Vm': Vm}) for j in range(len(neuron.states_names)): df[neuron.states_names[j]] = channels[j] meta = {'neuron': neuron.name, 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'tcomp': tcomp} # Export into to PKL file output_filepath = '{}/{}.pkl'.format(batch_dir, simcode) with open(output_filepath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.debug('simulation data exported to "%s"', output_filepath) # Detect spikes on Vm signal n_spikes, lat, sr = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) logger.debug('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': Astim, 'E': tstim * 1e3, 'F': PRF * 1e-3 if DC < 1 else 'N/A', 'G': DC, 'H': t.size, 'I': round(tcomp, 2), 'J': n_spikes, 'K': lat * 1e3 if isinstance(lat, float) else 'N/A', 'L': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.debug('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) return output_filepath def titrateEStim(solver, neuron, Astim, tstim, toffset, PRF=1.5e3, DC=1.0): """ Use a dichotomic recursive search to determine the threshold value of a specific electric stimulation parameter needed to obtain neural excitation, keeping all other parameters fixed. The titration parameter can be stimulation amplitude, duration or any variable for which the number of spikes is a monotonically increasing function. This function is called recursively until an accurate threshold is found. :param solver: solver instance :param neuron: neuron object :param Astim: injected current density amplitude (mA/m2) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :return: 5-tuple with the determined amplitude threshold, time profile, solution matrix, state vector and response latency """ # Determine titration type if isinstance(Astim, tuple): t_type = 'A' interval = Astim thr = TITRATION_ESTIM_DA_MAX maxval = TITRATION_ESTIM_A_MAX elif isinstance(tstim, tuple): t_type = 't' interval = tstim thr = TITRATION_DT_THR maxval = TITRATION_T_MAX elif isinstance(DC, tuple): t_type = 'DC' interval = DC thr = TITRATION_DDC_THR maxval = TITRATION_DC_MAX else: logger.error('Invalid titration type') return 0. t_var = ESTIM_params[t_type] # Check amplitude interval and define current value if interval[0] >= interval[1]: raise InputError('Invaid {} interval: {} (must be defined as [lb, ub])' .format(t_type, interval)) value = (interval[0] + interval[1]) / 2 # Define stimulation parameters if t_type == 'A': stim_params = [value, tstim, toffset, PRF, DC] elif t_type == 't': stim_params = [Astim, value, toffset, PRF, DC] elif t_type == 'DC': stim_params = [Astim, tstim, toffset, PRF, value] # Run simulation and detect spikes (t, y, states) = solver.run(neuron, *stim_params) n_spikes, latency, _ = detectSpikes(t, y[0, :], SPIKE_MIN_VAMP, SPIKE_MIN_DT) logger.debug('%.2f %s ---> %u spike%s detected', value * t_var['factor'], t_var['unit'], n_spikes, "s" if n_spikes > 1 else "") # If accurate threshold is found, return simulation results if (interval[1] - interval[0]) <= thr and n_spikes == 1: return (value, t, y, states, latency) # Otherwise, refine titration interval and iterate recursively else: if n_spikes == 0: if (maxval - interval[1]) <= thr: # if upper bound too close to max then stop logger.warning('no spikes detected within titration interval') return (np.nan, t, y, states, latency) new_interval = (value, interval[1]) else: new_interval = (interval[0], value) stim_params[t_var['index']] = new_interval return titrateEStim(solver, neuron, *stim_params) def runEStimBatch(batch_dir, log_filepath, neurons, stim_params): ''' Run batch E-STIM simulations of the system for various neuron types and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :return: list of full paths to the output files ''' mandatory_params = ['amps', 'durations', 'offsets', 'PRFs', 'DCs'] for mp in mandatory_params: if mp not in stim_params: raise InputError('Missing stimulation parameter field: "{}"'.format(mp)) # Define logging format ESTIM_CW_log = 'E-STIM simulation %u/%u: %s neuron, A = %.1f mA/m2, t = %.1f ms' ESTIM_PW_log = ('E-STIM simulation %u/%u: %s neuron, A = %.1f mA/m2, t = %.1f ms, ' 'PRF = %.2f kHz, DC = %.2f') logger.info("Starting E-STIM simulation batch") # Generate simulations queue sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DCs']) nqueue = sim_queue.shape[0] # Initialize solver solver = SolverElec() # Run simulations nsims = len(neurons) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for i in range(nqueue): simcount += 1 Astim, tstim, toffset, PRF, DC = sim_queue[i, :] if DC == 1.0: logger.info(ESTIM_CW_log, simcount, nsims, neuron.name, Astim, tstim * 1e3) else: logger.info(ESTIM_PW_log, simcount, nsims, neuron.name, Astim, tstim * 1e3, PRF * 1e-3, DC) try: output_filepath = runEStim(batch_dir, log_filepath, solver, neuron, Astim, tstim, toffset, PRF, DC) filepaths.append(output_filepath) except (Warning, AssertionError) as inst: logger.warning('Integration error: %s. Continue batch? (y/n)', extra={inst}) user_str = input() if user_str not in ['y', 'Y']: return filepaths return filepaths def titrateEStimBatch(batch_dir, log_filepath, neurons, stim_params): ''' Run batch electrical titrations of the system for various neuron types and stimulation parameters, to determine the threshold of a specific stimulus parameter for neural excitation. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :return: list of full paths to the output files ''' # Define logging format ESTIM_titration_log = '%s neuron - E-STIM titration %u/%u (%s)' logger.info("Starting E-STIM titration batch") # Determine titration parameter and titrations list if 'durations' not in stim_params: t_type = 't' sim_queue = createSimQueue(stim_params['amps'], [None], [TITRATION_T_OFFSET], stim_params['PRFs'], stim_params['DCs']) elif 'amps' not in stim_params: t_type = 'A' sim_queue = createSimQueue([None], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], stim_params['DCs']) elif 'DC' not in stim_params: t_type = 'DC' sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], [None]) nqueue = sim_queue.shape[0] t_var = ESTIM_params[t_type] # Create SolverElec instance solver = SolverElec() # Run titrations nsims = len(neurons) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for i in range(nqueue): simcount += 1 # Extract parameters Astim, tstim, toffset, PRF, DC = sim_queue[i, :] if Astim is None: Astim = (0., 2 * TITRATION_ESTIM_A_MAX) elif tstim is None: tstim = (0., 2 * TITRATION_T_MAX) elif DC is None: DC = (0., 2 * TITRATION_DC_MAX) curr_params = [Astim, tstim, PRF, DC] # Generate log str log_str = '' pnames = list(ESTIM_params.keys()) j = 0 for cp in curr_params: pn = pnames[j] pi = ESTIM_params[pn] if not isinstance(cp, tuple): if log_str: log_str += ', ' log_str += '{} = {:.2f} {}'.format(pn, pi['factor'] * cp, pi['unit']) j += 1 # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log logger.info(ESTIM_titration_log, neuron.name, simcount, nsims, log_str) # Run titration tstart = time.time() try: (output_thr, t, y, states, lat) = titrateEStim(solver, neuron, Astim, tstim, toffset, PRF, DC) Vm, *channels = y tcomp = time.time() - tstart logger.info('completed in %.2f s, threshold = %.2f %s', tcomp, output_thr * t_var['factor'], t_var['unit']) # Determine output variable if t_type == 'A': Astim = output_thr elif t_type == 't': tstim = output_thr elif t_type == 'DC': DC = output_thr # Define output naming if DC == 1.0: simcode = ESTIM_CW_code.format(neuron.name, Astim, tstim * 1e3) else: simcode = ESTIM_PW_code.format(neuron.name, Astim, tstim * 1e3, PRF * 1e-3, DC) # Store dataframe and metadata df = pd.DataFrame({'t': t, 'states': states, 'Vm': Vm}) for j in range(len(neuron.states_names)): df[neuron.states_names[j]] = channels[j] meta = {'neuron': neuron.name, 'Astim': Astim, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'tcomp': tcomp} # Export into to PKL file output_filepath = '{}/{}.pkl'.format(batch_dir, simcode) with open(output_filepath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.info('simulation data exported to "%s"', output_filepath) filepaths.append(output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': Astim, 'E': tstim * 1e3, 'F': PRF * 1e-3 if DC < 1 else 'N/A', 'G': DC, 'H': t.size, 'I': round(tcomp, 2), 'J': n_spikes, 'K': lat * 1e3 if isinstance(lat, float) else 'N/A', 'L': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except (Warning, AssertionError) as inst: logger.warning('Integration error: %s. Continue batch? (y/n)', extra={inst}) user_str = input() if user_str not in ['y', 'Y']: return filepaths return filepaths def runAStim(batch_dir, log_filepath, solver, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, int_method='effective'): ''' Run a single A-STIM simulation a given neuron for specific stimulation parameters, and save the results in a PKL file. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param solver: SolverUS instance :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param int_method: selected integration method :return: full path to the output file ''' if DC == 1.0: simcode = ASTIM_CW_code.format(neuron.name, solver.a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, int_method) else: simcode = ASTIM_PW_code.format(neuron.name, solver.a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DC, int_method) # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Run simulation tstart = time.time() (t, y, states) = solver.run(neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, int_method) Z, ng, Qm, *channels = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.debug('completed in %.2f seconds', tcomp) # Store dataframe and metadata df = pd.DataFrame({'t': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng, 'Qm': Qm, 'Vm': Qm * 1e3 / np.array([solver.Capct(ZZ) for ZZ in Z])}) for j in range(len(neuron.states_names)): df[neuron.states_names[j]] = channels[j] meta = {'neuron': neuron.name, 'a': solver.a, 'd': solver.d, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'tcomp': tcomp} # Export into to PKL file output_filepath = '{}/{}.pkl'.format(batch_dir, simcode) with open(output_filepath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.debug('simulation data exported to "%s"', output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.debug('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': solver.a * 1e9, 'E': solver.d * 1e6, 'F': Fdrive * 1e-3, 'G': Adrive * 1e-3, 'H': tstim * 1e3, 'I': PRF * 1e-3 if DC < 1 else 'N/A', 'J': DC, 'K': int_method, 'L': t.size, 'M': round(tcomp, 2), 'N': n_spikes, 'O': lat * 1e3 if isinstance(lat, float) else 'N/A', 'P': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.debug('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) return output_filepath def titrateAStim(solver, neuron, Fdrive, Adrive, tstim, toffset, PRF=1.5e3, DC=1.0, int_method='effective'): """ Use a dichotomic recursive search to determine the threshold value of a specific acoustic stimulation parameter needed to obtain neural excitation, keeping all other parameters fixed. The titration parameter can be stimulation amplitude, duration or any variable for which the number of spikes is a monotonically increasing function. This function is called recursively until an accurate threshold is found. :param solver: solver instance :param neuron: neuron object :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param tstim: duration of US stimulation (s) :param toffset: duration of the offset (s) :param PRF: pulse repetition frequency (Hz) :param DC: pulse duty cycle (-) :param int_method: selected integration method :return: 5-tuple with the determined amplitude threshold, time profile, solution matrix, state vector and response latency """ # Determine titration type if isinstance(Adrive, tuple): t_type = 'A' interval = Adrive thr = TITRATION_ASTIM_DA_MAX maxval = TITRATION_ASTIM_A_MAX elif isinstance(tstim, tuple): t_type = 't' interval = tstim thr = TITRATION_DT_THR maxval = TITRATION_T_MAX elif isinstance(DC, tuple): t_type = 'DC' interval = DC thr = TITRATION_DDC_THR maxval = TITRATION_DC_MAX else: logger.error('Invalid titration type') return 0. t_var = ASTIM_params[t_type] # Check amplitude interval and define current value if interval[0] >= interval[1]: raise InputError('Invaid {} interval: {} (must be defined as [lb, ub])' .format(t_type, interval)) value = (interval[0] + interval[1]) / 2 # Define stimulation parameters if t_type == 'A': stim_params = [Fdrive, value, tstim, toffset, PRF, DC] elif t_type == 't': stim_params = [Fdrive, Adrive, value, toffset, PRF, DC] elif t_type == 'DC': stim_params = [Fdrive, Adrive, tstim, toffset, PRF, value] # Run simulation and detect spikes (t, y, states) = solver.run(neuron, *stim_params, int_method) n_spikes, latency, _ = detectSpikes(t, y[2, :], SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.debug('%.2f %s ---> %u spike%s detected', value * t_var['factor'], t_var['unit'], n_spikes, "s" if n_spikes > 1 else "") # If accurate threshold is found, return simulation results if (interval[1] - interval[0]) <= thr and n_spikes == 1: return (value, t, y, states, latency) # Otherwise, refine titration interval and iterate recursively else: if n_spikes == 0: if (maxval - interval[1]) <= thr: # if upper bound too close to max then stop logger.warning('no spikes detected within titration interval') return (np.nan, t, y, states, latency) new_interval = (value, interval[1]) else: new_interval = (interval[0], value) stim_params[t_var['index']] = new_interval return titrateAStim(solver, neuron, *stim_params, int_method) def runAStimBatch(batch_dir, log_filepath, neurons, stim_params, a=default_diam, int_method='effective'): ''' Run batch simulations of the system for various neuron types, sonophore and stimulation parameters. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :param a: BLS structure diameter (m) :param int_method: selected integration method :return: list of full paths to the output files ''' mandatory_params = ['freqs', 'amps', 'durations', 'offsets', 'PRFs', 'DCs'] for mp in mandatory_params: if mp not in stim_params: raise InputError('Missing stimulation parameter field: "{}"'.format(mp)) # Define logging format ASTIM_CW_log = ('A-STIM %s simulation %u/%u: %s neuron, a = %.1f nm, f = %.2f kHz, ' 'A = %.2f kPa, t = %.2f ms') ASTIM_PW_log = ('A-STIM %s simulation %u/%u: %s neuron, a = %.1f nm, f = %.2f kHz, ' 'A = %.2f kPa, t = %.2f ms, PRF = %.2f kHz, DC = %.3f') logger.info("Starting A-STIM simulation batch") # Generate simulations queue sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DCs']) nqueue = sim_queue.shape[0] # Run simulations nsims = len(neurons) * len(stim_params['freqs']) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for Fdrive in stim_params['freqs']: # Initialize SolverUS solver = SolverUS(a, neuron, Fdrive) for i in range(nqueue): simcount += 1 Adrive, tstim, toffset, PRF, DC = sim_queue[i, :] # Log and define naming if DC == 1.0: logger.info(ASTIM_CW_log, int_method, simcount, nsims, neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3) else: logger.info(ASTIM_PW_log, int_method, simcount, nsims, neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DC) # Run simulation try: output_filepath = runAStim(batch_dir, log_filepath, solver, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC, int_method) filepaths.append(output_filepath) except (Warning, AssertionError) as inst: logger.warning('Integration error: %s. Continue batch? (y/n)', extra={inst}) user_str = input() if user_str not in ['y', 'Y']: return filepaths return filepaths def titrateAStimBatch(batch_dir, log_filepath, neurons, stim_params, a=default_diam, int_method='effective'): ''' Run batch acoustic titrations of the system for various neuron types, sonophore and stimulation parameters, to determine the threshold of a specific stimulus parameter for neural excitation. :param batch_dir: full path to output directory of batch :param log_filepath: full path log file of batch :param neurons: list of neurons names :param stim_params: dictionary containing sweeps for all stimulation parameters :param a: BLS structure diameter (m) :param int_method: selected integration method :return: list of full paths to the output files ''' # Define logging format ASTIM_titration_log = '%s neuron - A-STIM titration %u/%u (a = %.1f nm, %s)' logger.info("Starting A-STIM titration batch") # Define default parameters int_method = 'effective' # Determine titration parameter and titrations list if 'durations' not in stim_params: t_type = 't' sim_queue = createSimQueue(stim_params['amps'], [None], [TITRATION_T_OFFSET], stim_params['PRFs'], stim_params['DCs']) elif 'amps' not in stim_params: t_type = 'A' sim_queue = createSimQueue([None], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], stim_params['DCs']) elif 'DC' not in stim_params: t_type = 'DC' sim_queue = createSimQueue(stim_params['amps'], stim_params['durations'], [TITRATION_T_OFFSET] * len(stim_params['durations']), stim_params['PRFs'], [None]) nqueue = sim_queue.shape[0] t_var = ASTIM_params[t_type] # Run titrations nsims = len(neurons) * len(stim_params['freqs']) * nqueue simcount = 0 filepaths = [] for nname in neurons: neuron = getNeuronsDict()[nname]() for Fdrive in stim_params['freqs']: # Create SolverUS instance (modulus of embedding tissue depends on frequency) solver = SolverUS(a, neuron, Fdrive) for i in range(nqueue): simcount += 1 # Extract parameters Adrive, tstim, toffset, PRF, DC = sim_queue[i, :] if Adrive is None: Adrive = (0., 2 * TITRATION_ASTIM_A_MAX) elif tstim is None: tstim = (0., 2 * TITRATION_T_MAX) elif DC is None: DC = (0., 2 * TITRATION_DC_MAX) curr_params = [Fdrive, Adrive, tstim, PRF, DC] # Generate log str log_str = '' pnames = list(ASTIM_params.keys()) j = 0 for cp in curr_params: pn = pnames[j] pi = ASTIM_params[pn] if not isinstance(cp, tuple): if log_str: log_str += ', ' log_str += '{} = {:.2f} {}'.format(pn, pi['factor'] * cp, pi['unit']) j += 1 # Get date and time info date_str = time.strftime("%Y.%m.%d") daytime_str = time.strftime("%H:%M:%S") # Log logger.info(ASTIM_titration_log, neuron.name, simcount, nsims, a * 1e9, log_str) # Run titration tstart = time.time() try: (output_thr, t, y, states, lat) = titrateAStim(solver, neuron, Fdrive, Adrive, tstim, toffset, PRF, DC) Z, ng, Qm, *channels = y U = np.insert(np.diff(Z) / np.diff(t), 0, 0.0) tcomp = time.time() - tstart logger.info('completed in %.2f s, threshold = %.2f %s', tcomp, output_thr * t_var['factor'], t_var['unit']) # Determine output variable if t_type == 'A': Adrive = output_thr elif t_type == 't': tstim = output_thr elif t_type == 'DC': DC = output_thr # Define output naming if DC == 1.0: simcode = ASTIM_CW_code.format(neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, int_method) else: simcode = ASTIM_PW_code.format(neuron.name, a * 1e9, Fdrive * 1e-3, Adrive * 1e-3, tstim * 1e3, PRF * 1e-3, DC, int_method) # Store dataframe and metadata df = pd.DataFrame({'t': t, 'states': states, 'U': U, 'Z': Z, 'ng': ng, 'Qm': Qm, 'Vm': Qm * 1e3 / np.array([solver.Capct(ZZ) for ZZ in Z])}) for j in range(len(neuron.states_names)): df[neuron.states_names[j]] = channels[j] meta = {'neuron': neuron.name, 'a': solver.a, 'd': solver.d, 'Fdrive': Fdrive, 'Adrive': Adrive, 'phi': np.pi, 'tstim': tstim, 'toffset': toffset, 'PRF': PRF, 'DC': DC, 'tcomp': tcomp} # Export into to PKL file output_filepath = '{}/{}.pkl'.format(batch_dir, simcode) with open(output_filepath, 'wb') as fh: pickle.dump({'meta': meta, 'data': df}, fh) logger.debug('simulation data exported to "%s"', output_filepath) filepaths.append(output_filepath) # Detect spikes on Qm signal n_spikes, lat, sr = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) logger.info('%u spike%s detected', n_spikes, "s" if n_spikes > 1 else "") # Export key metrics to log file log = { 'A': date_str, 'B': daytime_str, 'C': neuron.name, 'D': solver.a * 1e9, 'E': solver.d * 1e6, 'F': Fdrive * 1e-3, 'G': Adrive * 1e-3, 'H': tstim * 1e3, 'I': PRF * 1e-3 if DC < 1 else 'N/A', 'J': DC, 'K': int_method, 'L': t.size, 'M': round(tcomp, 2), 'N': n_spikes, 'O': lat * 1e3 if isinstance(lat, float) else 'N/A', 'P': sr * 1e-3 if isinstance(sr, float) else 'N/A' } if xlslog(log_filepath, 'Data', log) == 1: logger.info('log exported to "%s"', log_filepath) else: logger.error('log export to "%s" aborted', log_filepath) except (Warning, AssertionError) as inst: logger.warning('Integration error: %s. Continue batch? (y/n)', extra={inst}) user_str = input() if user_str not in ['y', 'Y']: return filepaths return filepaths