diff --git a/PySONIC/neurons/cortical.py b/PySONIC/neurons/cortical.py index bfd7e25..4aded4e 100644 --- a/PySONIC/neurons/cortical.py +++ b/PySONIC/neurons/cortical.py @@ -1,390 +1,390 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-07-31 15:19:51 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-06-29 18:41:50 +# @Last Modified time: 2019-09-13 14:26:14 import numpy as np from ..core import PointNeuron class Cortical(PointNeuron): ''' Generic cortical neuron Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 1e-2 # Membrane capacitance (F/m2) # Reversal potentials (mV) ENa = 50.0 # Sodium EK = -90.0 # Potassium ECa = 120.0 # Calcium # ------------------------------ Gating states kinetics ------------------------------ @classmethod def alpham(cls, Vm): return 0.32 * cls.vtrap(13 - (Vm - cls.VT), 4) * 1e3 # s-1 @classmethod def betam(cls, Vm): return 0.28 * cls.vtrap((Vm - cls.VT) - 40, 5) * 1e3 # s-1 @classmethod def alphah(cls, Vm): return 0.128 * np.exp(-((Vm - cls.VT) - 17) / 18) * 1e3 # s-1 @classmethod def betah(cls, Vm): return 4 / (1 + np.exp(-((Vm - cls.VT) - 40) / 5)) * 1e3 # s-1 @classmethod def alphan(cls, Vm): return 0.032 * cls.vtrap(15 - (Vm - cls.VT), 5) * 1e3 # s-1 @classmethod def betan(cls, Vm): return 0.5 * np.exp(-((Vm - cls.VT) - 10) / 40) * 1e3 # s-1 @staticmethod def pinf(Vm): return 1.0 / (1 + np.exp(-(Vm + 35) / 10)) @classmethod def taup(cls, Vm): return cls.TauMax / (3.3 * np.exp((Vm + 35) / 20) + np.exp(-(Vm + 35) / 20)) # s # ------------------------------ States derivatives ------------------------------ @classmethod def derStates(cls): return { 'm': lambda Vm, x: cls.alpham(Vm) * (1 - x['m']) - cls.betam(Vm) * x['m'], 'h': lambda Vm, x: cls.alphah(Vm) * (1 - x['h']) - cls.betah(Vm) * x['h'], 'n': lambda Vm, x: cls.alphan(Vm) * (1 - x['n']) - cls.betan(Vm) * x['n'], 'p': lambda Vm, x: (cls.pinf(Vm) - x['p']) / cls.taup(Vm) } # ------------------------------ Steady states ------------------------------ @classmethod def steadyStates(cls): return { 'm': lambda Vm: cls.alpham(Vm) / (cls.alpham(Vm) + cls.betam(Vm)), 'h': lambda Vm: cls.alphah(Vm) / (cls.alphah(Vm) + cls.betah(Vm)), 'n': lambda Vm: cls.alphan(Vm) / (cls.alphan(Vm) + cls.betan(Vm)), 'p': lambda Vm: cls.pinf(Vm) } # ------------------------------ Membrane currents ------------------------------ @classmethod def iNa(cls, m, h, Vm): ''' Sodium current ''' return cls.gNabar * m**3 * h * (Vm - cls.ENa) # mA/m2 @classmethod def iKd(cls, n, Vm): ''' delayed-rectifier Potassium current ''' return cls.gKdbar * n**4 * (Vm - cls.EK) # mA/m2 @classmethod def iM(cls, p, Vm): ''' slow non-inactivating Potassium current ''' return cls.gMbar * p * (Vm - cls.EK) # mA/m2 @classmethod def iLeak(cls, Vm): ''' non-specific leakage current ''' return cls.gLeak * (Vm - cls.ELeak) # mA/m2 @classmethod def currents(cls): return { 'iNa': lambda Vm, x: cls.iNa(x['m'], x['h'], Vm), 'iKd': lambda Vm, x: cls.iKd(x['n'], Vm), 'iM': lambda Vm, x: cls.iM(x['p'], Vm), 'iLeak': lambda Vm, _: cls.iLeak(Vm) } class CorticalRS(Cortical): ''' Cortical regular spiking neuron Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Neuron name name = 'RS' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Vm0 = -71.9 # Membrane potential (mV) # Reversal potentials (mV) ELeak = -70.3 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 560.0 # Sodium gKdbar = 60.0 # Delayed-rectifier Potassium gMbar = 0.75 # Slow non-inactivating Potassium gLeak = 0.205 # Non-specific leakage # Additional parameters VT = -56.2 # Spike threshold adjustment parameter (mV) TauMax = 0.608 # Max. adaptation decay of slow non-inactivating Potassium current (s) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 'p': 'iM gate' } class CorticalFS(Cortical): ''' Cortical fast-spiking neuron Reference: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* ''' # Neuron name name = 'FS' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Vm0 = -71.4 # Membrane potential (mV) # Reversal potentials (mV) ELeak = -70.4 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 580.0 # Sodium gKdbar = 39.0 # Delayed-rectifier Potassium gMbar = 0.787 # Slow non-inactivating Potassium gLeak = 0.38 # Non-specific leakage # Additional parameters VT = -57.9 # Spike threshold adjustment parameter (mV) TauMax = 0.502 # Max. adaptation decay of slow non-inactivating Potassium current (s) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 'p': 'iM gate' } class CorticalLTS(Cortical): ''' Cortical low-threshold spiking neuron References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Huguenard, J.R., and McCormick, D.A. (1992). Simulation of the currents involved in rhythmic oscillations in thalamic relay neurons. J. Neurophysiol. 68, 1373–1383.* ''' # Neuron name name = 'LTS' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Vm0 = -54.0 # Membrane potential (mV) # Reversal potentials (mV) ELeak = -50.0 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 500.0 # Sodium gKdbar = 40.0 # Delayed-rectifier Potassium gMbar = 0.28 # Slow non-inactivating Potassium gCaTbar = 4.0 # Low-threshold Calcium gLeak = 0.19 # Non-specific leakage # Additional parameters VT = -50.0 # Spike threshold adjustment parameter (mV) TauMax = 4.0 # Max. adaptation decay of slow non-inactivating Potassium current (s) Vx = -7.0 # Voltage-dependence uniform shift factor at 36°C (mV) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 'p': 'iM gate', 's': 'iCaT activation gate', 'u': 'iCaT inactivation gate' } # ------------------------------ Gating states kinetics ------------------------------ @classmethod def sinf(cls, Vm): return 1.0 / (1.0 + np.exp(-(Vm + cls.Vx + 57.0) / 6.2)) @classmethod def taus(cls, Vm): x = np.exp(-(Vm + cls.Vx + 132.0) / 16.7) + np.exp((Vm + cls.Vx + 16.8) / 18.2) return 1.0 / 3.7 * (0.612 + 1.0 / x) * 1e-3 # s @classmethod def uinf(cls, Vm): return 1.0 / (1.0 + np.exp((Vm + cls.Vx + 81.0) / 4.0)) @classmethod def tauu(cls, Vm): if Vm + cls.Vx < -80.0: return 1.0 / 3.7 * np.exp((Vm + cls.Vx + 467.0) / 66.6) * 1e-3 # s else: return 1.0 / 3.7 * (np.exp(-(Vm + cls.Vx + 22) / 10.5) + 28.0) * 1e-3 # s # ------------------------------ States derivatives ------------------------------ @classmethod def derStates(cls): return {**super().derStates(), **{ 's': lambda Vm, x: (cls.sinf(Vm) - x['s']) / cls.taus(Vm), 'u': lambda Vm, x: (cls.uinf(Vm) - x['u']) / cls.tauu(Vm) }} # ------------------------------ Steady states ------------------------------ @classmethod def steadyStates(cls): return {**super().steadyStates(), **{ 's': lambda Vm: cls.sinf(Vm), 'u': lambda Vm: cls.uinf(Vm) }} # ------------------------------ Membrane currents ------------------------------ @classmethod def iCaT(cls, s, u, Vm): ''' low-threshold (T-type) Calcium current ''' return cls.gCaTbar * s**2 * u * (Vm - cls.ECa) # mA/m2 @classmethod def currents(cls): return {**super().currents(), **{ 'iCaT': lambda Vm, x: cls.iCaT(x['s'], x['u'], Vm) }} class CorticalIB(Cortical): ''' Cortical intrinsically bursting neuron References: *Pospischil, M., Toledo-Rodriguez, M., Monier, C., Piwkowska, Z., Bal, T., Frégnac, Y., Markram, H., and Destexhe, A. (2008). Minimal Hodgkin-Huxley type models for different classes of cortical and thalamic neurons. Biol Cybern 99, 427–441.* *Reuveni, I., Friedman, A., Amitai, Y., and Gutnick, M.J. (1993). Stepwise repolarization from Ca2+ plateaus in neocortical pyramidal cells: evidence for nonhomogeneous distribution of HVA Ca2+ channels in dendrites. J. Neurosci. 13, 4609–4621.* ''' # Neuron name name = 'IB' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Vm0 = -71.4 # Membrane potential (mV) # Reversal potentials (mV) ELeak = -70 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 500 # Sodium gKdbar = 50 # Delayed-rectifier Potassium gMbar = 0.3 # Slow non-inactivating Potassium gCaLbar = 1.0 # High-threshold Calcium gLeak = 0.1 # Non-specific leakage # Additional parameters VT = -56.2 # Spike threshold adjustment parameter (mV) TauMax = 0.608 # Max. adaptation decay of slow non-inactivating Potassium current (s) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate', 'p': 'iM gate', 'q': 'iCaL activation gate', 'r': 'iCaL inactivation gate' } # ------------------------------ Gating states kinetics ------------------------------ @classmethod def alphaq(cls, Vm): return 0.055 * cls.vtrap(-(Vm + 27), 3.8) * 1e3 # s-1 @staticmethod def betaq(Vm): return 0.94 * np.exp(-(Vm + 75) / 17) * 1e3 # s-1 @staticmethod def alphar(Vm): return 0.000457 * np.exp(-(Vm + 13) / 50) * 1e3 # s-1 @staticmethod def betar(Vm): return 0.0065 / (np.exp(-(Vm + 15) / 28) + 1) * 1e3 # s-1 # ------------------------------ States derivatives ------------------------------ @classmethod def derStates(cls): return {**super().derStates(), **{ 'q': lambda Vm, x: cls.alphaq(Vm) * (1 - x['q']) - cls.betaq(Vm) * x['q'], 'r': lambda Vm, x: cls.alphar(Vm) * (1 - x['r']) - cls.betar(Vm) * x['r'] }} # ------------------------------ Steady states ------------------------------ @classmethod def steadyStates(cls): return {**super().steadyStates(), **{ 'q': lambda Vm: cls.alphaq(Vm) / (cls.alphaq(Vm) + cls.betaq(Vm)), 'r': lambda Vm: cls.alphar(Vm) / (cls.alphar(Vm) + cls.betar(Vm)) }} # ------------------------------ Membrane currents ------------------------------ @classmethod def iCaL(cls, q, r, Vm): ''' high-threshold (L-type) Calcium current ''' return cls.gCaLbar * q**2 * r * (Vm - cls.ECa) # mA/m2 @classmethod def currents(cls): return {**super().currents(), **{ 'iCaL': lambda Vm, x: cls.iCaL(x['q'], x['r'], Vm) }} diff --git a/PySONIC/neurons/template.py b/PySONIC/neurons/template.py index 15adbbe..20f3388 100644 --- a/PySONIC/neurons/template.py +++ b/PySONIC/neurons/template.py @@ -1,114 +1,114 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2019-06-11 15:58:38 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-08-22 15:19:33 +# @Last Modified time: 2019-09-12 11:15:17 import numpy as np from ..core import PointNeuron class TemplateNeuron(PointNeuron): ''' Template neuron class ''' # Neuron name name = 'template' # ------------------------------ Biophysical parameters ------------------------------ # Resting parameters Cm0 = 1e-2 # Membrane capacitance (F/m2) Vm0 = -71.9 # Membrane potential (mV) # Reversal potentials (mV) ENa = 50.0 # Sodium EK = -90.0 # Potassium ELeak = -70.3 # Non-specific leakage # Maximal channel conductances (S/m2) gNabar = 560.0 # Sodium gKdbar = 60.0 # Delayed-rectifier Potassium gLeak = 0.205 # Non-specific leakage # Additional parameters VT = -56.2 # Spike threshold adjustment parameter (mV) # ------------------------------ States names & descriptions ------------------------------ states = { 'm': 'iNa activation gate', 'h': 'iNa inactivation gate', 'n': 'iKd gate' } # ------------------------------ Gating states kinetics ------------------------------ @classmethod def alpham(cls, Vm): return 0.32 * cls.vtrap(13 - (Vm - cls.VT), 4) * 1e3 # s-1 @classmethod def betam(cls, Vm): return 0.28 * cls.vtrap((Vm - cls.VT) - 40, 5) * 1e3 # s-1 @classmethod def alphah(cls, Vm): return 0.128 * np.exp(-((Vm - cls.VT) - 17) / 18) * 1e3 # s-1 @classmethod def betah(cls, Vm): return 4 / (1 + np.exp(-((Vm - cls.VT) - 40) / 5)) * 1e3 # s-1 @classmethod def alphan(cls, Vm): return 0.032 * cls.vtrap(15 - (Vm - cls.VT), 5) * 1e3 # s-1 @classmethod def betan(cls, Vm): return 0.5 * np.exp(-((Vm - cls.VT) - 10) / 40) * 1e3 # s-1 # ------------------------------ States derivatives ------------------------------ @classmethod def derStates(cls): return { 'm': lambda Vm, x: cls.alpham(Vm) * (1 - x['m']) - cls.betam(Vm) * x['m'], 'h': lambda Vm, x: cls.alphah(Vm) * (1 - x['h']) - cls.betah(Vm) * x['h'], 'n': lambda Vm, x: cls.alphan(Vm) * (1 - x['n']) - cls.betan(Vm) * x['n'] } # ------------------------------ Steady states ------------------------------ @classmethod def steadyStates(cls): return { 'm': lambda Vm: cls.alpham(Vm) / (cls.alpham(Vm) + cls.betam(Vm)), 'h': lambda Vm: cls.alphah(Vm) / (cls.alphah(Vm) + cls.betah(Vm)), 'n': lambda Vm: cls.alphan(Vm) / (cls.alphan(Vm) + cls.betan(Vm)) } # ------------------------------ Membrane currents ------------------------------ @classmethod def iNa(cls, m, h, Vm): ''' Sodium current ''' return cls.gNabar * m**3 * h * (Vm - cls.ENa) # mA/m2 @classmethod def iKd(cls, n, Vm): ''' delayed-rectifier Potassium current ''' return cls.gKdbar * n**4 * (Vm - cls.EK) # mA/m2 @classmethod def iLeak(cls, Vm): ''' non-specific leakage current ''' return cls.gLeak * (Vm - cls.ELeak) # mA/m2 @classmethod def currents(cls): return { 'iNa': lambda Vm, x: cls.iNa(x['m'], x['h'], Vm), 'iKd': lambda Vm, x: cls.iKd(x['n'], Vm), 'iLeak': lambda Vm, _: cls.iLeak(Vm) } diff --git a/PySONIC/parsers.py b/PySONIC/parsers.py index 5569f45..2136d2e 100644 --- a/PySONIC/parsers.py +++ b/PySONIC/parsers.py @@ -1,666 +1,678 @@ # -*- 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-09-06 16:07:23 +# @Last Modified time: 2019-09-12 16:24:53 import os import logging import pprint import numpy as np import matplotlib.pyplot as plt from argparse import ArgumentParser from .utils import Intensity2Pressure, selectDirDialog, OpenFilesDialog, isIterable from .neurons import getPointNeuron, CorticalRS from .plt import GroupedTimeSeries, CompTimeSeries DEFAULT_OUTPUT_FOLDER = os.path.abspath(os.path.split(__file__)[0] + '../../../../dump') 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 + ['all'], 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 os.path.abspath(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: if args['outputdir'] is not None and args['outputdir'] == 'dump': return DEFAULT_OUTPUT_FOLDER 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 @staticmethod def parsePlot(args, output): render_args = {} if 'spikes' in args: render_args['spikes'] = args['spikes'] if args['compare']: if args['plot'] == ['all']: logger.error('Specific variables must be specified for comparative plots') return for key in ['cmap', 'cscale']: if key in args: render_args[key] = args[key] for pltvar in args['plot']: comp_plot = CompTimeSeries(output, pltvar) comp_plot.render(**render_args) else: scheme_plot = GroupedTimeSeries(output, pltscheme=args['pltscheme']) scheme_plot.render(**render_args) plt.show() class TestParser(Parser): def __init__(self, valid_subsets): super().__init__() self.addProfiling() self.addSubset(valid_subsets) 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) 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(dep_key='save') self.addSave() self.addCompare() self.addCmap() self.addCscale() 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': CorticalRS.Cm0 * 1e2, # uF/m2 'Qm0': CorticalRS.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, nargs='+', help='Resting membrane capacitance (uF/cm2)') def addQm0(self): self.add_argument( '--Qm0', type=float, nargs='+', 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 @staticmethod def parseSimInputs(args): return [args[k] for k in ['freq', 'amp', 'charge']] class NeuronSimParser(SimParser): def __init__(self): super().__init__() self.defaults.update({ 'neuron': 'RS', 'tstim': 100.0, # ms 'toffset': 50. # ms }) self.factors.update({ 'tstim': 1e-3, 'toffset': 1e-3 }) self.addNeuron() self.addTstim() self.addToffset() 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)') class VClampParser(NeuronSimParser): def __init__(self): super().__init__() self.defaults.update({ 'vhold': -70.0, # mV 'vstep': 0.0 # mV }) self.factors.update({ 'vhold': 1., 'vstep': 1. }) self.addVhold() self.addVstep() def addVhold(self): self.add_argument( '--vhold', nargs='+', type=float, help='Held membrane potential (mV)') def addVstep(self): self.add_argument( '--vstep', nargs='+', type=float, help='Step membrane potential (mV)') + self.add_argument( + '--vsteprange', type=str, nargs='+', + help='Step membrane potential range {} (mV)'.format(self.dist_str)) + self.to_parse['vstep'] = self.parseVstep + + def parseVstep(self, args): + vstep_range, vstep = [args.pop(k) for k in ['vsteprange', 'vstep']] + if vstep_range is not None: + vsteps = self.getDistFromList(vstep_range) # mV + else: + vsteps = np.array(vstep) # mV + return vsteps def parse(self, args=None): if args is None: args = super().parse() for key in ['vhold', 'vstep', 'tstim', 'toffset']: args[key] = self.parse2array(args, key, factor=self.factors[key]) return args @staticmethod def parseSimInputs(args): return [args[k] for k in ['vhold', 'vstep', 'tstim', 'toffset']] class PWSimParser(NeuronSimParser): ''' Generic parser interface to run PW patterned simulations from the command line. ''' def __init__(self): super().__init__() self.defaults.update({ 'PRF': 100.0, # Hz 'DC': 100.0 # % }) self.factors.update({ 'PRF': 1., 'DC': 1e-2 }) self.allowed.update({ 'DC': range(101) }) self.addPRF() self.addDC() self.addTitrate() self.addSpikes() 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 @staticmethod def parseSimInputs(args): return [args[k] for k in ['amp', 'tstim', 'toffset', 'PRF', 'DC']] 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 addVext(self): self.add_argument( '--Vext', nargs='+', type=float, help='Extracellular potential (mV)') 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 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 @staticmethod def parseSimInputs(args): return [args['freq']] + PWSimParser.parseSimInputs(args) + [args[k] for k in ['fs', 'method']] diff --git a/PySONIC/plt/timeseries.py b/PySONIC/plt/timeseries.py index c4be5aa..f1d7834 100644 --- a/PySONIC/plt/timeseries.py +++ b/PySONIC/plt/timeseries.py @@ -1,493 +1,494 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2018-09-25 16:18:45 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-09-06 11:20:34 +# @Last Modified time: 2019-09-12 16:49:11 import numpy as np import matplotlib.pyplot as plt from ..postpro import detectSpikes, convertPeaksProperties from ..utils import * from .pltutils import * class TimeSeriesPlot(GenericPlot): ''' Generic interface to build a plot displaying temporal profiles of model simulations. ''' @classmethod def setTimeLabel(cls, ax, tplt, fs): return super().setXLabel(ax, tplt, fs) @classmethod def setYLabel(cls, ax, yplt, fs, grouplabel=None): if grouplabel is not None: yplt['label'] = grouplabel return super().setYLabel(ax, yplt, fs) def checkInputs(self, *args, **kwargs): return NotImplementedError @staticmethod def getStimStates(df): try: stimstate = df['stimstate'] except KeyError: stimstate = df['states'] return stimstate.values @classmethod def getStimPulses(cls, t, states): ''' Determine the onset and offset times of pulses from a stimulation vector. :param t: time vector (s). :param states: a vector of stimulation state (ON/OFF) at each instant in time. :return: 3-tuple with number of patches, timing of STIM-ON an STIM-OFF instants. ''' # Compute states derivatives and identify bounds indexes of pulses dstates = np.diff(states) ipulse_on = np.where(dstates > 0.0)[0] + 1 ipulse_off = np.where(dstates < 0.0)[0] + 1 if ipulse_off.size < ipulse_on.size: ioff = t.size - 1 if ipulse_off.size == 0: ipulse_off = np.array([ioff]) else: ipulse_off = np.insert(ipulse_off, ipulse_off.size - 1, ioff) # Get time instants for pulses ON and OFF tpulse_on = t[ipulse_on] tpulse_off = t[ipulse_off] return tpulse_on, tpulse_off @staticmethod def addLegend(ax, handles, labels, fs, color=None, ls=None): lh = ax.legend(handles, labels, loc=1, fontsize=fs, frameon=False) if color is not None: for l in lh.get_lines(): l.set_color(color) if ls: for l in lh.get_lines(): l.set_linestyle(ls) @classmethod def materializeSpikes(cls, ax, data, tplt, yplt, color, mode, add_to_legend=False): ispikes, properties = detectSpikes(data) t = data['t'].values Qm = data['Qm'].values if ispikes is not None: yoffset = 5 ax.plot(t[ispikes] * tplt['factor'], Qm[ispikes] * yplt['factor'] + yoffset, 'v', color=color, label='spikes' if add_to_legend else None) if mode == 'details': ileft = properties['left_bases'] iright = properties['right_bases'] properties = convertPeaksProperties(t, properties) ax.plot(t[ileft] * tplt['factor'], Qm[ileft] * yplt['factor'] - 5, '<', color=color, label='left-bases' if add_to_legend else None) ax.plot(t[iright] * tplt['factor'], Qm[iright] * yplt['factor'] - 10, '>', color=color, label='right-bases' if add_to_legend else None) ax.vlines( x=t[ispikes] * tplt['factor'], ymin=(Qm[ispikes] - properties['prominences']) * yplt['factor'], ymax=Qm[ispikes] * yplt['factor'], color=color, linestyles='dashed', label='prominences' if add_to_legend else '') ax.hlines( y=properties['width_heights'] * yplt['factor'], xmin=properties['left_ips'] * tplt['factor'], xmax=properties['right_ips'] * tplt['factor'], color=color, linestyles='dotted', label='half-widths' if add_to_legend else '') return add_to_legend @staticmethod def prepareTime(t, tplt): if tplt['onset'] > 0.0: - t = np.insert(t, 0, -tplt['onset']) + tonset = t.min() - 0.05 * np.ptp(t) + t = np.insert(t, 0, tonset) return t * tplt['factor'] @staticmethod def addPatches(ax, tpatch_on, tpatch_off, tplt, color='#8A8A8A'): for i in range(tpatch_on.size): ax.axvspan(tpatch_on[i] * tplt['factor'], tpatch_off[i] * tplt['factor'], edgecolor='none', facecolor=color, alpha=0.2) @staticmethod def plotInset(inset_ax, t, y, tplt, yplt, line, color, lw): inset_window = np.logical_and(t > (inset['xlims'][0] / tplt['factor']), t < (inset['xlims'][1] / tplt['factor'])) inset_ax.plot(t[inset_window] * tplt['factor'], y[inset_window] * yplt['factor'], linewidth=lw, linestyle=line, color=color) return inset_ax @staticmethod def addInsetPatches(ax, inset_ax, inset, tpatch_on, tpatch_off, tplt, color): ybottom, ytop = ax.get_ylim() cond_on = np.logical_and(tpatch_on > (inset['xlims'][0] / tfactor), tpatch_on < (inset['xlims'][1] / tfactor)) cond_off = np.logical_and(tpatch_off > (inset['xlims'][0] / tfactor), tpatch_off < (inset['xlims'][1] / tfactor)) cond_glob = np.logical_and(tpatch_on < (inset['xlims'][0] / tfactor), tpatch_off > (inset['xlims'][1] / tfactor)) cond_onoff = np.logical_or(cond_on, cond_off) cond = np.logical_or(cond_onoff, cond_glob) npatches_inset = np.sum(cond) for i in range(npatches_inset): inset_ax.add_patch(Rectangle((tpatch_on[cond][i] * tfactor, ybottom), (tpatch_off[cond][i] - tpatch_on[cond][i]) * tfactor, ytop - ybottom, color=color, alpha=0.1)) class CompTimeSeries(ComparativePlot, TimeSeriesPlot): ''' Interface to build a comparative plot displaying profiles of a specific output variable across different model simulations. ''' def __init__(self, filepaths, varname): ''' Constructor. :param filepaths: list of full paths to output data files to be compared :param varname: name of variable to extract and compare ''' ComparativePlot.__init__(self, filepaths, varname) def checkPatches(self, patches): greypatch = False if patches == 'none': patches = [False] * len(self.filepaths) elif patches == 'all': patches = [True] * len(self.filepaths) elif patches == 'one': patches = [True] + [False] * (len(self.filepaths) - 1) greypatch = True elif isinstance(patches, list): if len(patches) != len(self.filepaths): raise ValueError( 'Invalid patches ({}): not matching number of compared files ({})'.format( len(patches), len(self.filepaths))) if not all(isinstance(p, bool) for p in patches): raise TypeError('Invalid patch sequence: all list items must be boolean typed') else: raise ValueError( 'Invalid patches: must be either "none", all", "one", or a boolean list') return patches, greypatch def checkInputs(self, lines, labels, colors, patches): self.checkLabels(labels) lines = self.checkLines(lines) colors = self.checkColors(colors) patches, greypatch = self.checkPatches(patches) return lines, labels, colors, patches, greypatch @staticmethod def createBackBone(figsize): fig, ax = plt.subplots(figsize=figsize) ax.set_zorder(0) return fig, ax @classmethod def postProcess(cls, ax, tplt, yplt, fs, meta, prettify): cls.removeSpines(ax) if 'bounds' in yplt: ax.set_ylim(*yplt['bounds']) cls.setTimeLabel(ax, tplt, fs) cls.setYLabel(ax, yplt, fs) if prettify: cls.prettify(ax, xticks=(0, meta['tstim'] * tplt['factor'])) cls.setTickLabelsFontSize(ax, fs) def render(self, figsize=(11, 4), fs=10, lw=2, labels=None, colors=None, lines=None, patches='one', inset=None, frequency=1, spikes='none', cmap=None, cscale='lin', trange=None, prettify=False): ''' Render plot. :param figsize: figure size (x, y) :param fs: labels fontsize :param lw: linewidth :param labels: list of labels to use in the legend :param colors: list of colors to use for each curve :param lines: list of linestyles :param patches: string indicating whether/how to mark stimulation periods with rectangular patches :param inset: string indicating whether/how to mark an inset zooming on a particular region of the graph :param frequency: frequency at which to plot samples :param spikes: string indicating how to show spikes ("none", "marks" or "details") :param cmap: color map to use for colobar-based comparison (if not None) :param cscale: color scale to use for colobar-based comparison :param trange: optional lower and upper bounds to time axis :return: figure handle ''' lines, labels, colors, patches, greypatch = self.checkInputs( lines, labels, colors, patches) fcodes = [] fig, ax = self.createBackBone(figsize) if inset is not None: inset_ax = self.addInset(fig, ax, inset) # Loop through data files handles, comp_values, full_labels = [], [], [] tmin, tmax = np.inf, -np.inf for j, filepath in enumerate(self.filepaths): # Load data try: data, meta = self.getData(filepath, frequency, trange) except ValueError as err: continue if 'tcomp' in meta: meta.pop('tcomp') full_labels.append(self.figtitle(meta)) # Extract model model = self.getModel(meta) fcodes.append(model.filecode(meta)) # Check consistency of sim types and check differing inputs comp_values = self.checkConsistency(meta, comp_values) # Extract time and stim pulses t = data['t'].values stimstate = self.getStimStates(data) tpatch_on, tpatch_off = self.getStimPulses(t, stimstate) tplt = self.getTimePltVar(model.tscale) t = self.prepareTime(t, tplt) # Extract y-variable pltvars = model.getPltVars() if self.varname not in pltvars: raise KeyError( 'Unknown plot variable: "{}". Possible plot variables are: {}'.format( self.varname, ', '.join(['"{}"'.format(p) for p in pltvars.keys()]))) yplt = pltvars[self.varname] y = extractPltVar(model, yplt, data, meta, t.size, self.varname) # Plot time series handles.append(ax.plot(t, y, linewidth=lw, linestyle=lines[j], color=colors[j])[0]) # Optional: add spikes if self.varname == 'Qm' and spikes != 'none': self.materializeSpikes(ax, data, tplt, yplt, colors[j], spikes) # Plot optional inset if inset is not None: inset_ax = self.plotInset(inset_ax, t, y, tplt, yplt, lines[j], colors[j], lw) # Add optional STIM-ON patches if patches[j]: ybottom, ytop = ax.get_ylim() color = '#8A8A8A' if greypatch else handles[j].get_color() self.addPatches(ax, tpatch_on, tpatch_off, tplt, color) if inset is not None: self.addInsetPatches(ax, inset_ax, inset, tpatch_on, tpatch_off, tplt, color) tmin, tmax = min(tmin, t.min()), max(tmax, t.max()) # Determine labels if self.comp_ref_key is not None: self.comp_info = model.inputs().get(self.comp_ref_key, None) comp_values, comp_labels = self.getCompLabels(comp_values) labels = self.chooseLabels(labels, comp_labels, full_labels) # - split labels by space merge_keys = ['neuron'] common_label = self.getCommonLabel(full_labels.copy(), sep=' ', merge_keys=merge_keys) ax.set_title(common_label, fontsize=fs) # Post-process figure self.postProcess(ax, tplt, yplt, fs, meta, prettify) ax.set_xlim(tmin, tmax) fig.tight_layout() if inset is not None: self.materializeInset(ax, inset_ax, inset) # Add labels or colorbar legend if cmap is not None: if not self.is_unique_comp: raise ValueError('Colormap mode unavailable for multiple differing parameters') if self.comp_info is None: raise ValueError('Colormap mode unavailable for qualitative comparisons') self.addCmap( fig, cmap, handles, comp_values, self.comp_info, fs, prettify, zscale=cscale) else: self.addLegend(ax, handles, labels, fs) # Add window title based on common pattern common_fcode = self.getCommonLabel(fcodes.copy(), sep='_') fig.canvas.set_window_title(common_fcode) return fig class GroupedTimeSeries(TimeSeriesPlot): ''' Interface to build a plot displaying profiles of several output variables arranged into specific schemes. ''' def __init__(self, filepaths, pltscheme=None): ''' Constructor. :param filepaths: list of full paths to output data files to be compared :param varname: name of variable to extract and compare ''' super().__init__(filepaths) self.pltscheme = pltscheme @staticmethod def createBackBone(pltscheme): naxes = len(pltscheme) if naxes == 1: fig, ax = plt.subplots(figsize=(11, 4)) axes = [ax] else: fig, axes = plt.subplots(naxes, 1, figsize=(11, min(3 * naxes, 9))) return fig, axes @classmethod def postProcess(cls, axes, tplt, fs, meta, prettify): for ax in axes: cls.removeSpines(ax) # if prettify: # cls.prettify(ax, xticks=(0, meta['tstim'] * tplt['factor']), yfmt=None) cls.setTickLabelsFontSize(ax, fs) for ax in axes[:-1]: ax.set_xticklabels([]) cls.setTimeLabel(axes[-1], tplt, fs) def render(self, fs=10, lw=2, labels=None, colors=None, lines=None, patches='one', save=False, outputdir=None, fig_ext='png', frequency=1, spikes='none', trange=None, prettify=False): ''' Render plot. :param fs: labels fontsize :param lw: linewidth :param labels: list of labels to use in the legend :param colors: list of colors to use for each curve :param lines: list of linestyles :param patches: boolean indicating whether to mark stimulation periods with rectangular patches :param save: boolean indicating whether or not to save the figure(s) :param outputdir: path to output directory in which to save figure(s) :param fig_ext: string indcating figure extension ("png", "pdf", ...) :param frequency: frequency at which to plot samples :param spikes: string indicating how to show spikes ("none", "marks" or "details") :param trange: optional lower and upper bounds to time axis :return: figure handle(s) ''' figs = [] for filepath in self.filepaths: # Load data and extract model try: data, meta = self.getData(filepath, frequency, trange) except ValueError as err: continue model = self.getModel(meta) # Extract time and stim pulses t = data['t'].values stimstate = self.getStimStates(data) tpatch_on, tpatch_off = self.getStimPulses(t, stimstate) tplt = self.getTimePltVar(model.tscale) t = self.prepareTime(t, tplt) # Check plot scheme if provided, otherwise generate it pltvars = model.getPltVars() if self.pltscheme is not None: for key in list(sum(list(self.pltscheme.values()), [])): if key not in pltvars: raise KeyError('Unknown plot variable: "{}"'.format(key)) pltscheme = self.pltscheme else: pltscheme = model.getPltScheme() # Create figure fig, axes = self.createBackBone(pltscheme) # Loop through each subgraph for ax, (grouplabel, keys) in zip(axes, pltscheme.items()): ax_legend_spikes = False # Extract variables to plot nvars = len(keys) ax_pltvars = [pltvars[k] for k in keys] if nvars == 1: ax_pltvars[0]['color'] = 'k' ax_pltvars[0]['ls'] = '-' # Set y-axis unit and bounds self.setYLabel(ax, ax_pltvars[0].copy(), fs, grouplabel=grouplabel) if 'bounds' in ax_pltvars[0]: ax_min = min([ap['bounds'][0] for ap in ax_pltvars]) ax_max = max([ap['bounds'][1] for ap in ax_pltvars]) ax.set_ylim(ax_min, ax_max) # Plot time series icolor = 0 for yplt, name in zip(ax_pltvars, pltscheme[grouplabel]): color = yplt.get('color', 'C{}'.format(icolor)) y = extractPltVar(model, yplt, data, meta, t.size, name) ax.plot(t, y, yplt.get('ls', '-'), c=color, lw=lw, label='$\\rm {}$'.format(yplt['label'])) if 'color' not in yplt: icolor += 1 # Optional: add spikes if name == 'Qm' and spikes != 'none': ax_legend_spikes = self.materializeSpikes( ax, data, tplt, yplt, color, spikes, add_to_legend=True) # Add legend if nvars > 1 or 'gate' in ax_pltvars[0]['desc'] or ax_legend_spikes: ax.legend(fontsize=fs, loc=7, ncol=nvars // 4 + 1, frameon=False) # Set x-limits and add optional patches for ax in axes: ax.set_xlim(t.min(), t.max()) if patches != 'none': self.addPatches(ax, tpatch_on, tpatch_off, tplt) # Post-process figure self.postProcess(axes, tplt, fs, meta, prettify) axes[0].set_title(self.figtitle(meta), fontsize=fs) fig.tight_layout() fig.canvas.set_window_title(model.filecode(meta)) # Save figure if needed (automatic or checked) if save: filecode = model.filecode(meta) if outputdir is None: outputdir = os.path.split(filepath)[0] plt_filename = '{}/{}.{}'.format(outputdir, filecode, fig_ext) plt.savefig(plt_filename) logger.info('Saving figure as "{}"'.format(plt_filename)) plt.close() figs.append(fig) return figs if __name__ == '__main__': # example of use filepaths = OpenFilesDialog('pkl')[0] comp_plot = CompTimeSeries(filepaths, 'Qm') fig = comp_plot.render( lines=['-', '--'], labels=['60 kPa', '80 kPa'], patches='one', colors=['r', 'g'], xticks=[0, 100], yticks=[-80, +50], inset={'xcoords': [5, 40], 'ycoords': [-35, 45], 'xlims': [57.5, 60.5], 'ylims': [10, 35]} ) scheme_plot = GroupedTimeSeries(filepaths) figs = scheme_plot.render() plt.show()