diff --git a/PySONIC/parsers.py b/PySONIC/parsers.py index 15b1378..2ae5da1 100644 --- a/PySONIC/parsers.py +++ b/PySONIC/parsers.py @@ -1,560 +1,577 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2019-06-04 18:24:29 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-17 21:40:39 +# @Last Modified time: 2019-06-17 21:51:43 import logging import pprint import numpy as np from argparse import ArgumentParser from .utils import Intensity2Pressure, selectDirDialog, OpenFilesDialog, isIterable from .neurons import getPointNeuron class Parser(ArgumentParser): ''' Generic parser interface. ''' dist_str = '[scale min max n]' def __init__(self): super().__init__() self.pp = pprint.PrettyPrinter(indent=4) self.defaults = {} self.allowed = {} self.factors = {} self.to_parse = {} self.addPlot() self.addVerbose() def pprint(self, args): self.pp.pprint(args) def getDistribution(self, xmin, xmax, nx, scale='lin'): if scale == 'log': xmin, xmax = np.log10(xmin), np.log10(xmax) return {'lin': np.linspace, 'log': np.logspace}[scale](xmin, xmax, nx) def getDistFromList(self, xlist): if not isinstance(xlist, list): raise TypeError('Input must be a list') if len(xlist) != 4: raise ValueError('List must contain exactly 4 arguments ([type, min, max, n])') scale = xlist[0] if scale not in ('log', 'lin'): raise ValueError('Unknown distribution type (must be "lin" or "log")') xmin, xmax = [float(x) for x in xlist[1:-1]] if xmin >= xmax: raise ValueError('Specified minimum higher or equal than specified maximum') nx = int(xlist[-1]) if nx < 2: raise ValueError('Specified number must be at least 2') return self.getDistribution(xmin, xmax, nx, scale=scale) def addVerbose(self): self.add_argument( '-v', '--verbose', default=False, action='store_true', help='Increase verbosity') self.to_parse['loglevel'] = self.parseLogLevel def addPlot(self): self.add_argument( '-p', '--plot', type=str, nargs='+', help='Variables to plot') self.to_parse['pltscheme'] = self.parsePltScheme def addMPI(self): self.add_argument( '--mpi', default=False, action='store_true', help='Use multiprocessing') def addTest(self): self.add_argument( '--test', default=False, action='store_true', help='Run test configuration') def addSave(self): self.add_argument( '-s', '--save', default=False, action='store_true', help='Save output figure(s)') def addFigureExtension(self): self.add_argument( '--figext', type=str, default='png', help='Figure type (extension)') def addCompare(self, desc='Comparative graph'): self.add_argument( '--compare', default=False, action='store_true', help=desc) def addSamplingRate(self): self.add_argument( '--sr', type=int, default=1, help='Sampling rate for plot') def addSpikes(self): self.add_argument( '--spikes', type=str, default='none', help='How to indicate spikes on charge profile ("none", "marks" or "details")') def addNColumns(self): self.add_argument( '--ncol', type=int, default=1, help='Number of columns in figure') def addNLevels(self): self.add_argument( '--nlevels', type=int, default=10, help='Number of levels') def addHideOutput(self): self.add_argument( '--hide', default=False, action='store_true', help='Hide output') def addTimeRange(self, default=None): self.add_argument( '--trange', type=float, nargs=2, default=default, help='Time lower and upper bounds (ms)') self.to_parse['trange'] = self.parseTimeRange def addPotentialBounds(self, default=None): self.add_argument( '--Vbounds', type=float, nargs=2, default=default, help='Membrane potential lower and upper bounds (mV)') def addFiringRateBounds(self, default): self.add_argument( '--FRbounds', type=float, nargs=2, default=default, help='Firing rate lower and upper bounds (Hz)') def addFiringRateScale(self, default='lin'): self.add_argument( '--FRscale', type=str, choices=('lin', 'log'), default=default, help='Firing rate scale for plot ("lin" or "log")') def addCmap(self, default=None): self.add_argument( '--cmap', type=str, default=default, help='Colormap name') def addCscale(self, default='lin'): self.add_argument( '--cscale', type=str, default=default, choices=('lin', 'log'), help='Color scale ("lin" or "log")') def addInputDir(self, dep_key=None): self.inputdir_dep_key = dep_key self.add_argument( '-i', '--inputdir', type=str, help='Input directory') self.to_parse['inputdir'] = self.parseInputDir def addOutputDir(self, dep_key=None): self.outputdir_dep_key = dep_key self.add_argument( '-o', '--outputdir', type=str, help='Output directory') self.to_parse['outputdir'] = self.parseOutputDir def addInputFiles(self, dep_key=None): self.inputfiles_dep_key = dep_key self.add_argument( '-i', '--inputfiles', type=str, help='Input files') self.to_parse['inputfiles'] = self.parseInputFiles def addPatches(self): self.add_argument( '--patches', type=str, default='one', help='Stimulus patching mode ("none", "one", all", or a boolean list)') self.to_parse['patches'] = self.parsePatches def addThresholdCurve(self): self.add_argument( '--threshold', default=False, action='store_true', help='Show threshold amplitudes') def addNeuron(self): self.add_argument( '-n', '--neuron', type=str, nargs='+', help='Neuron name (string)') self.to_parse['neuron'] = self.parseNeuron def parseNeuron(self, args): return [getPointNeuron(n) for n in args['neuron']] def addInteractive(self): self.add_argument( '--interactive', default=False, action='store_true', help='Make interactive') def addLabels(self): self.add_argument( '--labels', type=str, nargs='+', default=None, help='Labels') def addRelativeTimeBounds(self): self.add_argument( '--rel_tbounds', type=float, nargs='+', default=None, help='Relative time lower and upper bounds') def addPretty(self): self.add_argument( '--pretty', default=False, action='store_true', help='Make figure pretty') def addSubset(self, choices): self.add_argument( '--subset', type=str, nargs='+', default='all', choices=choices, help='Run specific subset(s)') + self.subset_choices = choices + self.to_parse['subset'] = self.parseSubset + + def parseSubset(self, args): + if args['subset'] == ['all']: + return self.subset_choices + else: + return args['subset'] def parseTimeRange(self, args): if args['trange'] is None: return None return np.array(args['trange']) * 1e-3 def parsePatches(self, args): if args['patches'] not in ('none', 'one', 'all'): return eval(args['patches']) else: return args['patches'] def parseInputFiles(self, args): if self.inputfiles_dep_key is not None and not args[self.inputfiles_dep_key]: return None elif args['inputfiles'] is None: return OpenFilesDialog('pkl')[0] def parseDir(self, key, args, title, dep_key=None): if dep_key is not None and args[dep_key] is False: return None try: if args[key] is not None: return args[key] else: return selectDirDialog(title=title) except ValueError: raise ValueError('No {} selected'.format(key)) def parseInputDir(self, args): return self.parseDir( 'inputdir', args, 'Select input directory', self.inputdir_dep_key) def parseOutputDir(self, args): if hasattr(self, 'outputdir') and self.outputdir is not None: return self.outputdir else: return self.parseDir( 'outputdir', args, 'Select output directory', self.outputdir_dep_key) def parseLogLevel(self, args): return logging.DEBUG if args.pop('verbose') else logging.INFO def parsePltScheme(self, args): if args['plot'] is None or args['plot'] == ['all']: return None else: return {x: [x] for x in args['plot']} def restrict(self, args, keys): if sum([args[x] is not None for x in keys]) > 1: raise ValueError( 'You must provide only one of the following arguments: {}'.format(', '.join(keys))) def parse2array(self, args, key, factor=1): return np.array(args[key]) * factor def parse(self): args = vars(super().parse_args()) for k, v in self.defaults.items(): if k in args and args[k] is None: args[k] = v if isIterable(v) else [v] for k, parse_method in self.to_parse.items(): args[k] = parse_method(args) return args class TestParser(Parser): def __init__(self, valid_subsets): super().__init__() self.addProfiling() self.addSubset(valid_subsets + ['all']) def addProfiling(self): self.add_argument( '--profile', default=False, action='store_true', help='Run with profiling') +class FigureParser(Parser): + + def __init__(self, valid_subsets): + super().__init__() + self.addSubset(valid_subsets + ['all']) + self.addSave() + self.addOutputDir(dep_key='save') + + class PlotParser(Parser): def __init__(self): super().__init__() self.addHideOutput() self.addInputFiles() self.addOutputDir(dep_key='save') self.addSave() self.addFigureExtension() self.addCmap() self.addPretty() self.addTimeRange() self.addCscale() self.addLabels() class TimeSeriesParser(PlotParser): def __init__(self): super().__init__() self.addSpikes() self.addSamplingRate() self.addCompare() self.addPatches() class SimParser(Parser): ''' Generic simulation parser interface. ''' def __init__(self, outputdir=None): super().__init__() self.outputdir = outputdir self.addMPI() self.addOutputDir() def parse(self): args = super().parse() return args class MechSimParser(SimParser): ''' Parser to run mechanical simulations from the command line. ''' def __init__(self, *args, **kwargs): super().__init__(*args, **kwargs) self.defaults.update({ 'radius': 32.0, # nm 'embedding': 0., # um 'Cm0': getPointNeuron('RS').Cm0 * 1e2, # uF/m2 'Qm0': getPointNeuron('RS').Qm0 * 1e5, # nC/m2 'freq': 500.0, # kHz 'amp': 100.0, # kPa 'charge': 0., # nC/cm2 'fs': 100. # % }) self.factors.update({ 'radius': 1e-9, 'embedding': 1e-6, 'Cm0': 1e-2, 'Qm0': 1e-5, 'freq': 1e3, 'amp': 1e3, 'charge': 1e-5, 'fs': 1e-2 }) self.addRadius() self.addEmbedding() self.addCm0() self.addQm0() self.addFdrive() self.addAdrive() self.addCharge() self.addFs() def addRadius(self): self.add_argument( '-a', '--radius', nargs='+', type=float, help='Sonophore radius (nm)') def addEmbedding(self): self.add_argument( '--embedding', nargs='+', type=float, help='Embedding depth (um)') def addCm0(self): self.add_argument( '--Cm0', type=float, help='Resting membrane capacitance (uF/cm2)') def addQm0(self): self.add_argument( '--Qm0', type=float, help='Resting membrane charge density (nC/cm2)') def addFdrive(self): self.add_argument( '-f', '--freq', nargs='+', type=float, help='US frequency (kHz)') def addAdrive(self): self.add_argument( '-A', '--amp', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)') self.add_argument( '--Arange', type=str, nargs='+', help='Amplitude range {} (kPa)'.format(self.dist_str)) self.add_argument( '-I', '--intensity', nargs='+', type=float, help='Acoustic intensity (W/cm2)') self.add_argument( '--Irange', type=str, nargs='+', help='Intensity range {} (W/cm2)'.format(self.dist_str)) self.to_parse['amp'] = self.parseAmp def addAscale(self, default='lin'): self.add_argument( '--Ascale', type=str, choices=('lin', 'log'), default=default, help='Amplitude scale for plot ("lin" or "log")') def addCharge(self): self.add_argument( '-Q', '--charge', nargs='+', type=float, help='Membrane charge density (nC/cm2)') def addFs(self): self.add_argument( '--fs', nargs='+', type=float, help='Sonophore coverage fraction (%%)') self.add_argument( '--spanFs', default=False, action='store_true', help='Span Fs from 1 to 100%%') self.to_parse['fs'] = self.parseFs def parseAmp(self, args): params = ['Irange', 'Arange', 'intensity', 'amp'] self.restrict(args, params[:-1]) Irange, Arange, Int, Adrive = [args.pop(k) for k in params] if Irange is not None: amps = Intensity2Pressure(self.getDistFromList(Irange) * 1e4) # Pa elif Int is not None: amps = Intensity2Pressure(np.array(Int) * 1e4) # Pa elif Arange is not None: amps = self.getDistFromList(Arange) * self.factors['amp'] # Pa else: amps = np.array(Adrive) * self.factors['amp'] # Pa return amps def parseFs(self, args): if args.pop('spanFs', False): return np.arange(1, 101) * self.factors['fs'] # (-) else: return np.array(args['fs']) * self.factors['fs'] # (-) def parse(self): args = super().parse() for key in ['radius', 'embedding', 'Cm0', 'Qm0', 'freq', 'charge']: args[key] = self.parse2array(args, key, factor=self.factors[key]) return args class PWSimParser(SimParser): ''' Generic parser interface to run PW patterned simulations from the command line. ''' def __init__(self): super().__init__() self.defaults.update({ 'neuron': 'RS', 'tstim': 100.0, # ms 'toffset': 50., # ms 'PRF': 100.0, # Hz 'DC': 100.0 # % }) self.factors.update({ 'tstim': 1e-3, 'toffset': 1e-3, 'PRF': 1., 'DC': 1e-2 }) self.allowed.update({ 'DC': range(101) }) self.addNeuron() self.addTstim() self.addToffset() self.addPRF() self.addDC() self.addTitrate() self.addSpikes() def addTstim(self): self.add_argument( '-t', '--tstim', nargs='+', type=float, help='Stimulus duration (ms)') def addToffset(self): self.add_argument( '--toffset', nargs='+', type=float, help='Offset duration (ms)') def addPRF(self): self.add_argument( '--PRF', nargs='+', type=float, help='PRF (Hz)') def addDC(self): self.add_argument( '--DC', nargs='+', type=float, help='Duty cycle (%%)') self.add_argument( '--spanDC', default=False, action='store_true', help='Span DC from 1 to 100%%') self.to_parse['DC'] = self.parseDC def addTitrate(self): self.add_argument( '--titrate', default=False, action='store_true', help='Perform titration') def parseAmp(self, args): return NotImplementedError def parseDC(self, args): if args.pop('spanDC'): return np.arange(1, 101) * self.factors['DC'] # (-) else: return np.array(args['DC']) * self.factors['DC'] # (-) def parse(self, args=None): if args is None: args = super().parse() for key in ['tstim', 'toffset', 'PRF']: args[key] = self.parse2array(args, key, factor=self.factors[key]) return args class EStimParser(PWSimParser): ''' Parser to run E-STIM simulations from the command line. ''' def __init__(self): super().__init__() self.defaults.update({'amp': 10.0}) # mA/m2 self.factors.update({'amp': 1.}) self.addAstim() def addAstim(self): self.add_argument( '-A', '--amp', nargs='+', type=float, help='Amplitude of injected current density (mA/m2)') self.add_argument( '--Arange', type=str, nargs='+', help='Amplitude range {} (mA/m2)'.format(self.dist_str)) self.to_parse['amp'] = self.parseAmp def parseAmp(self, args): if args.pop('titrate'): return None Arange, Astim = [args.pop(k) for k in ['Arange', 'amp']] if Arange is not None: amps = self.getDistFromList(Arange) * self.factors['amp'] # mA/m2 else: amps = np.array(Astim) * self.factors['amp'] # mA/m2 return amps def parse(self): args = super().parse() return args class AStimParser(PWSimParser, MechSimParser): ''' Parser to run A-STIM simulations from the command line. ''' def __init__(self): MechSimParser.__init__(self) PWSimParser.__init__(self) self.defaults.update({'method': 'sonic'}) self.allowed.update({'method': ['full', 'hybrid', 'sonic']}) self.addMethod() def addMethod(self): self.add_argument( '-m', '--method', nargs='+', type=str, help='Numerical integration method ({})'.format(', '.join(self.allowed['method']))) self.to_parse['method'] = self.parseMethod def parseMethod(self, args): for item in args['method']: if item not in self.allowed['method']: raise ValueError('Unknown method type: "{}"'.format(item)) return args['method'] def parseAmp(self, args): if args.pop('titrate'): return None return MechSimParser.parseAmp(self, args) def parse(self): args = PWSimParser.parse(self, args=MechSimParser.parse(self)) for k in ['Cm0', 'Qm0', 'embedding', 'charge']: del args[k] return args diff --git a/paper figures/deprecated/figQSS.py b/paper figures/deprecated/figQSS.py index cdbdad0..b47fe59 100644 --- a/paper figures/deprecated/figQSS.py +++ b/paper figures/deprecated/figQSS.py @@ -1,296 +1,281 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-09-28 16:13:34 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-14 08:09:50 +# @Last Modified time: 2019-06-17 22:08:28 ''' Subpanels of the QSS approximation figure. ''' import os import logging import numpy as np import matplotlib.pyplot as plt import matplotlib import matplotlib.cm as cm from argparse import ArgumentParser from PySONIC.core import NeuronalBilayerSonophore from PySONIC.utils import logger, selectDirDialog from PySONIC.neurons import getPointNeuron +from PySONIC.parsers import FigureParser # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def plotQSSvars_vs_Adrive(neuron, a, Fdrive, PRF, DC, fs=8, markers=['-', '--', '.-'], title=None): pneuron = getPointNeuron(neuron) # Determine spiking threshold Vthr = pneuron.VT # mV Qthr = pneuron.Cm0 * Vthr * 1e-3 # C/m2 # Get QSS variables for each amplitude at threshold charge nbls = NeuronalBilayerSonophore(a, pneuron, Fdrive) Aref, _, Vmeff, QS_states = nbls.quasiSteadyStates(Fdrive, charges=Qthr, DCs=DC) # Compute US-ON and US-OFF ionic currents currents_on = pneuron.currents(Vmeff, QS_states) currents_off = pneuron.currents(pneuron.VT, QS_states) iNet_on = sum(currents_on.values()) iNet_off = sum(currents_off.values()) # Retrieve list of ionic currents names, with iLeak first ckeys = list(currents_on.keys()) ckeys.insert(0, ckeys.pop(ckeys.index('iLeak'))) # Compute quasi-steady ON, OFF and net charge variations, and threshold amplitude dQ_on = -iNet_on * DC / PRF dQ_off = -iNet_off * (1 - DC) / PRF dQ_net = dQ_on + dQ_off Athr = np.interp(0, dQ_net, Aref, left=0., right=np.nan) # Create figure fig, axes = plt.subplots(4, 1, figsize=(4, 6)) axes[-1].set_xlabel('Amplitude (kPa)', fontsize=fs) for ax in axes: for skey in ['top', 'right']: ax.spines[skey].set_visible(False) ax.set_xscale('log') ax.set_xlim(1e1, 1e2) ax.set_xticks([1e1, 1e2]) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(minor=True): item.set_visible(False) figname = '{} neuron thr dynamics {:.1f}nC_cm2 {:.0f}% DC'.format( pneuron.name, Qthr * 1e5, DC * 1e2) fig.suptitle(figname, fontsize=fs) # Subplot 1: Vmeff ax = axes[0] ax.set_ylabel('Effective potential (mV)', fontsize=fs) Vbounds = (-120, -40) ax.set_ylim(Vbounds) ax.set_yticks([Vbounds[0], pneuron.Vm0, Vbounds[1]]) ax.set_yticklabels(['{:.0f}'.format(Vbounds[0]), '$V_{m0}$', '{:.0f}'.format(Vbounds[1])]) ax.plot(Aref * 1e-3, Vmeff, '--', color='C0', label='ON') ax.plot(Aref * 1e-3, pneuron.VT * np.ones(Aref.size), ':', color='C0', label='OFF') ax.axhline(pneuron.Vm0, linewidth=0.5, color='k') # Subplot 2: quasi-steady states ax = axes[1] ax.set_ylabel('Quasi-steady states', fontsize=fs) ax.set_yticks([0, 0.5, 0.6]) ax.set_yticklabels(['0', '0.5', '1']) ax.set_ylim([-0.05, 0.65]) d = .01 f = 1.03 xcut = ax.get_xlim()[0] for ycut in [0.54, 0.56]: ax.plot([xcut / f, xcut * f], [ycut - d, ycut + d], color='k', clip_on=False) for label, QS_state in zip(pneuron.states, QS_states): if label == 'h': QS_state -= 0.4 ax.plot(Aref * 1e-3, QS_state, label=label) # Subplot 3: currents ax = axes[2] ax.set_ylabel('QSS Currents (mA/m2)', fontsize=fs) Ibounds = (-10, 10) ax.set_ylim(Ibounds) ax.set_yticks([Ibounds[0], 0.0, Ibounds[1]]) for i, key in enumerate(ckeys): c = 'C{}'.format(i) if isinstance(currents_off[key], float): currents_off[key] = np.ones(Aref.size) * currents_off[key] ax.plot(Aref * 1e-3, currents_on[key], '--', label=key, c=c) ax.plot(Aref * 1e-3, currents_off[key], ':', c=c) ax.plot(Aref * 1e-3, iNet_on, '--', color='k', label='iNet') ax.plot(Aref * 1e-3, iNet_off, ':', color='k') ax.axhline(0, color='k', linewidth=0.5) # Subplot 4: charge variations and activation threshold ax = axes[3] ax.set_ylabel('$\\rm \Delta Q_{QS}\ (nC/cm^2)$', fontsize=fs) dQbounds = (-0.06, 0.1) ax.set_ylim(dQbounds) ax.set_yticks([dQbounds[0], 0.0, dQbounds[1]]) ax.plot(Aref * 1e-3, dQ_on, '--', color='C0', label='ON') ax.plot(Aref * 1e-3, dQ_off, ':', color='C0', label='OFF') ax.plot(Aref * 1e-3, dQ_net, color='C0', label='Net') ax.plot([Athr * 1e-3] * 2, [ax.get_ylim()[0], 0], linestyle='--', color='k') ax.plot([Athr * 1e-3], [0], 'o', c='k') ax.axhline(0, color='k', linewidth=0.5) fig.tight_layout() fig.subplots_adjust(right=0.8) for ax in axes: ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5)) if title is not None: fig.canvas.set_window_title(title) return fig def plotQSSdQ_vs_Adrive(neuron, a, Fdrive, PRF, DCs, fs=8, title=None): pneuron = getPointNeuron(neuron) # Determine spiking threshold Vthr = pneuron.VT # mV Qthr = pneuron.Cm0 * Vthr * 1e-3 # C/m2 # Get QSS variables for each amplitude and DC at threshold charge nbls = NeuronalBilayerSonophore(a, pneuron, Fdrive) Aref, _, Vmeff, QS_states = nbls.quasiSteadyStates(Fdrive, charges=Qthr, DCs=DCs) dQnet = np.empty((DCs.size, Aref.size)) Athr = np.empty(DCs.size) for i, DC in enumerate(DCs): # Compute US-ON and US-OFF net membrane current from QSS variables iNet_on = pneuron.iNet(Vmeff, QS_states[:, :, i]) iNet_off = pneuron.iNet(Vthr, QS_states[:, :, i]) # Compute the pulse average net current along the amplitude space iNet_avg = iNet_on * DC + iNet_off * (1 - DC) dQnet[i, :] = -iNet_avg / PRF # Find the threshold amplitude that cancels the pulse average net current Athr[i] = np.interp(0, -iNet_avg, Aref, left=0., right=np.nan) # Create figure fig, ax = plt.subplots(figsize=(4, 2)) figname = '{} neuron thr vs DC'.format(pneuron.name, Qthr * 1e5) fig.suptitle(figname, fontsize=fs) for key in ['top', 'right']: ax.spines[key].set_visible(False) ax.set_xscale('log') for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(minor=True): item.set_visible(False) ax.set_xlabel('Amplitude (kPa)', fontsize=fs) ax.set_ylabel('$\\rm \Delta Q_{QS}\ (nC/cm^2)$', fontsize=fs) ax.set_xlim(1e1, 1e2) ax.axhline(0., linewidth=0.5, color='k') ax.set_ylim(-0.06, 0.12) ax.set_yticks([-0.05, 0.0, 0.10]) ax.set_yticklabels(['-0.05', '0', '0.10']) norm = matplotlib.colors.LogNorm(DCs.min(), DCs.max()) sm = cm.ScalarMappable(norm=norm, cmap='viridis') sm._A = [] for i, DC in enumerate(DCs): ax.plot(Aref * 1e-3, dQnet[i, :], c=sm.to_rgba(DC), label='{:.0f}% DC'.format(DC * 1e2)) ax.plot([Athr[i] * 1e-3] * 2, [ax.get_ylim()[0], 0], linestyle='--', c=sm.to_rgba(DC)) ax.plot([Athr[i] * 1e-3], [0], 'o', c=sm.to_rgba(DC)) fig.tight_layout() fig.subplots_adjust(right=0.8) ax.legend(loc='center right', fontsize=fs, frameon=False, bbox_to_anchor=(1.3, 0.5)) if title is not None: fig.canvas.set_window_title(title) return fig def plotQSSAthr_vs_DC(neurons, a, Fdrive, DCs_dense, DCs_sparse, fs=8, title=None): fig, ax = plt.subplots(figsize=(3, 3)) ax.set_title('Rheobase amplitudes', fontsize=fs) ax.set_xlabel('Duty cycle (%)', fontsize=fs) ax.set_ylabel('$\\rm A_T\ (kPa)$', fontsize=fs) for key in ['top', 'right']: ax.spines[key].set_visible(False) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) ax.set_xticks([25, 50, 75, 100]) ax.set_yscale('log') ax.set_ylim([10, 600]) norm = matplotlib.colors.LogNorm(DCs_sparse.min(), DCs_sparse.max()) sm = cm.ScalarMappable(norm=norm, cmap='viridis') sm._A = [] for i, neuron in enumerate(neurons): pneuron = getPointNeuron(neuron) nbls = NeuronalBilayerSonophore(a, pneuron) Athrs_dense = nbls.findRheobaseAmps(DCs_dense, Fdrive, pneuron.VT)[0] * 1e-3 # kPa Athrs_sparse = nbls.findRheobaseAmps(DCs_sparse, Fdrive, pneuron.VT)[0] * 1e-3 # kPa ax.plot(DCs_dense * 1e2, Athrs_dense, label='{} neuron'.format(pneuron.name)) for DC, Athr in zip(DCs_sparse, Athrs_sparse): ax.plot(DC * 1e2, Athr, 'o', label='{:.0f}% DC'.format(DC * 1e2) if i == len(neurons) - 1 else None, c=sm.to_rgba(DC)) ax.legend(fontsize=fs, frameon=False) fig.tight_layout() if title is not None: fig.canvas.set_window_title(title) return fig def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') - ap.add_argument('-o', '--outdir', type=str, help='Output directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - figset = args.figset - if figset == 'all': - figset = ['a', 'b', 'c', 'e'] + parser = FigureParser(['a', 'b', 'c', 'd', 'e']) + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters a = 32e-9 # m Fdrive = 500e3 # Hz PRF = 100.0 # Hz DC = 0.5 DCs_sparse = np.array([5, 15, 50, 75, 95]) / 1e2 DCs_dense = np.arange(1, 101) / 1e2 # Figures figs = [] if 'a' in figset: figs += [ plotQSSvars_vs_Adrive('RS', a, Fdrive, PRF, DC, title=figbase + 'a RS'), plotQSSvars_vs_Adrive('LTS', a, Fdrive, PRF, DC, title=figbase + 'a LTS') ] if 'b' in figset: figs += [ plotQSSdQ_vs_Adrive('RS', a, Fdrive, PRF, DCs_sparse, title=figbase + 'b RS'), plotQSSdQ_vs_Adrive('LTS', a, Fdrive, PRF, DCs_sparse, title=figbase + 'b LTS') ] if 'c' in figset: figs.append(plotQSSAthr_vs_DC(['RS', 'LTS'], a, Fdrive, DCs_dense, DCs_sparse, title=figbase + 'c')) - if args.save: - try: - outdir = selectDirDialog() if args.outdir is None else args.outdir - except ValueError as err: - logger.error(err) - return + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(outdir, figname), transparent=True) + fig.savefig(os.path.join(args['outputdir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig2.py b/paper figures/fig2.py index de1f490..2075e2e 100644 --- a/paper figures/fig2.py +++ b/paper figures/fig2.py @@ -1,331 +1,313 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-06-06 18:38:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-17 14:24:17 +# @Last Modified time: 2019-06-17 21:55:40 ''' Sub-panels of the model optimization figure. ''' import os -import logging import numpy as np import matplotlib import matplotlib.pyplot as plt from matplotlib.ticker import FormatStrFormatter from matplotlib.patches import Rectangle -from argparse import ArgumentParser -from PySONIC.utils import logger, rescale, si_format, selectDirDialog +from PySONIC.utils import logger, rescale, si_format from PySONIC.plt import GroupedTimeSeries, cm2inch from PySONIC.constants import NPC_DENSE from PySONIC.neurons import getPointNeuron from PySONIC.core import BilayerSonophore, NeuronalBilayerSonophore +from PySONIC.parsers import FigureParser # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def PmApprox(bls, Z, fs=12, lw=2): fig, ax = plt.subplots(figsize=cm2inch(7, 7)) for key in ['right', 'top']: ax.spines[key].set_visible(False) for key in ['bottom', 'left']: ax.spines[key].set_linewidth(2) ax.spines['bottom'].set_position('zero') ax.set_xlabel('Z (nm)', fontsize=fs) ax.set_ylabel('Pressure (kPa)', fontsize=fs, labelpad=-10) ax.set_xticks([0, bls.a * 1e9]) ax.set_xticklabels(['0', 'a']) ax.tick_params(axis='x', which='major', length=25, pad=5) ax.set_yticks([0]) ax.set_ylim([-10, 50]) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) ax.plot(Z * 1e9, bls.v_PMavg(Z, bls.v_curvrad(Z), bls.surface(Z)) * 1e-3, c='g', label='$P_m$') ax.plot(Z * 1e9, bls.PMavgpred(Z) * 1e-3, '--', c='r', label='$\~P_m$') ax.axhline(y=0, color='k') ax.legend(fontsize=fs, frameon=False) fig.tight_layout() fig.canvas.set_window_title(figbase + 'a') return fig def recasting(nbls, Fdrive, Adrive, fs=12, lw=2, ps=15): # Run effective simulation data, _ = nbls.simulate(Fdrive, Adrive, 5 / Fdrive, 0., method='full') t, Qm, Vm = [data[key].values for key in ['t', 'Qm', 'Vm']] t *= 1e6 # us Qm *= 1e5 # nC/cm2 Qrange = (Qm.min(), Qm.max()) dQ = Qrange[1] - Qrange[0] # Create figure and axes fig, axes = plt.subplots(1, 2, figsize=cm2inch(17, 5)) for ax in axes: ax.set_xticks([]) ax.set_yticks([]) # Plot Q-trace and V-trace ax = axes[0] for key in ['top', 'right']: ax.spines[key].set_visible(False) for key in ['bottom', 'left']: ax.spines[key].set_position(('axes', -0.03)) ax.spines[key].set_linewidth(2) ax.plot(t, Vm, label='Vm', c='dimgrey', linewidth=lw) ax.plot(t, Qm, label='Qm', c='k', linewidth=lw) ax.add_patch(Rectangle( (t[0], Qrange[0] - 5), t[-1], dQ + 10, fill=False, edgecolor='k', linestyle='--', linewidth=1 )) ax.yaxis.set_tick_params(width=2) ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f')) # ax.set_xlim((t.min(), t.max())) ax.set_xticks([]) ax.set_xlabel('{}s'.format(si_format((t.max()), space=' ')), fontsize=fs) ax.set_ylabel('$\\rm nC/cm^2$ - mV', fontsize=fs, labelpad=-15) ax.set_yticks(ax.get_ylim()) for item in ax.get_yticklabels(): item.set_fontsize(fs) # Plot inset on Q-trace ax = axes[1] for key in ['top', 'right', 'bottom', 'left']: ax.spines[key].set_linewidth(1) ax.spines[key].set_linestyle('--') ax.plot(t, Vm, label='Vm', c='dimgrey', linewidth=lw) ax.plot(t, Qm, label='Qm', c='k', linewidth=lw) ax.set_xlim((t.min(), t.max())) ax.set_xticks([]) ax.set_yticks([]) delta = 0.05 ax.set_ylim(Qrange[0] - delta * dQ, Qrange[1] + delta * dQ) fig.canvas.set_window_title(figbase + 'b') return fig def mechSim(bls, Fdrive, Adrive, Qm, fs=12, lw=2, ps=15): # Run mechanical simulation data, _ = bls.simulate(Fdrive, Adrive, Qm) t, Z, ng = [data[key].values for key in ['t', 'Z', 'ng']] # Create figure fig, ax = plt.subplots(figsize=cm2inch(7, 7)) fig.suptitle('Mechanical simulation', fontsize=12) for skey in ['bottom', 'left', 'right', 'top']: ax.spines[skey].set_visible(False) ax.set_xticks([]) ax.set_yticks([]) # Plot variables and labels t_plot = np.insert(t, 0, -1e-6) * 1e6 Pac = Adrive * np.sin(2 * np.pi * Fdrive * t + np.pi) # Pa yvars = {'P_A': Pac * 1e-3, 'Z': Z * 1e9, 'n_g': ng * 1e22} colors = {'P_A': 'k', 'Z': 'C0', 'n_g': 'C5'} dy = 1.2 for i, ykey in enumerate(yvars.keys()): y = yvars[ykey] y_plot = rescale(np.insert(y, 0, y[0])) - dy * i ax.plot(t_plot, y_plot, color=colors[ykey], linewidth=lw) ax.text(t_plot[0] - 0.1, y_plot[0], '$\mathregular{{{}}}$'.format(ykey), fontsize=fs, horizontalalignment='right', verticalalignment='center', color=colors[ykey]) # Acoustic pressure annotations ax.annotate(s='', xy=(1.5, 1.1), xytext=(3.5, 1.1), arrowprops=dict(arrowstyle='<|-|>', color='k')) ax.text(2.5, 1.12, '1/f', fontsize=fs, color='k', horizontalalignment='center', verticalalignment='bottom') ax.annotate(s='', xy=(1.5, -0.1), xytext=(1.5, 1), arrowprops=dict(arrowstyle='<|-|>', color='k')) ax.text(1.55, 0.4, '2A', fontsize=fs, color='k', horizontalalignment='left', verticalalignment='center') # Periodic stabilization patch ax.add_patch(Rectangle((2, -2 * dy - 0.1), 2, 2 * dy, color='dimgrey', alpha=0.3)) ax.text(3, -2 * dy - 0.2, 'limit cycle', fontsize=fs, color='dimgrey', horizontalalignment='center', verticalalignment='top') # Z_last patch ax.add_patch(Rectangle((2, -dy - 0.1), 2, dy, edgecolor='k', facecolor='none', linestyle='--')) # ngeff annotations c = plt.get_cmap('tab20').colors[11] ax.text(t_plot[-1] + 0.1, y_plot[-1], '$\mathregular{n_{g,eff}}$', fontsize=fs, color=c, horizontalalignment='left', verticalalignment='center') ax.scatter([t_plot[-1]], [y_plot[-1]], color=c, s=ps) fig.canvas.set_window_title(figbase + 'c mechsim') return fig def cycleAveraging(bls, pneuron, Fdrive, Adrive, Qm, fs=12, lw=2, ps=15): # Run mechanical simulation data, _ = bls.simulate(Fdrive, Adrive, Qm) t, Z, ng = [data[key].values for key in ['t', 'Z', 'ng']] # Compute variables evolution over last acoustic cycle t_last = t[-NPC_DENSE:] * 1e6 # us Z_last = Z[-NPC_DENSE:] # m Cm = bls.v_Capct(Z_last) * 1e2 # uF/m2 Vm = Qm / Cm * 1e5 # mV yvars = { 'C_m': Cm, # uF/cm2 'V_m': Vm, # mV '\\alpha_m': pneuron.alpham(Vm) * 1e3, # ms-1 '\\beta_m': pneuron.betam(Vm) * 1e3, # ms-1 'p_\\infty / \\tau_p': pneuron.pinf(Vm) / pneuron.taup(Vm) * 1e3, # ms-1 '(1-p_\\infty) / \\tau_p': (1 - pneuron.pinf(Vm)) / pneuron.taup(Vm) * 1e3 # ms-1 } # Determine colors violets = plt.get_cmap('Paired').colors[8:10][::-1] oranges = plt.get_cmap('Paired').colors[6:8][::-1] colors = { 'C_m': ['k', 'dimgrey'], 'V_m': plt.get_cmap('tab20').colors[14:16], '\\alpha_m': violets, '\\beta_m': oranges, 'p_\\infty / \\tau_p': violets, '(1-p_\\infty) / \\tau_p': oranges } # Create figure and axes fig, axes = plt.subplots(6, 1, figsize=cm2inch(4, 15)) fig.suptitle('Cycle-averaging', fontsize=fs) for ax in axes: ax.set_xticks([]) ax.set_yticks([]) for skey in ['bottom', 'left', 'right', 'top']: ax.spines[skey].set_visible(False) # Plot variables for ax, ykey in zip(axes, yvars.keys()): ax.set_xticks([]) ax.set_yticks([]) for skey in ['bottom', 'left', 'right', 'top']: ax.spines[skey].set_visible(False) y = yvars[ykey] ax.plot(t_last, y, color=colors[ykey][0], linewidth=lw) ax.plot([t_last[0], t_last[-1]], [np.mean(y)] * 2, '--', color=colors[ykey][1]) ax.scatter([t_last[-1]], [np.mean(y)], s=ps, color=colors[ykey][1]) ax.text(t_last[0] - 0.1, y[0], '$\mathregular{{{}}}$'.format(ykey), fontsize=fs, horizontalalignment='right', verticalalignment='center', color=colors[ykey][0]) fig.canvas.set_window_title(figbase + 'c cycleavg') return fig def Qsolution(nbls, Fdrive, Adrive, tstim, toffset, PRF, DC, fs=12, lw=2, ps=15): # Run effective simulation data, _ = nbls.simulate(Fdrive, Adrive, tstim, toffset, PRF, DC, method='sonic') t, Qm, states = [data[key].values for key in ['t', 'Qm', 'stimstate']] t *= 1e3 # ms Qm *= 1e5 # nC/cm2 - _, tpulse_on, tpulse_off = GroupedTimeSeries.getStimPulses(_, t, states) + tpulse_on, tpulse_off = GroupedTimeSeries.getStimPulses(_, t, states) # Add small onset t = np.insert(t, 0, -5.0) Qm = np.insert(Qm, 0, Qm[0]) # Create figure and axes fig, ax = plt.subplots(figsize=cm2inch(12, 6)) ax.set_xticks([]) ax.set_yticks([]) for key in ['top', 'right']: ax.spines[key].set_visible(False) for key in ['bottom', 'left']: ax.spines[key].set_position(('axes', -0.03)) ax.spines[key].set_linewidth(2) # Plot Q-trace and stimulation pulses ax.plot(t, Qm, label='Qm', c='k', linewidth=lw) for ton, toff in zip(tpulse_on, tpulse_off): ax.axvspan(ton, toff, edgecolor='none', facecolor='#8A8A8A', alpha=0.2) ax.yaxis.set_tick_params(width=2) ax.yaxis.set_major_formatter(FormatStrFormatter('%.0f')) ax.set_xlim((t.min(), t.max())) ax.set_xticks([]) ax.set_xlabel('{}s'.format(si_format((t.max()) * 1e-3, space=' ')), fontsize=fs) ax.set_ylabel('$\\rm nC/cm^2$', fontsize=fs, labelpad=-15) ax.set_yticks(ax.get_ylim()) for item in ax.get_yticklabels(): item.set_fontsize(fs) ax.legend(fontsize=fs, frameon=False) fig.canvas.set_window_title(figbase + 'e Qtrace') return fig def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') - ap.add_argument('-o', '--outdir', type=str, help='Output directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - figset = args.figset - if figset == 'all': - figset = ['a', 'b', 'c', 'e'] - + parser = FigureParser(['a', 'b', 'c', 'e']) + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters pneuron = getPointNeuron('RS') a = 32e-9 # m Fdrive = 500e3 # Hz Adrive = 100e3 # Pa PRF = 100. # Hz DC = 0.5 tstim = 150e-3 # s toffset = 100e-3 # s Qm = -71.9e-5 # C/cm2 bls = BilayerSonophore(a, pneuron.Cm0, pneuron.Qm0) nbls = NeuronalBilayerSonophore(a, pneuron) # Figures figs = [] if 'a' in figset: figs.append(PmApprox(bls, np.linspace(-0.4 * bls.Delta_, bls.a, 1000))) if 'b' in figset: figs.append(recasting(nbls, Fdrive, Adrive)) if 'c' in figset: figs += [ mechSim(bls, Fdrive, Adrive, Qm), cycleAveraging(bls, pneuron, Fdrive, Adrive, Qm) ] if 'e' in figset: figs.append(Qsolution(nbls, Fdrive, Adrive, tstim, toffset, PRF, DC)) - if args.save: - try: - outdir = selectDirDialog() if args.outdir is None else args.outdir - except ValueError as err: - logger.error(err) - return + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(outdir, figname), transparent=True) + fig.savefig(os.path.join(args['outputdir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig4.py b/paper figures/fig4.py index 540037e..34ac908 100644 --- a/paper figures/fig4.py +++ b/paper figures/fig4.py @@ -1,85 +1,66 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-02-15 15:59:37 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-14 08:07:29 +# @Last Modified time: 2019-06-17 21:55:47 ''' Sub-panels of the effective variables figure. ''' import os import matplotlib import matplotlib.pyplot as plt -from argparse import ArgumentParser -import logging from PySONIC.plt import plotEffectiveVariables -from PySONIC.utils import logger, selectDirDialog +from PySONIC.utils import logger from PySONIC.neurons import getPointNeuron +from PySONIC.parsers import FigureParser # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - ap.add_argument('-o', '--outdir', type=str, help='Output directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - figset = args.figset - if figset == 'all': - figset = ['a', 'b', 'c'] - + parser = FigureParser(['a', 'b', 'c']) + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters pneuron = getPointNeuron('RS') a = 32e-9 # m Fdrive = 500e3 # Hz Adrive = 50e3 # Pa # Generate figures figs = [] if 'a' in figset: fig = plotEffectiveVariables(pneuron, a=a, Fdrive=Fdrive, cmap='Oranges', zscale='log') fig.canvas.set_window_title(figbase + 'a') figs.append(fig) if 'b' in figset: fig = plotEffectiveVariables(pneuron, a=a, Adrive=Adrive, cmap='Greens', zscale='log') fig.canvas.set_window_title(figbase + 'b') figs.append(fig) if 'c' in figset: fig = plotEffectiveVariables( pneuron, Fdrive=Fdrive, Adrive=Adrive, cmap='Blues', zscale='log') fig.canvas.set_window_title(figbase + 'c') figs.append(fig) - if args.save: - try: - outdir = selectDirDialog() if args.outdir is None else args.outdir - except ValueError as err: - logger.error(err) - return + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(outdir, figname), transparent=True) + fig.savefig(os.path.join(args['outputdir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig5.py b/paper figures/fig5.py index ccea95d..8bc17ef 100644 --- a/paper figures/fig5.py +++ b/paper figures/fig5.py @@ -1,445 +1,432 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-06-06 18:38:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-17 14:14:12 +# @Last Modified time: 2019-06-17 22:00:32 ''' Sub-panels of the NICE and SONIC accuracies comparative figure. ''' import os -import logging import numpy as np import matplotlib import matplotlib.pyplot as plt -from argparse import ArgumentParser from PySONIC.utils import * from PySONIC.neurons import * from PySONIC.plt import CompTimeSeries, cm2inch +from PySONIC.parsers import FigureParser from utils import * # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def plotSpikingMetrics(xvar, xlabel, metrics_dict, logscale=False, spikeamp=True, colors=None, fs=8, lw=2, ps=4, figsize=cm2inch(7.25, 5.8)): ''' Plot the evolution of spiking metrics as function of a specific stimulation parameter. ''' ls = {'full': 'o-', 'sonic': 'o--'} cdefault = {'full': 'silver', 'sonic': 'k'} # Create figure fig, axes = plt.subplots(3, 1, figsize=figsize) ibase = 0 if spikeamp else 1 axes[ibase].set_ylabel('Latency\n (ms)', fontsize=fs, rotation=0, ha='right', va='center') axes[ibase + 1].set_ylabel( 'Firing\n rate (Hz)', fontsize=fs, rotation=0, ha='right', va='center') if spikeamp: axes[2].set_ylabel('Spike amp.\n ($\\rm nC/cm^2$)', fontsize=fs, rotation=0, ha='right', va='center') for ax in axes: ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) if logscale: ax.set_xscale('log') for item in ax.get_yticklabels(): item.set_fontsize(fs) for ax in axes[:-1]: ax.spines['bottom'].set_visible(False) ax.set_xticks([]) plt.setp(ax.get_xticklabels(minor=True), visible=False) ax.get_xaxis().set_tick_params(which='minor', size=0) ax.get_xaxis().set_tick_params(which='minor', width=0) axes[-1].set_xlabel(xlabel, fontsize=fs) if not logscale: axes[-1].set_xticks([min(xvar), max(xvar)]) for item in axes[-1].get_xticklabels(): item.set_fontsize(fs) # Plot metrics for each neuron for i, neuron in enumerate(metrics_dict.keys()): full_metrics = metrics_dict[neuron]['full'] sonic_metrics = metrics_dict[neuron]['sonic'] c = colors[neuron] if colors is not None else cdefault # Latency rf = 10 ax = axes[ibase] ax.plot(xvar, full_metrics['latencies (ms)'].values, ls['full'], color=c['full'], linewidth=lw, markersize=ps) ax.plot(xvar, sonic_metrics['latencies (ms)'].values, ls['sonic'], color=c['sonic'], linewidth=lw, markersize=ps, label=neuron) # Firing rate rf = 10 ax = axes[ibase + 1] ax.errorbar(xvar, full_metrics['mean firing rates (Hz)'].values, yerr=full_metrics['std firing rates (Hz)'].values, fmt=ls['full'], color=c['full'], linewidth=lw, markersize=ps) ax.errorbar(xvar, sonic_metrics['mean firing rates (Hz)'].values, yerr=sonic_metrics['std firing rates (Hz)'].values, fmt=ls['sonic'], color=c['sonic'], linewidth=lw, markersize=ps) # Spike amplitudes if spikeamp: ax = axes[2] rf = 10 ax.errorbar(xvar, full_metrics['mean spike amplitudes (nC/cm2)'].values, yerr=full_metrics['std spike amplitudes (nC/cm2)'].values, fmt=ls['full'], color=c['full'], linewidth=lw, markersize=ps) ax.errorbar(xvar, sonic_metrics['mean spike amplitudes (nC/cm2)'].values, yerr=sonic_metrics['std spike amplitudes (nC/cm2)'].values, fmt=ls['sonic'], color=c['sonic'], linewidth=lw, markersize=ps) # Adapt axes y-limits rf = 10 for ax in axes: ax.set_ylim([np.floor(ax.get_ylim()[0] / rf) * rf, np.ceil(ax.get_ylim()[1] / rf) * rf]) ax.set_yticks([max(ax.get_ylim()[0], 0), ax.get_ylim()[1]]) # Legend if len(metrics_dict.keys()) > 1: leg = axes[0].legend(fontsize=fs, frameon=False, bbox_to_anchor=(0., 0.9, 1., .102), loc=8, ncol=2, borderaxespad=0.) for l in leg.get_lines(): l.set_linestyle('-') fig.subplots_adjust(hspace=.3, bottom=0.2, left=0.35, right=0.95, top=0.95) return fig def Qprofiles_vs_amp(neuron, a, Fdrive, CW_Athrs, tstim, toffset, inputdir): ''' Comparison of resulting charge profiles for CW stimuli at sub-threshold, threshold and supra-threshold amplitudes. ''' Athr = CW_Athrs[neuron].loc[Fdrive * 1e-3] # kPa amps = np.array([Athr - 5., Athr, Athr + 20.]) * 1e3 # Pa subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], amps, [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], amps, [tstim], [toffset], [None], [1.], 'full')) regimes = ['AT - 5 kPa', 'AT', 'AT + 20 kPa'] comp_plot = CompTimeSeries(sum([[x, y] for x, y in zip(full_fpaths, sonic_fpaths)], []), 'Qm') fig = comp_plot.render( labels=sum([['', x] for x in regimes], []), lines=['-', '--'] * len(regimes), colors=plt.get_cmap('Paired').colors[:2 * len(regimes)], fs=8, patches='one', xticks=[0, 250], yticks=[getPointNeuron(neuron).Vm0, 25], straightlegend=True, figsize=cm2inch(12.5, 5.8) ) fig.axes[0].get_xaxis().set_label_coords(0.5, -0.05) fig.subplots_adjust(bottom=0.2, right=0.95, top=0.95) fig.canvas.set_window_title(figbase + 'a Qprofiles') return fig def spikemetrics_vs_amp(neuron, a, Fdrive, amps, tstim, toffset, inputdir): ''' Comparison of spiking metrics for CW stimuli at various supra-threshold amplitudes. ''' subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], amps, [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], amps, [tstim], [toffset], [None], [1.], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} metrics_files = {x: '{}_spikemetrics_vs_amplitude_{}.csv'.format(neuron, x) for x in ['full', 'sonic']} metrics_fpaths = {key: os.path.join(inputdir, value) for key, value in metrics_files.items()} xlabel = 'Amplitude (kPa)' metrics = getSpikingMetrics( subdir, neuron, amps * 1e-3, xlabel, data_fpaths, metrics_fpaths) fig = plotSpikingMetrics(amps * 1e-3, xlabel, {neuron: metrics}, logscale=True) fig.canvas.set_window_title(figbase + 'a spikemetrics') return fig def Qprofiles_vs_freq(neuron, a, freqs, CW_Athrs, tstim, toffset, inputdir): ''' Comparison of resulting charge profiles for supra-threshold CW stimuli at low and high US frequencies. ''' subdir = os.path.join(inputdir, neuron) sonic_fpaths, full_fpaths = [], [] for Fdrive in freqs: Athr = CW_Athrs[neuron].loc[Fdrive * 1e-3] # kPa Adrive = (Athr + 20.) * 1e3 # Pa sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'full')) comp_plot = CompTimeSeries(sum([[x, y] for x, y in zip(full_fpaths, sonic_fpaths)], []), 'Qm') fig = comp_plot.render( labels=sum([['', '{}Hz'.format(si_format(f))] for f in freqs], []), lines=['-', '--'] * len(freqs), colors=plt.get_cmap('Paired').colors[6:10], fs=8, patches='one', xticks=[0, 250], yticks=[getPointNeuron(neuron).Vm0, 25], straightlegend=True, figsize=cm2inch(12.5, 5.8), inset={'xcoords': [5, 40], 'ycoords': [-35, 45], 'xlims': [57.5, 58.5], 'ylims': [10, 35]} ) fig.axes[0].get_xaxis().set_label_coords(0.5, -0.05) fig.subplots_adjust(bottom=0.2, right=0.95, top=0.95) fig.canvas.set_window_title(figbase + 'b Qprofiles') return fig def spikemetrics_vs_freq(neuron, a, freqs, CW_Athrs, tstim, toffset, inputdir): ''' Comparison of spiking metrics for supra-threshold CW stimuli at various US frequencies. ''' subdir = os.path.join(inputdir, neuron) sonic_fpaths, full_fpaths = [], [] for Fdrive in freqs: Athr = CW_Athrs[neuron].loc[Fdrive * 1e-3] # kPa Adrive = (Athr + 20.) * 1e3 # Pa sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} metrics_files = {x: '{}_spikemetrics_vs_frequency_{}.csv'.format(neuron, x) for x in ['full', 'sonic']} metrics_fpaths = {key: os.path.join(inputdir, value) for key, value in metrics_files.items()} xlabel = 'Frequency (kHz)' metrics = getSpikingMetrics( subdir, neuron, freqs * 1e-3, xlabel, data_fpaths, metrics_fpaths) fig = plotSpikingMetrics(freqs * 1e-3, xlabel, {neuron: metrics}, logscale=True) fig.canvas.set_window_title(figbase + 'b spikemetrics') return fig def Qprofiles_vs_radius(neuron, radii, Fdrive, CW_Athrs, tstim, toffset, inputdir): ''' Comparison of resulting charge profiles for supra-threshold CW stimuli for small and large sonophore radii. ''' subdir = os.path.join(inputdir, neuron) sonic_fpaths, full_fpaths = [], [] for a in radii: Athr = CW_Athrs[neuron].loc[a * 1e9] # kPa Adrive = (Athr + 20.) * 1e3 # Pa sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'full')) tmp = plt.get_cmap('Paired').colors colors = tmp[2:4] + tmp[10:12] comp_plot = CompTimeSeries(sum([[x, y] for x, y in zip(full_fpaths, sonic_fpaths)], []), 'Qm') fig = comp_plot.render( labels=sum([['', '{:.0f} nm'.format(a * 1e9)] for a in radii], []), lines=['-', '--'] * len(radii), colors=colors, fs=8, patches='one', xticks=[0, 250], yticks=[getPointNeuron(neuron).Vm0, 25], straightlegend=True, figsize=cm2inch(12.5, 5.8) ) fig.axes[0].get_xaxis().set_label_coords(0.5, -0.05) fig.subplots_adjust(bottom=0.2, right=0.95, top=0.95) fig.canvas.set_window_title(figbase + 'c Qprofiles') return fig def spikemetrics_vs_radius(neuron, radii, Fdrive, CW_Athrs, tstim, toffset, inputdir): ''' Comparison of spiking metrics for supra-threshold CW stimuli with various sonophore diameters. ''' subdir = os.path.join(inputdir, neuron) sonic_fpaths, full_fpaths = [], [] for a in radii: Athr = CW_Athrs[neuron].loc[np.round(a * 1e9, 1)] # kPa Adrive = (Athr + 20.) * 1e3 # Pa sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} metrics_files = {x: '{}_spikemetrics_vs_radius_{}.csv'.format(neuron, x) for x in ['full', 'sonic']} metrics_fpaths = {key: os.path.join(inputdir, value) for key, value in metrics_files.items()} xlabel = 'Sonophore radius (nm)' metrics = getSpikingMetrics( subdir, neuron, radii * 1e9, xlabel, data_fpaths, metrics_fpaths) fig = plotSpikingMetrics(radii * 1e9, xlabel, {neuron: metrics}, logscale=True) fig.canvas.set_window_title(figbase + 'c spikemetrics') return fig def Qprofiles_vs_DC(neurons, a, Fdrive, Adrive, tstim, toffset, PRF, DC, inputdir): ''' Comparison of resulting charge profiles for PW stimuli at 5% duty cycle for different neuron types. ''' sonic_fpaths, full_fpaths = [], [] for neuron in neurons: subdir = os.path.join(inputdir, neuron) sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [PRF], [DC], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [PRF], [DC], 'full')) colors = list(plt.get_cmap('Paired').colors[:6]) del colors[2:4] comp_plot = CompTimeSeries(sum([[x, y] for x, y in zip(full_fpaths, sonic_fpaths)], []), 'Qm') fig = comp_plot.render( labels=sum([['', '{}, {:.0f}% DC'.format(x, DC * 1e2)] for x in neurons], []), lines=['-', '--'] * len(neurons), colors=colors, fs=8, patches='one', xticks=[0, 250], yticks=[min(getPointNeuron(n).Vm0 for n in neurons), 50], straightlegend=True, figsize=cm2inch(12.5, 5.8) ) fig.axes[0].get_xaxis().set_label_coords(0.5, -0.05) fig.subplots_adjust(bottom=0.2, right=0.95, top=0.95) fig.canvas.set_window_title(figbase + 'd Qprofiles') return fig def spikemetrics_vs_DC(neurons, a, Fdrive, Adrive, tstim, toffset, PRF, DCs, inputdir): ''' Comparison of spiking metrics for PW stimuli at various duty cycle for different neuron types. ''' metrics_dict = {} xlabel = 'Duty cycle (%)' colors = list(plt.get_cmap('Paired').colors[:6]) del colors[2:4] colors_dict = {} for i, neuron in enumerate(neurons): subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [PRF], DCs, 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [PRF], DCs, 'full')) metrics_files = {x: '{}_spikemetrics_vs_DC_{}.csv'.format(neuron, x) for x in ['full', 'sonic']} metrics_fpaths = {key: os.path.join(inputdir, value) for key, value in metrics_files.items()} sonic_fpaths = sonic_fpaths[1:] + [sonic_fpaths[0]] full_fpaths = full_fpaths[1:] + [full_fpaths[0]] data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} metrics_dict[neuron] = getSpikingMetrics( subdir, neuron, DCs * 1e2, xlabel, data_fpaths, metrics_fpaths) colors_dict[neuron] = {'full': colors[2 * i], 'sonic': colors[2 * i + 1]} fig = plotSpikingMetrics(DCs * 1e2, xlabel, metrics_dict, spikeamp=False, colors=colors_dict) fig.canvas.set_window_title(figbase + 'd spikemetrics') return fig def Qprofiles_vs_PRF(neuron, a, Fdrive, Adrive, tstim, toffset, PRFs, DC, inputdir): ''' Comparison of resulting charge profiles for PW stimuli at 5% duty cycle with different pulse repetition frequencies. ''' subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], PRFs, [DC], 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], PRFs, [DC], 'full')) patches = [False, True] * len(PRFs) patches[-1] = False comp_plot = CompTimeSeries(sum([[x, y] for x, y in zip(full_fpaths, sonic_fpaths)], []), 'Qm') fig = comp_plot.render( labels=sum([['', '{}Hz PRF'.format(si_format(PRF, space=' '))] for PRF in PRFs], []), lines=['-', '--'] * len(PRFs), colors=plt.get_cmap('Paired').colors[4:12], fs=8, patches=patches, xticks=[0, 250], yticks=[getPointNeuron(neuron).Vm0, 50], straightlegend=True, figsize=cm2inch(12.5, 5.8) ) fig.axes[0].get_xaxis().set_label_coords(0.5, -0.05) fig.subplots_adjust(bottom=0.2, right=0.95, top=0.95) fig.canvas.set_window_title(figbase + 'e Qprofiles') return fig def spikemetrics_vs_PRF(neuron, a, Fdrive, Adrive, tstim, toffset, PRFs, DC, inputdir): ''' Comparison of spiking metrics for PW stimuli at 5% duty cycle with different pulse repetition frequencies. ''' xlabel = 'PRF (Hz)' subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], PRFs, [DC], 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], PRFs, [DC], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} metrics_files = {x: '{}_spikemetrics_vs_PRF_{}.csv'.format(neuron, x) for x in ['full', 'sonic']} metrics_fpaths = {key: os.path.join(inputdir, value) for key, value in metrics_files.items()} metrics = getSpikingMetrics( subdir, neuron, PRFs, xlabel, data_fpaths, metrics_fpaths) fig = plotSpikingMetrics(PRFs, xlabel, {neuron: metrics}, spikeamp=False, logscale=True) fig.canvas.set_window_title(figbase + 'e spikemetrics') return fig def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - ap.add_argument('-i', '--inputdir', type=str, help='Input directory') - ap.add_argument('-f', '--figset', type=str, help='Figure set', default='a') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - try: - inputdir = selectDirDialog() if args.inputdir is None else args.inputdir - except ValueError as err: - logger.error(err) - return - figset = args.figset - - logger.info('Generating panel {} of {}'.format(figset, figbase)) + + parser = FigureParser(['a', 'b', 'c', 'd', 'e']) + parser.addInputDir() + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] + inputdir = args['inputdir'] + logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters radii = np.array([16, 22.6, 32, 45.3, 64]) * 1e-9 # m a = 32e-9 # m tstim = 150e-3 # s toffset = 100e-3 # s freqs = np.array([20e3, 100e3, 500e3, 1e6, 2e6, 3e6, 4e6]) # Hz Fdrive = 500e3 # Hz amps = np.array([50, 100, 300, 600]) * 1e3 # Pa Adrive = 100e3 # Pa PRFs_sparse = np.array([1e1, 1e2, 1e3, 1e4]) # Hz PRFs_dense = sum([[x, 2 * x, 5 * x] for x in PRFs_sparse[:-1]], []) + [PRFs_sparse[-1]] # Hz PRF = 100 # Hz DCs = np.array([5, 10, 25, 50, 75, 100]) * 1e-2 DC = 0.05 # Get threshold amplitudes if needed if 'a' in figset or 'b' in figset: CW_Athr_vs_Fdrive = getCWtitrations_vs_Fdrive( ['RS'], a, freqs, tstim, toffset, os.path.join(inputdir, 'CW_Athrs_vs_freqs.csv')) if 'c' in figset: CW_Athr_vs_radius = getCWtitrations_vs_radius( ['RS'], radii, Fdrive, tstim, toffset, os.path.join(inputdir, 'CW_Athrs_vs_radius.csv')) # Generate figures figs = [] - if figset == 'a': + if 'a' in figset: figs.append(Qprofiles_vs_amp('RS', a, Fdrive, CW_Athr_vs_Fdrive, tstim, toffset, inputdir)) figs.append(spikemetrics_vs_amp('RS', a, Fdrive, amps, tstim, toffset, inputdir)) - if figset == 'b': + if 'b' in figset: figs.append(Qprofiles_vs_freq( 'RS', a, [freqs.min(), freqs.max()], CW_Athr_vs_Fdrive, tstim, toffset, inputdir)) - figs.append(spikemetrics_vs_freq('RS', a, freqs, CW_Athr_vs_Fdrive, tstim, toffset, inputdir)) - if figset == 'c': + figs.append(spikemetrics_vs_freq( + 'RS', a, freqs, CW_Athr_vs_Fdrive, tstim, toffset, inputdir)) + if 'c' in figset: figs.append(Qprofiles_vs_radius( 'RS', [radii.min(), radii.max()], Fdrive, CW_Athr_vs_radius, tstim, toffset, inputdir)) figs.append(spikemetrics_vs_radius( 'RS', radii, Fdrive, CW_Athr_vs_radius, tstim, toffset, inputdir)) - if figset == 'd': + if 'd' in figset: figs.append(Qprofiles_vs_DC( ['RS', 'LTS'], a, Fdrive, Adrive, tstim, toffset, PRF, DC, inputdir)) figs.append(spikemetrics_vs_DC( ['RS', 'LTS'], a, Fdrive, Adrive, tstim, toffset, PRF, DCs, inputdir)) - if figset == 'e': + if 'e' in figset: figs.append(Qprofiles_vs_PRF( 'LTS', a, Fdrive, Adrive, tstim, toffset, PRFs_sparse, DC, inputdir)) figs.append(spikemetrics_vs_PRF( 'LTS', a, Fdrive, Adrive, tstim, toffset, PRFs_dense, DC, inputdir)) - if args.save: + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(inputdir, figname), transparent=True) + fig.savefig(os.path.join(args['outputdir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig6.py b/paper figures/fig6.py index cc2c27f..1b5de53 100644 --- a/paper figures/fig6.py +++ b/paper figures/fig6.py @@ -1,305 +1,288 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-06-06 18:38:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-14 08:08:28 +# @Last Modified time: 2019-06-17 22:01:42 ''' Sub-panels of the NICE and SONIC computation times comparative figure. ''' import os -import logging import numpy as np import matplotlib import matplotlib.pyplot as plt -from argparse import ArgumentParser from PySONIC.utils import * from PySONIC.neurons import * from PySONIC.plt import cm2inch +from PySONIC.parsers import FigureParser from utils import * # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] time_indicators = [1, 60, 60**2, 60**2 * 24, 60**2 * 24 * 7] time_indicators_labels = ['1 s', '1 min', '1 hour', '1 day', '1 week'] def comptime_vs_amp(neuron, a, Fdrive, amps, tstim, toffset, inputdir, fs=8, lw=2, ps=4): ''' Comparative plot of computation times for different acoustic amplitudes. ''' # Get filepaths xlabel = 'Amplitude (kPa)' subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], amps, [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], amps, [tstim], [toffset], [None], [1.], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} # Extract computation times (s) comptimes_fpath = os.path.join(inputdir, '{}_comptimes_vs_amps.csv'.format(neuron)) comptimes = getCompTimesQuant( inputdir, neuron, amps * 1e-3, xlabel, data_fpaths, comptimes_fpath) tcomp_lookup = getLookupsCompTime(neuron) # Extract threshold excitation amplitude CW_Athr_vs_Fdrive = getCWtitrations_vs_Fdrive( ['RS'], a, [Fdrive], tstim, toffset, os.path.join(inputdir, 'CW_Athrs_vs_freqs.csv')) Athr = CW_Athr_vs_Fdrive.loc[Fdrive * 1e-3, 'RS'] # Plot comparative profiles of computation times vs. amplitude fig, ax = plt.subplots(figsize=cm2inch(5.5, 5.8)) plt.subplots_adjust(bottom=0.2, left=0.25, right=0.95, top=0.95) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.set_xlabel(xlabel, fontsize=fs, labelpad=1) ax.set_ylabel('Computation time (s)', fontsize=fs) ax.set_xscale('log') ax.set_yscale('log') ax.set_ylim((1e-1, 1e6)) for y, lbl in zip(time_indicators, time_indicators_labels): ax.axhline(y, linewidth=0.5, linestyle='--', c='k') ax.text(amps.max() * 1e-3, 1.2 * y, lbl, horizontalalignment='right', fontsize=fs) ax.get_yaxis().set_tick_params(which='minor', size=0) ax.get_yaxis().set_tick_params(which='minor', width=0) ax.axhline(tcomp_lookup, color='k', linewidth=lw) ax.axvline(Athr, linestyle='--', color='k', linewidth=1) colors = ['silver', 'dimgrey'] for i, key in enumerate(comptimes): ax.plot(amps * 1e-3, comptimes[key], 'o--', color=colors[i], linewidth=lw, label=key, markersize=ps) for item in ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(): item.set_fontsize(fs) fig.canvas.set_window_title(figbase + 'a') return fig def comptime_vs_freq(neuron, a, freqs, CW_Athrs, tstim, toffset, inputdir, fs=8, lw=2, ps=4): ''' Comparative plot of computation times for different US frequencies. ''' # Get filepaths xlabel = 'Frequency (kHz)' subdir = os.path.join(inputdir, neuron) sonic_fpaths, full_fpaths = [], [] for Fdrive in freqs: Athr = CW_Athrs[neuron].loc[Fdrive * 1e-3] # kPa Adrive = (Athr + 20.) * 1e3 # Pa sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} # Extract computation times (s) comptimes_fpath = os.path.join(inputdir, '{}_comptimes_vs_freqs.csv'.format(neuron)) comptimes = getCompTimesQuant( inputdir, neuron, freqs * 1e-3, xlabel, data_fpaths, comptimes_fpath) tcomp_lookup = getLookupsCompTime(neuron) # Plot comparative profiles of computation time vs. frequency fig, ax = plt.subplots(figsize=cm2inch(5.5, 5.8)) plt.subplots_adjust(bottom=0.2, left=0.25, right=0.95, top=0.95) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.set_xlabel(xlabel, fontsize=fs, labelpad=1) ax.set_ylabel('Computation time (s)', fontsize=fs) ax.set_xscale('log') ax.set_yscale('log') ax.set_ylim((1e-1, 1e6)) for y, lbl in zip(time_indicators, time_indicators_labels): ax.axhline(y, linewidth=0.5, linestyle='--', c='k') ax.text(freqs.max() * 1e-3, 1.2 * y, lbl, horizontalalignment='right', fontsize=fs) ax.get_yaxis().set_tick_params(which='minor', size=0) ax.get_yaxis().set_tick_params(which='minor', width=0) ax.axhline(tcomp_lookup, color='k', linewidth=lw) colors = ['silver', 'dimgrey'] for i, key in enumerate(comptimes): ax.plot(freqs * 1e-3, comptimes[key], 'o--', color=colors[i], linewidth=lw, label=key, markersize=ps) for item in ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(): item.set_fontsize(fs) fig.canvas.set_window_title(figbase + 'b') return fig def comptime_vs_radius(neuron, radii, Fdrive, CW_Athrs, tstim, toffset, inputdir, fs=8, lw=2, ps=4): ''' Comparative plot of computation times for different sonophore radii. ''' # Get filepaths xlabel = 'Sonophore radius (nm)' subdir = os.path.join(inputdir, neuron) sonic_fpaths, full_fpaths = [], [] for a in radii: Athr = CW_Athrs[neuron].loc[np.round(a * 1e9, 1)] # kPa Adrive = (Athr + 20.) * 1e3 # Pa sonic_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'sonic')) full_fpaths += getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [None], [1.], 'full')) data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} # Extract computation times (s) comptimes_fpath = os.path.join(inputdir, '{}_comptimes_vs_radius.csv'.format(neuron)) comptimes = getCompTimesQuant( inputdir, neuron, radii * 1e9, xlabel, data_fpaths, comptimes_fpath) tcomp_lookup = getLookupsCompTime(neuron) # Plot comparative profiles of computation time vs. frequency fig, ax = plt.subplots(figsize=cm2inch(5.5, 5.8)) plt.subplots_adjust(bottom=0.2, left=0.25, right=0.95, top=0.95) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.set_xlabel(xlabel, fontsize=fs, labelpad=1) ax.set_ylabel('Computation time (s)', fontsize=fs) ax.set_xscale('log') ax.set_yscale('log') ax.set_ylim((1e-1, 1e6)) for y, lbl in zip(time_indicators, time_indicators_labels): ax.axhline(y, linewidth=0.5, linestyle='--', c='k') ax.text(radii.max() * 1e9, 1.2 * y, lbl, horizontalalignment='right', fontsize=fs) ax.get_yaxis().set_tick_params(which='minor', size=0) ax.get_yaxis().set_tick_params(which='minor', width=0) ax.axhline(tcomp_lookup, color='k', linewidth=lw) colors = ['silver', 'dimgrey'] for i, key in enumerate(comptimes): ax.plot(radii * 1e9, comptimes[key], 'o--', color=colors[i], linewidth=lw, label=key, markersize=ps) for item in ax.get_yticklabels(): item.set_fontsize(fs) for item in ax.get_xticklabels(): item.set_fontsize(fs) fig.canvas.set_window_title(figbase + 'c') return fig def comptime_vs_DC(neurons, a, Fdrive, Adrive, tstim, toffset, PRF, DCs, inputdir, fs=8, lw=2, ps=4): ''' Comparative plot of computation times for different dity cycles and neuron types. ''' xlabel = 'Duty cycle (%)' colors = list(plt.get_cmap('Paired').colors[:6]) del colors[2:4] # Create figure fig, ax = plt.subplots(figsize=cm2inch(5.5, 5.8)) plt.subplots_adjust(bottom=0.2, left=0.25, right=0.95, top=0.95) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) ax.set_xlabel(xlabel, fontsize=fs, labelpad=-7) ax.set_ylabel('Computation time (s)', fontsize=fs) ax.set_xticks([DCs.min() * 1e2, DCs.max() * 1e2]) ax.set_yscale('log') ax.set_ylim((1e-1, 1e6)) for y, lbl in zip(time_indicators, time_indicators_labels): ax.axhline(y, linewidth=0.5, linestyle='--', c='k') ax.text(DCs.max() * 1e2, 1.2 * y, lbl, horizontalalignment='right', fontsize=fs) ax.get_yaxis().set_tick_params(which='minor', size=0) ax.get_yaxis().set_tick_params(which='minor', width=0) # Loop through neurons for i, neuron in enumerate(neurons): # Get filepaths subdir = os.path.join(inputdir, neuron) sonic_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [PRF], DCs, 'sonic')) full_fpaths = getSims(subdir, neuron, a, nbls.simQueue( [Fdrive], [Adrive], [tstim], [toffset], [PRF], DCs, 'full')) sonic_fpaths = sonic_fpaths[1:] + [sonic_fpaths[0]] full_fpaths = full_fpaths[1:] + [full_fpaths[0]] data_fpaths = {'full': full_fpaths, 'sonic': sonic_fpaths} # Extract computation times (s) comptimes_fpath = os.path.join(inputdir, '{}_comptimes_vs_DC.csv'.format(neuron)) comptimes = getCompTimesQuant( inputdir, neuron, DCs * 1e2, xlabel, data_fpaths, comptimes_fpath) # Plot ax.plot(DCs * 1e2, comptimes['full'], 'o--', color=colors[2 * i], linewidth=lw, markersize=ps) ax.plot(DCs * 1e2, comptimes['sonic'], 'o--', color=colors[2 * i + 1], linewidth=lw, markersize=ps, label=neuron) fig.canvas.set_window_title(figbase + 'd') return fig def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - ap.add_argument('-i', '--inputdir', type=str, help='Input directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - try: - inputdir = selectDirDialog() if args.inputdir is None else args.inputdir - except ValueError as err: - logger.error(err) - return - figset = args.figset - if figset == 'all': - figset = ['a', 'b', 'c', 'd'] - + parser = FigureParser(['a', 'b', 'c', 'd']) + parser.addInputDir() + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] + inputdir = args['inputdir'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters a = 32e-9 # m radii = np.array([16, 22.6, 32, 45.3, 64]) * 1e-9 # nm tstim = 150e-3 # s toffset = 100e-3 # s freqs = np.array([20e3, 100e3, 500e3, 1e6, 2e6, 3e6, 4e6]) # Hz Fdrive = 500e3 # Hz CW_Athr_vs_Fdrive = getCWtitrations_vs_Fdrive( ['RS'], a, freqs, tstim, toffset, os.path.join(inputdir, 'CW_Athrs_vs_freqs.csv')) Athr = CW_Athr_vs_Fdrive['RS'].loc[Fdrive * 1e-3] amps1 = np.array([Athr - 5, Athr, Athr + 20]) * 1e3 amps2 = np.array([50, 100, 300, 600]) * 1e3 # Pa amps = np.sort(np.hstack([amps1, amps2])) CW_Athr_vs_radius = getCWtitrations_vs_radius( ['RS'], radii, Fdrive, tstim, toffset, os.path.join(inputdir, 'CW_Athrs_vs_radius.csv')) Adrive = 100e3 # Pa PRF = 100 # Hz DCs = np.array([5, 10, 25, 50, 75, 100]) * 1e-2 # Generate figures figs = [] if 'a' in figset: figs.append(comptime_vs_amp('RS', a, Fdrive, amps, tstim, toffset, inputdir)) if 'b' in figset: figs.append(comptime_vs_freq('RS', a, freqs, CW_Athr_vs_Fdrive, tstim, toffset, inputdir)) if 'c' in figset: figs.append(comptime_vs_radius( 'RS', radii, Fdrive, CW_Athr_vs_radius, tstim, toffset, inputdir)) if 'd' in figset: figs.append(comptime_vs_DC( ['RS', 'LTS'], a, Fdrive, Adrive, tstim, toffset, PRF, DCs, inputdir)) - if args.save: + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(inputdir, figname), transparent=True) + fig.savefig(os.path.join(args['outpudir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig7.py b/paper figures/fig7.py index 4bfd2a9..12cff4f 100644 --- a/paper figures/fig7.py +++ b/paper figures/fig7.py @@ -1,151 +1,133 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-09-26 09:51:43 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-17 09:26:49 +# @Last Modified time: 2019-06-17 22:03:38 ''' Sub-panels of (duty-cycle x amplitude) US activation maps and related Q-V traces. ''' import os import numpy as np -import logging import matplotlib import matplotlib.pyplot as plt -from argparse import ArgumentParser from PySONIC.core import NeuronalBilayerSonophore -from PySONIC.utils import logger, selectDirDialog, si_format +from PySONIC.utils import logger, si_format from PySONIC.plt import ActivationMap from PySONIC.neurons import getPointNeuron +from PySONIC.parsers import FigureParser # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def plotMapAndTraces(inputdir, pneuron, a, Fdrive, tstim, amps, PRF, DCs, FRbounds, insets, tbounds, Vbounds, prefix): # Activation map mapcode = '{} {}Hz PRF{}Hz 1s'.format(pneuron.name, *si_format([Fdrive, PRF, tstim], space='')) subdir = os.path.join(inputdir, mapcode) actmap = ActivationMap(subdir, pneuron, a, Fdrive, tstim, PRF, amps, DCs) mapfig = actmap.render(FRbounds=FRbounds, thresholds=True) mapfig.canvas.set_window_title('{} map {}'.format(prefix, mapcode)) ax = mapfig.axes[0] DC_insets, A_insets = zip(*insets) ax.scatter(DC_insets, A_insets, s=80, facecolors='none', edgecolors='k', linestyle='--') # Related inset traces tracefigs = [] nbls = NeuronalBilayerSonophore(a, pneuron) for inset in insets: DC = inset[0] * 1e-2 Adrive = inset[1] * 1e3 fname = '{}.pkl'.format(nbls.filecode( Fdrive, actmap.correctAmp(Adrive), tstim, 0., PRF, DC, 'sonic')) fpath = os.path.join(subdir, fname) tracefig = actmap.plotQVeff(fpath, tbounds=tbounds, ybounds=Vbounds) figcode = '{} VQ trace {} {:.1f}kPa {:.0f}%DC'.format( prefix, pneuron.name, Adrive * 1e-3, DC * 1e2) tracefig.canvas.set_window_title(figcode) tracefigs.append(tracefig) return mapfig, tracefigs def panel(inputdir, pneurons, a, tstim, PRF, amps, DCs, FRbounds, tbounds, Vbounds, insets, prefix): mapfigs, tracefigs = [], [] for pn in pneurons: out = plotMapAndTraces( inputdir, pn, a, 500e3, tstim, amps, PRF, DCs, FRbounds, insets[pn.name], tbounds, Vbounds, prefix) mapfigs.append(out[0]) tracefigs += out[1] return mapfigs + tracefigs def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - ap.add_argument('-i', '--inputdir', type=str, help='Input directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - try: - inputdir = selectDirDialog() if args.inputdir is None else args.inputdir - except ValueError as err: - logger.error(err) - return - figset = args.figset - if figset == 'all': - figset = ['a', 'b', 'c'] - - logger.info('Generating panel {} of {}'.format(figset, figbase)) + parser = FigureParser(['a', 'b', 'c']) + parser.addInputDir() + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] + inputdir = args['inputdir'] + logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters pneurons = [getPointNeuron(n) for n in ['RS', 'LTS']] a = 32e-9 # m tstim = 1.0 # s amps = np.logspace(np.log10(10), np.log10(600), num=30) * 1e3 # Pa DCs = np.arange(1, 101) * 1e-2 FRbounds = (1e0, 1e3) # Hz tbounds = (0, 240e-3) # s Vbounds = -150, 50 # mV # Generate figures try: - figs = [] if 'a' in figset: PRF = 1e1 insets = { 'RS': [(28, 127.0), (37, 168.4)], 'LTS': [(8, 47.3), (30, 146.2)] } figs += panel(inputdir, pneurons, a, tstim, PRF, amps, DCs, FRbounds, tbounds, Vbounds, insets, figbase + 'a') if 'b' in figset: PRF = 1e2 insets = { 'RS': [(51, 452.4), (56, 452.4)], 'LTS': [(13, 193.9), (43, 257.2)] } figs += panel(inputdir, pneurons, a, tstim, PRF, amps, DCs, FRbounds, tbounds, Vbounds, insets, figbase + 'b') if 'c' in figset: PRF = 1e3 insets = { 'RS': [(40, 110.2), (64, 193.9)], 'LTS': [(10, 47.3), (53, 168.4)] } figs += panel(inputdir, pneurons, a, tstim, PRF, amps, DCs, FRbounds, tbounds, Vbounds, insets, figbase + 'c') except Exception as e: logger.error(e) - quit() + return - if args.save: + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(inputdir, figname), transparent=True) + fig.savefig(os.path.join(args['outputdir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig8.py b/paper figures/fig8.py index 580bc09..09f3515 100644 --- a/paper figures/fig8.py +++ b/paper figures/fig8.py @@ -1,152 +1,134 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-11-27 17:57:45 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-14 08:09:16 +# @Last Modified time: 2019-06-17 22:05:08 ''' Sub-panels of threshold curves for various sonophore radii and US frequencies. ''' import os -import logging import numpy as np import pandas as pd import matplotlib import matplotlib.pyplot as plt -from argparse import ArgumentParser -from PySONIC.neurons import getPointNeuron -from PySONIC.utils import logger, si_format, selectDirDialog +from PySONIC.utils import logger, si_format from PySONIC.plt import cm2inch +from PySONIC.parsers import FigureParser # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def getThresholdAmplitudes(root, neuron, a, Fdrive, tstim, PRF): subfolder = '{} {:.0f}nm {}Hz PRF{}Hz {}s'.format( neuron, a * 1e9, *si_format([Fdrive, PRF, tstim], 0, space='') ) fname = 'log_ASTIM.xlsx' fpath = os.path.join(root, subfolder, fname) df = pd.read_excel(fpath, sheet_name='Data') DCs = df['Duty factor'].values Athrs = df['Adrive (kPa)'].values iDCs = np.argsort(DCs) DCs = DCs[iDCs] Athrs = Athrs[iDCs] return DCs, Athrs def plotThresholdAmps(root, neurons, radii, freqs, PRF, tstim, fs=10, colors=None, figsize=None): ''' Plot threshold excitation amplitudes of several neurons determined by titration procedures, as a function of duty cycle, for various combinations of sonophore radius and US frequency. :param neurons: list of neuron names :param radii: list of sonophore radii (m) :param freqs: list US frequencies (Hz) :param PRF: pulse repetition frequency used for titration procedures (Hz) :param tstim: stimulus duration used for titration procedures :return: figure handle ''' if figsize is None: figsize = cm2inch(8, 7) linestyles = ['--', ':', '-.'] assert len(freqs) <= len(linestyles), 'too many frequencies' fig, ax = plt.subplots(figsize=figsize) ax.set_xlabel('Duty cycle (%)', fontsize=fs) ax.set_ylabel('Amplitude (kPa)', fontsize=fs) for item in ax.get_xticklabels() + ax.get_yticklabels(): item.set_fontsize(fs) ax.set_yscale('log') ax.set_xlim([0, 100]) ax.set_ylim([10, 600]) linestyles = ['-', '--'] for neuron, ls in zip(neurons, linestyles): icolor = 0 for i, a in enumerate(radii): for j, Fdrive in enumerate(freqs): if colors is None: color = 'C{}'.format(icolor) else: color = colors[icolor] DCs, Athrs = getThresholdAmplitudes(root, neuron, a, Fdrive, tstim, PRF) lbl = '{} neuron, {:.0f} nm, {}Hz, {}Hz PRF'.format( neuron, a * 1e9, *si_format([Fdrive, PRF], 0, space=' ')) ax.plot(DCs * 1e2, Athrs, ls, c=color, label=lbl) icolor += 1 ax.legend(fontsize=fs - 5, frameon=False) fig.tight_layout() return fig def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - ap.add_argument('-i', '--inputdir', type=str, help='Input directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - try: - inputdir = selectDirDialog() if args.inputdir is None else args.inputdir - except ValueError as err: - logger.error(err) - return - figset = args.figset - if figset == 'all': - figset = ['a', 'b'] - + parser = FigureParser(['a', 'b']) + parser.addInputDir() + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] + inputdir = args['inputdir'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters neurons = ['RS', 'LTS'] radii = np.array([16, 32, 64]) * 1e-9 # m a = radii[1] freqs = np.array([20, 500, 4000]) * 1e3 # Hz Fdrive = freqs[1] PRFs = np.array([1e1, 1e2, 1e3]) # Hz PRF = PRFs[1] tstim = 1 # s colors = plt.get_cmap('tab20c').colors # Generate figures figs = [] if 'a' in figset: fig = plotThresholdAmps(inputdir, neurons, radii, [Fdrive], PRF, tstim, fs=12, colors=colors[:3][::-1]) fig.canvas.set_window_title(figbase + 'a') figs.append(fig) if 'b' in figset: fig = plotThresholdAmps(inputdir, neurons, [a], freqs, PRF, tstim, fs=12, colors=colors[8:11][::-1]) fig.canvas.set_window_title(figbase + 'b') figs.append(fig) - if args.save: + if args['save']: for fig in figs: figname = '{}.pdf'.format(fig.canvas.get_window_title()) - fig.savefig(os.path.join(inputdir, figname), transparent=True) + fig.savefig(os.path.join(args['outpudir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/paper figures/fig9.py b/paper figures/fig9.py index 3065f40..2cf203c 100644 --- a/paper figures/fig9.py +++ b/paper figures/fig9.py @@ -1,111 +1,94 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-12-09 12:06:01 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-17 14:24:31 +# @Last Modified time: 2019-06-17 22:06:53 ''' Sub-panels of SONIC model validation on an STN neuron (response to CW sonication). ''' import os -import logging import numpy as np import matplotlib import matplotlib.pyplot as plt -from argparse import ArgumentParser from PySONIC.core import NeuronalBilayerSonophore from PySONIC.neurons import getPointNeuron -from PySONIC.utils import logger, selectDirDialog, Intensity2Pressure +from PySONIC.utils import logger, Intensity2Pressure from PySONIC.plt import CompTimeSeries, GroupedTimeSeries +from PySONIC.parsers import FigureParser # Plot parameters matplotlib.rcParams['pdf.fonttype'] = 42 matplotlib.rcParams['ps.fonttype'] = 42 matplotlib.rcParams['font.family'] = 'arial' # Figure basename figbase = os.path.splitext(__file__)[0] def main(): - ap = ArgumentParser() - - # Runtime options - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - ap.add_argument('-i', '--inputdir', type=str, help='Input directory') - ap.add_argument('-f', '--figset', type=str, nargs='+', help='Figure set', default='all') - ap.add_argument('-s', '--save', default=False, action='store_true', - help='Save output figures as pdf') - - args = ap.parse_args() - loglevel = logging.DEBUG if args.verbose is True else logging.INFO - logger.setLevel(loglevel) - try: - inputdir = selectDirDialog() if args.inputdir is None else args.inputdir - except ValueError as err: - logger.error(err) - return - figset = args.figset - if figset is 'all': - figset = ['a', 'b'] - + parser = FigureParser(['a', 'b']) + parser.addInputDir() + args = parser.parse() + logger.setLevel(args['loglevel']) + figset = args['subset'] + inputdir = args['inputdir'] logger.info('Generating panels {} of {}'.format(figset, figbase)) # Parameters pneuron = getPointNeuron('STN') a = 32e-9 # m Fdrive = 500e3 # Hz tstim = 1 # s toffset = 0. # s PRF = 1e2 DC = 1. nbls = NeuronalBilayerSonophore(a, pneuron) # Range of intensities intensities = pneuron.getLowIntensities() # W/m2 # Levels depicted with individual traces subset_intensities = [112, 114, 123] # W/m2 # convert to amplitudes and get filepaths amplitudes = Intensity2Pressure(intensities) # Pa fnames = ['{}.pkl'.format(nbls.filecode(Fdrive, A, tstim, toffset, PRF, DC, 'sonic')) for A in amplitudes] fpaths = [os.path.join(inputdir, 'STN', fn) for fn in fnames] # Generate figures figs = [] if 'a' in figset: comp_plot = CompTimeSeries(fpaths, 'FR') fig = comp_plot.render( patches='none', cmap='Oranges', tbounds=(0, tstim) ) fig.canvas.set_window_title(figbase + 'a') figs.append(fig) if 'b' in figset: isubset = [np.argwhere(intensities == x)[0][0] for x in subset_intensities] subset_amplitudes = amplitudes[isubset] titles = ['{:.2f} kPa ({:.0f} W/m2)'.format(A * 1e-3, I) for A, I in zip(subset_amplitudes, subset_intensities)] figtraces = GroupedTimeSeries([fpaths[i] for i in isubset], pltscheme={'Q_m': ['Qm']})() for fig, title in zip(figtraces, titles): fig.axes[0].set_title(title) fig.canvas.set_window_title(figbase + 'b {}'.format(title)) figs.append(fig) - if args.save: + if args['save']: for fig in figs: s = fig.canvas.get_window_title() s = s.replace('(', '- ').replace('/', '_').replace(')', '') figname = '{}.pdf'.format(s) - fig.savefig(os.path.join(inputdir, figname), transparent=True) + fig.savefig(os.path.join(args['outputdir'], figname), transparent=True) else: plt.show() if __name__ == '__main__': main() diff --git a/tests/test_basic.py b/tests/test_basic.py index 4a39ec3..d4acffe 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,186 +1,184 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-06-14 18:37:45 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-17 21:37:43 +# @Last Modified time: 2019-06-17 21:51:05 ''' Test the basic functionalities of the package. ''' import os import time import cProfile import pstats from PySONIC.core import BilayerSonophore, NeuronalBilayerSonophore from PySONIC.utils import logger from PySONIC.neurons import getPointNeuron from PySONIC.parsers import TestParser def execute(func_str, globals, locals, is_profiled): ''' Execute function with or without profiling. ''' if is_profiled: pfile = 'tmp.stats' cProfile.runctx(func_str, globals, locals, pfile) stats = pstats.Stats(pfile) os.remove(pfile) stats.strip_dirs() stats.sort_stats('cumulative') stats.print_stats() else: eval(func_str, globals, locals) def test_MECH(is_profiled=False): ''' Mechanical simulation. ''' logger.info('Test: running MECH simulation') # Create BLS instance a = 32e-9 # m Qm0 = -80e-5 # membrane resting charge density (C/m2) Cm0 = 1e-2 # membrane resting capacitance (F/m2) bls = BilayerSonophore(a, Cm0, Qm0) # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa Qm = 50e-5 # C/m2 # Run simulation execute('bls.simulate(Fdrive, Adrive, Qm)', globals(), locals(), is_profiled) def test_ESTIM(is_profiled=False): ''' Electrical simulation ''' logger.info('Test: running ESTIM simulation') # Initialize neuron pneuron = getPointNeuron('RS') # Stimulation parameters Astim = 10.0 # mA/m2 tstim = 100e-3 # s toffset = 50e-3 # s # Run simulation execute('pneuron.simulate(Astim, tstim, toffset)', globals(), locals(), is_profiled) def test_ASTIM_sonic(is_profiled=False): ''' Effective acoustic simulation ''' logger.info('Test: ASTIM sonic simulation') # Default parameters a = 32e-9 # m pneuron = getPointNeuron('RS') nbls = NeuronalBilayerSonophore(a, pneuron) Fdrive = 500e3 # Hz Adrive = 100e3 # Pa tstim = 50e-3 # s toffset = 10e-3 # s # test error 1: sonophore radius outside of lookup range try: nbls = NeuronalBilayerSonophore(100e-9, pneuron) nbls.simulate(Fdrive, Adrive, tstim, toffset, method='sonic') except ValueError: logger.debug('Out of range radius: OK') # test error 2: frequency outside of lookups range try: nbls = NeuronalBilayerSonophore(a, pneuron) nbls.simulate(10e3, Adrive, tstim, toffset, method='sonic') except ValueError: logger.debug('Out of range frequency: OK') # test error 3: amplitude outside of lookups range try: nbls = NeuronalBilayerSonophore(a, pneuron) nbls.simulate(Fdrive, 1e6, tstim, toffset, method='sonic') except ValueError: logger.debug('Out of range amplitude: OK') # Run simulation execute("nbls.simulate(Fdrive, Adrive, tstim, toffset, method='sonic')", globals(), locals(), is_profiled) def test_ASTIM_dense(is_profiled=False): ''' Classic acoustic simulation ''' logger.info('Test: running ASTIM detailed simulation') # Initialize sonic neuron a = 32e-9 # m pneuron = getPointNeuron('RS') nbls = NeuronalBilayerSonophore(a, pneuron) # Stimulation parameters Fdrive = 500e3 # Hz Adrive = 100e3 # Pa tstim = 1e-6 # s toffset = 1e-6 # s # Run simulation execute("nbls.simulate(Fdrive, Adrive, tstim, toffset, method='full')", globals(), locals(), is_profiled) def test_ASTIM_hybrid(is_profiled=False): ''' Hybrid acoustic simulation ''' logger.info('Test: running ASTIM hybrid simulation') # Initialize sonic neuron a = 32e-9 # m pneuron = getPointNeuron('RS') nbls = NeuronalBilayerSonophore(a, pneuron) # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa tstim = 1e-3 # s toffset = 1e-3 # s # Run simulation execute("nbls.simulate(Fdrive, Adrive, tstim, toffset, method='hybrid')", globals(), locals(), is_profiled) def main(): test_funcs = { 'MECH': test_MECH, 'ESTIM': test_ESTIM, 'ASTIM_sonic': test_ASTIM_sonic, 'ASTIM_dense': test_ASTIM_dense, 'ASTIM_hybrid': test_ASTIM_hybrid } parser = TestParser(list(test_funcs.keys())) args = parser.parse() logger.setLevel(args['loglevel']) # Parse arguments if args['profile'] and args['subset'] == 'all': logger.error('profiling can only be run on individual tests') return + print(args['subset']) + # Run test - if args['subset'] == ['all']: - t0 = time.time() - for k, test_func in test_funcs.items(): - test_func(args['profile']) - tcomp = time.time() - t0 - logger.info('All tests completed in %.0f s', tcomp) - else: - for s in args['subset']: - test_funcs[s](args['profile']) + t0 = time.time() + for s in args['subset']: + test_funcs[s](args['profile']) + tcomp = time.time() - t0 + logger.info('tests completed in %.0f s', tcomp) if __name__ == '__main__': main() diff --git a/tests/test_values.py b/tests/test_values.py deleted file mode 100644 index 0ca5dc1..0000000 --- a/tests/test_values.py +++ /dev/null @@ -1,213 +0,0 @@ -# -*- coding: utf-8 -*- -# @Author: Theo Lemaire -# @Email: theo.lemaire@epfl.ch -# @Date: 2017-06-14 18:37:45 -# @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-13 23:26:05 - -''' Run functionalities of the package and test validity of outputs. ''' - -import sys -import logging -from argparse import ArgumentParser -import numpy as np - -from PySONIC.utils import logger -from PySONIC.core import BilayerSonophore, NeuronalBilayerSonophore -from PySONIC.neurons import getNeuronsDict - -from PySONIC.constants import * - -# Set logging level -logger.setLevel(logging.INFO) - - -def test_MECH(): - ''' Maximal negative and positive deflections of the BLS structure for a specific - sonophore size, resting membrane properties and stimulation parameters. ''' - - logger.info('Starting test: Mechanical simulation') - - # Create BLS instance - a = 32e-9 # m - Cm0 = 1e-2 # membrane resting capacitance (F/m2) - Qm0 = -80e-5 # membrane resting charge density (C/m2) - bls = BilayerSonophore(a, Cm0, Qm0) - - # Run mechanical simulation - Fdrive = 350e3 # Hz - Adrive = 100e3 # Pa - Qm = 50e-5 # C/m2 - data, _ = bls.simulate(Fdrive, Adrive, Qm) - - # Check validity of deflection extrema - Zlast = data.loc[-NPC_DENSE:, 'Z'].values - Zmin, Zmax = (Zlast.min(), Zlast.max()) - logger.info('Zmin = %.2f nm, Zmax = %.2f nm', Zmin * 1e9, Zmax * 1e9) - Zmin_ref, Zmax_ref = (-0.116e-9, 5.741e-9) - assert np.abs(Zmin - Zmin_ref) < 1e-12, 'Unexpected sonophore compression amplitude' - assert np.abs(Zmax - Zmax_ref) < 1e-12, 'Unexpected sonophore expansion amplitude' - - logger.info('Passed test: Mechanical simulation') - - -def test_resting_potential(): - ''' Neurons membrane potential in free conditions should stabilize to their - specified resting potential value. ''' - - conv_err_msg = ('{} neuron membrane potential in free conditions does not converge to ' - 'stable value (gap after 20s: {:.2e} mV)') - value_err_msg = ('{} neuron steady-state membrane potential in free conditions differs ' - 'significantly from specified resting potential (gap = {:.2f} mV)') - - logger.info('Starting test: neurons resting potential') - - for Neuron in getNeuronsDict().values(): - - # Simulate each neuron in free conditions - pneuron = Neuron() - - logger.info('%s neuron simulation in free conditions', pneuron.name) - - data, _ = pneuron.simulate(Astim=0.0, tstim=20.0, toffset=0.0) - Vm_free = data['Vm'].values - - # Check membrane potential convergence - Vm_free_last, Vm_free_beforelast = (Vm_free[-1], Vm_free[-2]) - Vm_free_conv = Vm_free_last - Vm_free_beforelast - assert np.abs(Vm_free_conv) < 1e-5, conv_err_msg.format(pneuron.name, Vm_free_conv) - - # Check membrane potential convergence to resting potential - Vm_free_diff = Vm_free_last - pneuron.Vm0 - assert np.abs(Vm_free_diff) < 0.1, value_err_msg.format(pneuron.name, Vm_free_diff) - - logger.info('Passed test: neurons resting potential') - - -def test_ESTIM(): - ''' Threshold E-STIM amplitude and needed to obtain an action potential and response latency - should match reference values. ''' - - Athr_err_msg = ('{} neuron threshold amplitude for excitation does not match reference value' - '(gap = {:.2f} mA/m2)') - latency_err_msg = ('{} neuron latency for excitation at threshold amplitude does not match ' - 'reference value (gap = {:.2f} ms)') - - logger.info('Starting test: E-STIM titration') - - # Stimulation parameters - tstim = 100e-3 # s - toffset = 50e-3 # s - - # Reference values - Athr_refs = {'FS': 6.91, 'LTS': 1.54, 'RS': 5.03, 'RE': 3.61, 'TC': 4.05, - 'LeechT': 4.66, 'LeechP': 13.72, 'IB': 3.08} - - for Neuron in getNeuronsDict().values(): - - # Perform titration for each neuron - pneuron = Neuron() - logger.info('%s neuron titration', pneuron.name) - Athr = pneuron.titrate(tstim, toffset) - - # Check threshold amplitude - Athr_diff = Athr - Athr_refs[pneuron.name] - assert np.abs(Athr_diff) < 0.1, Athr_err_msg.format(pneuron.name, Athr_diff) - - logger.info('Passed test: E-STIM titration') - - -def test_ASTIM(): - ''' Threshold A-STIM amplitude and needed to obtain an action potential and response latency - should match reference values. ''' - - Athr_err_msg = ('{} neuron threshold amplitude for excitation does not match reference value' - '(gap = {:.2f} kPa)') - latency_err_msg = ('{} neuron latency for excitation at threshold amplitude does not match ' - 'reference value (gap = {:.2f} ms)') - - logger.info('Starting test: A-STIM titration') - - # Sonophore radius - a = 32e-9 # m - - # Stimulation parameters - Fdrive = 350e3 # Hz - tstim = 50e-3 # s - toffset = 30e-3 # s - - # Reference values - Athr_refs = {'FS': 38.96e3, 'LTS': 24.90e3, 'RS': 50.90e3, 'RE': 46.36e3, 'TC': 23.14e3, - 'LeechT': 21.02e3, 'LeechP': 22.23e3, 'IB': 91.26e3} - - # Titration for each neuron - for Neuron in getNeuronsDict().values(): - - # Initialize sonic neuron - pneuron = Neuron() - nbls = NeuronalBilayerSonophore(a, pneuron) - logger.info('%s neuron titration', pneuron.name) - - # Perform titration - Athr = nbls.titrate(Fdrive, tstim, toffset, method='sonic') - - # Check threshold amplitude - Athr_diff = (Athr - Athr_refs[pneuron.name]) * 1e-3 - assert np.abs(Athr_diff) < 0.1, Athr_err_msg.format(pneuron.name, Athr_diff) - - logger.info('Passed test: A-STIM titration') - - -def test_all(): - logger.info('Starting tests') - test_MECH() - test_resting_potential() - test_ESTIM() - test_ASTIM() - logger.info('All tests successfully passed') - - - -def main(): - - # Define valid test sets - valid_testsets = [ - 'MECH', - 'resting_potential', - 'ESTIM', - 'ASTIM', - 'all' - ] - - # Define argument parser - ap = ArgumentParser() - - ap.add_argument('-t', '--testset', type=str, default='all', choices=valid_testsets, - help='Specific test set') - ap.add_argument('-v', '--verbose', default=False, action='store_true', - help='Increase verbosity') - - # Parse arguments - args = ap.parse_args() - if args.verbose: - logger.setLevel(logging.DEBUG) - else: - logger.setLevel(logging.INFO) - - # Run test - try: - if args.testset == 'all': - test_all() - else: - possibles = globals().copy() - possibles.update(locals()) - method = possibles.get('test_{}'.format(args.testset)) - method() - sys.exit(0) - except AssertionError as e: - logger.error(e) - sys.exit(1) - - -if __name__ == '__main__': - main()