diff --git a/PointNICE/plt/pltutils.py b/PointNICE/plt/pltutils.py index 279a0c8..3db09d6 100644 --- a/PointNICE/plt/pltutils.py +++ b/PointNICE/plt/pltutils.py @@ -1,851 +1,851 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-23 14:55:37 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-03-13 00:44:39 +# @Last Modified time: 2018-03-13 12:27:43 ''' 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.pyplot as plt from matplotlib.patches import Rectangle import matplotlib.cm as cm import pandas as pd from .. import neurons from ..utils import getNeuronsDict, getLookupDir, rescale from ..bls import BilayerSonophore from .pltvars import pltvars # Get package logger logger = logging.getLogger('PointNICE') # Define global variables neuron = None bls = None # 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 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(yvars, filepaths, labels=None, fs=15, show_patches=True): ''' Compare profiles of several specific output variables of NICE simulations. :param yvars: list of variables names 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 show_patches: boolean indicating whether to indicate periods of stimulation with colored rectangular patches ''' # check labels if given if labels: assert len(labels) == len(filepaths), 'labels do not match number of compared files' assert all(isinstance(x, str) for x in labels), 'labels must be string typed' nvars = len(yvars) # y variables plotting information for key in yvars: assert key in pltvars, 'unknown plot variable "{}"'.format(key) y_pltvars = [pltvars[key] for key in yvars] # Dictionary of neurons neurons_dict = getNeuronsDict() # Initialize figure and axes if nvars == 1: _, ax = plt.subplots(figsize=(11, 4)) axes = [ax] else: _, axes = plt.subplots(nvars, 1, figsize=(11, min(3 * nvars, 9))) for i in range(nvars): ax = axes[i] pltvar = y_pltvars[i] if 'min' in pltvar and 'max' in pltvar: ax.set_ylim(pltvar['min'], pltvar['max']) if pltvar['unit']: ax.set_ylabel('${}\ ({})$'.format(pltvar['label'], pltvar['unit']), fontsize=fs) else: ax.set_ylabel('${}$'.format(pltvar['label']), fontsize=fs) if i < nvars - 1: ax.get_xaxis().set_ticklabels([]) else: for tick in ax.xaxis.get_major_ticks(): tick.label.set_fontsize(fs) ax.locator_params(axis='y', nbins=2) for tick in ax.yaxis.get_major_ticks(): tick.label.set_fontsize(fs) # Loop through data files j = 0 aliases = {} for filepath in filepaths: pkl_filename = ntpath.basename(filepath) # 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) assert sim_type in ['MECH', 'ASTIM', 'ESTIM'], 'invalid stimulation type' if j == 0: sim_type_ref = sim_type else: assert sim_type == sim_type_ref, 'Error: comparing 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 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] + 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) # Extract variables to plot vrs = [] for i in range(nvars): pltvar = y_pltvars[i] 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[yvars[i]].values if var.size == t.size - 1: var = np.insert(var, 0, var[0]) vrs.append(var) # 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, Fdrive * 1e-3, meta['Adrive'] * 1e-3, meta['tstim'] * 1e3) else: label = 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': label = MECH_title.format(a * 1e9, Fdrive * 1e-3, meta['Adrive'] * 1e-3) # Plotting handles = [axes[i].plot(t * t_plt['factor'], vrs[i] * y_pltvars[i]['factor'], linewidth=2, label=label) for i in range(nvars)] if show_patches: k = 0 # stimulation patches for ax in axes: handle = handles[k] (ybottom, ytop) = ax.get_ylim() la = [] 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=handle[0].get_color(), alpha=0.2))) aliases[handle[0]] = la k += 1 j += 1 # set x-axis label axes[-1].set_xlabel('${}\ ({})$'.format(t_plt['label'], t_plt['unit']), fontsize=fs) plt.tight_layout() iLegends = [] for k in range(nvars): axes[k].legend(loc=7, fontsize=fs) iLegends.append(InteractiveLegend(axes[k].legend_, aliases)) plt.show() 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: assert key in pltvars, '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) assert sim_type in ['MECH', 'ASTIM', 'ESTIM'], 'invalid stimulation 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] + 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) assert os.path.isfile(lookup_path), ('No lookup file available for {} ' 'neuron type').format(neuron.name) # 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) assert frange[0] <= Fdrive <= frange[1], \ 'Fdrive must be within [{:.1f}, {:.1f}] kHz'.format(*[f * 1e-3 for f in frange]) # 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/plt/pltvars.py b/PointNICE/plt/pltvars.py index 9744db4..140b5f7 100644 --- a/PointNICE/plt/pltvars.py +++ b/PointNICE/plt/pltvars.py @@ -1,572 +1,572 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-21 14:33:36 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-12-13 11:26:44 +# @Last Modified time: 2018-03-13 12:27:12 ''' Dictionary of plotting settings for output variables of the model. ''' pltvars = { 't_ms': { 'desc': 'time', 'label': 'time', 'unit': 'ms', 'factor': 1e3, 'onset': 3e-3 }, 't_us': { 'desc': 'time', 'label': 'time', 'unit': 'us', 'factor': 1e6, 'onset': 1e-6 }, 'Z': { 'desc': 'leaflets deflection', 'label': 'Z', 'unit': 'nm', 'factor': 1e9, 'min': -1.0, 'max': 10.0 }, 'ng': { 'desc': 'gas content', 'label': 'gas', 'unit': '10^{-22}\ mol', 'factor': 1e22, 'min': 1.0, 'max': 15.0 }, 'Pac': { 'desc': 'acoustic pressure', 'label': 'P_{AC}', 'unit': 'kPa', 'factor': 1e-3, - 'alias': 'bls.Pacoustic(t, data["Adrive"] * states, Fdrive)' + 'alias': 'bls.Pacoustic(t, meta["Adrive"] * states, Fdrive)' }, 'Pmavg': { 'desc': 'average intermolecular pressure', 'label': 'P_M', 'unit': 'kPa', 'factor': 1e-3, - 'alias': 'bls.PMavgpred(data["Z"])' + 'alias': 'bls.PMavgpred(df["Z"].values)' }, 'Telastic': { 'desc': 'leaflet elastic tension', 'label': 'T_E', 'unit': 'mN/m', 'factor': 1e3, - 'alias': 'bls.TEleaflet(data["Z"])' + 'alias': 'bls.TEleaflet(df["Z"].values)' }, 'Qm': { 'desc': 'charge density', 'label': 'Q_m', 'unit': 'nC/cm^2', 'factor': 1e5, 'min': -100, 'max': 50 }, 'Cm': { 'desc': 'membrane capacitance', 'label': 'C_m', 'unit': 'uF/cm^2', 'factor': 1e2, 'min': 0.0, 'max': 1.5, - 'alias': 'np.array([bls.Capct(ZZ) for ZZ in data["Z"]])' + 'alias': 'np.array([bls.Capct(ZZ) for ZZ in df["Z"].values])' }, 'Vm': { 'desc': 'membrane potential', 'label': 'V_m', 'unit': 'mV', 'factor': 1, }, 'm': { 'desc': 'iNa activation gate opening', 'label': 'm-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'h': { 'desc': 'iNa inactivation gate opening', 'label': 'h-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'm2h': { 'desc': 'iNa relative conductance', 'label': 'm^2h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'data["m"]**2 * data["h"]' + 'alias': 'df["m"].values**2 * df["h"].values' }, 'm3h': { 'desc': 'iNa relative conductance', 'label': 'm^3h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'data["m"]**3 * data["h"]' + 'alias': 'df["m"].values**3 * df["h"].values' }, 'm4h': { 'desc': 'iNa relative conductance', 'label': 'm^4h', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'data["m"]**4 * data["h"]' + 'alias': 'df["m"].values**4 * df["h"].values' }, 'n': { 'desc': 'iK activation gate opening', 'label': 'n-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'p': { 'desc': 'iM activation gate opening', 'label': 'p-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 's': { 'desc': 'iCa activation gates opening', 'label': 's-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'u': { 'desc': 'iCa inactivation gates opening', 'label': 'u-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 's2u': { 'desc': 'iT relative conductance', 'label': 's^2u', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'data["s"]**2 * data["u"]' + 'alias': 'df["s"].values**2 * df["u"].values' }, 'c': { 'desc': 'iKCa activation gates opening', 'label': 'c-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'a': { 'desc': 'iKA activation gates opening', 'label': 'a-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'b': { 'desc': 'iKA inactivation gates opening', 'label': 'b-gate', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'ab': { 'desc': 'iKA relative conductance', 'label': 'ab', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': 'data["a"] * data["b"]' + 'alias': 'df["a"].values * df["b"].values' }, 'O': { 'desc': 'iH activation gate opening', 'label': 'O', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'OL': { 'desc': 'iH activation gate locked-opening', 'label': 'O_L', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1, - 'alias': '1 - data["O"] - data["C"]' + 'alias': '1 - df["O"].values - df["C"].values' }, 'O + 2OL': { 'desc': 'iH activation gate relative conductance', 'label': 'O\ +\ 2O_L', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 2.1, - 'alias': 'data["O"] + 2 * (1 - data["O"] - data["C"])' + 'alias': 'df["O"].values + 2 * (1 - df["O"].values - df["C"].values)' }, 'P0': { 'desc': 'iH regulating factor activation', 'label': 'P_0', 'unit': None, 'factor': 1, 'min': -0.1, 'max': 1.1 }, 'C_Ca': { 'desc': 'sumbmembrane Ca2+ concentration', 'label': '[Ca^{2+}]_i', 'unit': 'uM', 'factor': 1e6 # 'min': 0, # 'max': 150.0 }, 'C_Na': { 'desc': 'sumbmembrane Na+ concentration', 'label': '[Na^+]_i', 'unit': 'uM', 'factor': 1e6 # 'min': 0, # 'max': 150.0 }, 'C_Na_arb': { 'key': 'C_Na', 'desc': 'submembrane Na+ concentration', 'label': '[Na^+]', 'unit': 'arb.', 'factor': 1 }, 'C_Na_arb_activation': { 'key': 'A_Na', 'desc': 'Na+ dependent PumpNa current activation', 'label': 'A_{Na^+}', 'unit': 'arb', 'factor': 1 }, 'C_Ca_arb': { 'key': 'C_Ca', 'desc': 'submembrane Ca2+ concentration', 'label': '[Ca^{2+}]', 'unit': 'arb.', 'factor': 1 }, 'C_Ca_arb_activation': { 'key': 'A_Ca', 'desc': 'Ca2+ dependent Potassium current activation', 'label': 'A_{Ca^{2+}}', 'unit': 'arb', 'factor': 1 }, 'VL': { 'constant': 'neuron.VL', 'desc': 'non-specific leakage current resting potential', 'label': 'V_L', 'unit': 'mV', 'factor': 1e0 }, 'iL': { 'desc': 'leakage current', 'label': 'I_L', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currL(data["Vm"])' + 'alias': 'neuron.currL(df["Vm"].values)' }, 'iNa': { 'desc': 'Sodium current', 'label': 'I_{Na}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currNa(data["m"], data["h"], data["Vm"])' + 'alias': 'neuron.currNa(df["m"].values, df["h"].values, df["Vm"].values)' }, 'iNa2': { 'desc': 'Sodium current', 'label': 'I_{Na}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currNa(data["m"], data["h"], data["Vm"], data["C_Na"])' + 'alias': 'neuron.currNa(df["m"].values, df["h"].values, df["Vm"].values, df["C_Na"].values)' }, 'iK': { 'desc': 'delayed-rectifier Potassium current', 'label': 'I_K', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currK(data["n"], data["Vm"])' + 'alias': 'neuron.currK(df["n"].values, df["Vm"].values)' }, 'iM': { 'desc': 'slow non-inactivating Potassium current', 'label': 'I_M', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currM(data["p"], data["Vm"])' + 'alias': 'neuron.currM(df["p"].values, df["Vm"].values)' }, 'iA': { 'desc': 'transient Potassium current', 'label': 'I_A', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currA(data["a"], data["b"], data["Vm"])' + 'alias': 'neuron.currA(df["a"].values, df["b"].values, df["Vm"].values)' }, 'iT': { 'desc': 'low-threshold Calcium current', 'label': 'I_T', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currCa(data["s"], data["u"], data["Vm"])' + 'alias': 'neuron.currCa(df["s"].values, df["u"].values, df["Vm"].values)' }, 'iTs': { 'desc': 'low-threshold Calcium current', 'label': 'I_{TS}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currCa(data["s"], data["u"], data["Vm"])' + 'alias': 'neuron.currCa(df["s"].values, df["u"].values, df["Vm"].values)' }, 'iCa': { 'desc': 'leech Calcium current', 'label': 'I_{Ca}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currCa(data["s"], data["Vm"])' + 'alias': 'neuron.currCa(df["s"].values, df["Vm"].values)' }, 'iCa2': { 'desc': 'leech Calcium current', 'label': 'I_{Ca}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currCa(data["s"], data["Vm"], data["C_Ca"])' + 'alias': 'neuron.currCa(df["s"].values, df["Vm"].values, df["C_Ca"].values)' }, 'iH': { 'desc': 'hyperpolarization-activated cationic current', 'label': 'I_h', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currH(data["O"], data["C"], data["Vm"])' + 'alias': 'neuron.currH(df["O"].values, df["C"].values, df["Vm"].values)' }, 'iKL': { 'desc': 'leakage Potassium current', 'label': 'I_{KL}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currKL(data["Vm"])' + 'alias': 'neuron.currKL(df["Vm"].values)' }, 'iKCa': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{KCa}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currKCa(data["A_Ca"], data["Vm"])' + 'alias': 'neuron.currKCa(df["A_Ca"].values, df["Vm"].values)' }, 'iKCa2': { 'desc': 'Calcium-activated Potassium current', 'label': 'I_{KCa}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currKCa(data["c"], data["Vm"])' + 'alias': 'neuron.currKCa(df["c"].values, df["Vm"].values)' }, 'iPumpNa': { 'desc': 'Outward current mimicking the activity of the NaK-ATPase pump', 'label': 'I_{PumpNa}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currPumpNa(data["A_Na"], data["Vm"])' + 'alias': 'neuron.currPumpNa(df["A_Na"].values, df["Vm"].values)' }, 'iPumpNa2': { 'desc': 'Outward current mimicking the activity of the NaK-ATPase pump', 'label': 'I_{PumpNa}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currPumpNa(data["C_Na"])' + 'alias': 'neuron.currPumpNa(df["C_Na"].values)' }, 'iPumpCa2': { 'desc': 'Outward current describing the removal of Ca2+ from the intracellular space', 'label': 'I_{PumpCa}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currPumpCa(data["C_Ca"])' + 'alias': 'neuron.currPumpCa(df["C_Ca"].values)' }, 'iNet': { 'desc': 'net current', 'label': 'I_{net}', 'unit': 'A/m^2', 'factor': 1e-3, - 'alias': 'neuron.currNet(data["Vm"], neuron_states)' + 'alias': 'neuron.currNet(df["Vm"].values, neuron_states)' }, 'Veff': { 'key': 'V', 'desc': 'effective membrane potential', 'label': 'V_{m, eff}', 'unit': 'mV', 'factor': 1e0 }, 'alpham': { 'desc': 'iNa m-gate activation rate', 'label': '\\alpha_{m,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betam': { 'desc': 'iNa m-gate inactivation rate', 'label': '\\beta_{m,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphah': { 'desc': 'iNa h-gate activation rate', 'label': '\\alpha_{h,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betah': { 'desc': 'iNa h-gate inactivation rate', 'label': '\\beta_{h,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphan': { 'desc': 'iK n-gate activation rate', 'label': '\\alpha_{n,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betan': { 'desc': 'iK n-gate inactivation rate', 'label': '\\beta_{n,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphap': { 'desc': 'iM p-gate activation rate', 'label': '\\alpha_{p,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betap': { 'desc': 'iM p-gate inactivation rate', 'label': '\\beta_{p,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphas': { 'desc': 'iT s-gate activation rate', 'label': '\\alpha_{s,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betas': { 'desc': 'iT s-gate inactivation rate', 'label': '\\beta_{s,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphau': { 'desc': 'iT u-gate activation rate', 'label': '\\alpha_{u,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betau': { 'desc': 'iT u-gate inactivation rate', 'label': '\\beta_{u,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'alphao': { 'desc': 'iH channels activation rate (between closed and open forms)', 'label': '\\alpha_{O,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 }, 'betao': { 'desc': 'iH channels inactivation rate (between closed and open forms)', 'label': '\\beta_{O,\ eff}', 'unit': 'ms^-1', 'factor': 1e-3 } } diff --git a/postpro/postpro_rmse_charge.py b/postpro/postpro_rmse_charge.py index 68e54f8..0938889 100644 --- a/postpro/postpro_rmse_charge.py +++ b/postpro/postpro_rmse_charge.py @@ -1,73 +1,77 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-01 16:35:43 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-09-10 17:20:19 +# @Last Modified time: 2018-03-13 12:35:30 """ Compute RMSE between charge profiles of NICE output. """ import sys import pickle import ntpath import numpy as np from PointNICE.utils import OpenFilesDialog, rmse # Define options pkl_root = "../../Output/test Elec/" t_offset = 10e-3 # s # Select data files (PKL) pkl_filepaths, pkl_dir = OpenFilesDialog('pkl') # Quit if no file selected if not pkl_filepaths: print('error: no input file') sys.exit(1) # Quit if more than 2 files if len(pkl_filepaths) > 2: print('error: cannot compare more than 2 methods') sys.exit(1) # Load data from file 1 pkl_filename = ntpath.basename(pkl_filepaths[0]) print('Loading data from "' + pkl_filename + '"') -with open(pkl_filepaths[0], 'rb') as pkl_file: - data = pickle.load(pkl_file) +with open(pkl_filepaths[0], 'rb') as fh: + frame = pickle.load(fh) + df = frame['data'] + meta = frame['meta'] -t1 = data['t'] -tstim1 = data['tstim'] -toffset1 = data['toffset'] -f1 = data['Fdrive'] -A1 = data['Adrive'] -Q1 = data['Qm'] * 1e2 # nC/cm2 -states1 = data['states'] +tstim1 = meta['tstim'] +toffset1 = meta['toffset'] +f1 = meta['Fdrive'] +A1 = meta['Adrive'] +t1 = df['t'].values +Q1 = df['Qm'].values * 1e2 # nC/cm2 +states1 = df['states'].values # Load data from file 2 pkl_filename = ntpath.basename(pkl_filepaths[1]) print('Loading data from "' + pkl_filename + '"') -with open(pkl_filepaths[1], 'rb') as pkl_file: - data = pickle.load(pkl_file) +with open(pkl_filepaths[1], 'rb') as fh: + frame = pickle.load(fh) + df = frame['data'] + meta = frame['meta'] -t2 = data['t'] -tstim2 = data['tstim'] -toffset2 = data['toffset'] -f2 = data['Fdrive'] -A2 = data['Adrive'] -Q2 = data['Qm'] * 1e2 # nC/cm2 -states2 = data['states'] +tstim2 = meta['tstim'] +toffset2 = meta['toffset'] +f2 = meta['Fdrive'] +A2 = meta['Adrive'] +t2 = df['t'].values +Q2 = df['Qm'].values * 1e2 # nC/cm2 +states2 = df['states'].values if tstim1 != tstim2 or f1 != f2 or A1 != A2 or toffset1 != toffset2: print('error: different stimulation conditions') else: print('comparing charge profiles') tcomp = np.arange(0, tstim1 + toffset1, 1e-3) # every ms Qcomp1 = np.interp(tcomp, t1, Q1) Qcomp2 = np.interp(tcomp, t2, Q2) Q_rmse = rmse(Qcomp1, Qcomp2) print('rmse = {:.5f} nC/cm2'.format(Q_rmse * 1e5))