diff --git a/PySONIC/utils.py b/PySONIC/utils.py index ee9d70e..dadea23 100644 --- a/PySONIC/utils.py +++ b/PySONIC/utils.py @@ -1,304 +1,304 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-09-19 22:30:46 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-25 17:17:40 +# @Last Modified time: 2018-09-25 18:03:00 """ Definition of generic utility functions used in other modules """ import operator import os import pickle import tkinter as tk from tkinter import filedialog import numpy as np import colorlog from scipy.interpolate import interp1d def setLogger(): log_formatter = colorlog.ColoredFormatter( '%(log_color)s %(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:', reset=True, log_colors={ 'DEBUG': 'green', 'INFO': 'white', 'WARNING': 'yellow', 'ERROR': 'red', 'CRITICAL': 'red,bg_white', }, style='%' ) log_handler = colorlog.StreamHandler() log_handler.setFormatter(log_formatter) color_logger = colorlog.getLogger('PySONIC') color_logger.addHandler(log_handler) return color_logger # Get package logger logger = setLogger() # SI units prefixes si_prefixes = { 'y': 1e-24, # yocto 'z': 1e-21, # zepto 'a': 1e-18, # atto 'f': 1e-15, # femto 'p': 1e-12, # pico 'n': 1e-9, # nano 'u': 1e-6, # micro 'm': 1e-3, # mili '': 1e0, # None 'k': 1e3, # kilo 'M': 1e6, # mega 'G': 1e9, # giga 'T': 1e12, # tera 'P': 1e15, # peta 'E': 1e18, # exa 'Z': 1e21, # zetta 'Y': 1e24, # yotta } def si_format(x, precision=0, space=''): ''' Format a float according to the SI unit system, with the appropriate prefix letter. ''' if isinstance(x, float) or isinstance(x, int) or isinstance(x, np.float) or\ isinstance(x, np.int32) or isinstance(x, np.int64): if x == 0: factor = 1e0 prefix = '' else: sorted_si_prefixes = sorted(si_prefixes.items(), key=operator.itemgetter(1)) vals = [tmp[1] for tmp in sorted_si_prefixes] # vals = list(si_prefixes.values()) ix = np.searchsorted(vals, np.abs(x)) - 1 if np.abs(x) == vals[ix + 1]: ix += 1 factor = vals[ix] prefix = sorted_si_prefixes[ix][0] # prefix = list(si_prefixes.keys())[ix] return '{{:.{}f}}{}{}'.format(precision, space, prefix).format(x / factor) elif isinstance(x, list) or isinstance(x, tuple): return [si_format(item, precision, space) for item in x] elif isinstance(x, np.ndarray) and x.ndim == 1: return [si_format(float(item), precision, space) for item in x] else: print(type(x)) def pow10_format(number, precision=2): ''' Format a number in power of 10 notation. ''' ret_string = '{0:.{1:d}e}'.format(number, precision) a, b = ret_string.split("e") a = float(a) b = int(b) return '{}10^{{{}}}'.format('{} * '.format(a) if a != 1. else '', b) def rmse(x1, x2): """ Compute the root mean square error between two 1D arrays """ return np.sqrt(((x1 - x2) ** 2).mean()) def rsquared(x1, x2): ''' compute the R-squared coefficient between two 1D arrays ''' residuals = x1 - x2 ss_res = np.sum(residuals**2) ss_tot = np.sum((x1 - np.mean(x1))**2) return 1 - (ss_res / ss_tot) def Pressure2Intensity(p, rho=1075.0, c=1515.0): """ Return the spatial peak, pulse average acoustic intensity (ISPPA) associated with the specified pressure amplitude. :param p: pressure amplitude (Pa) :param rho: medium density (kg/m3) :param c: speed of sound in medium (m/s) :return: spatial peak, pulse average acoustic intensity (W/m2) """ return p**2 / (2 * rho * c) def Intensity2Pressure(I, rho=1075.0, c=1515.0): """ Return the pressure amplitude associated with the specified spatial peak, pulse average acoustic intensity (ISPPA). :param I: spatial peak, pulse average acoustic intensity (W/m2) :param rho: medium density (kg/m3) :param c: speed of sound in medium (m/s) :return: pressure amplitude (Pa) """ return np.sqrt(2 * rho * c * I) def OpenFilesDialog(filetype, dirname=''): """ Open a FileOpenDialogBox to select one or multiple file. The default directory and file type are given. :param dirname: default directory :param filetype: default file type :return: tuple of full paths to the chosen filenames """ root = tk.Tk() root.withdraw() filenames = filedialog.askopenfilenames(filetypes=[(filetype + " files", '.' + filetype)], initialdir=dirname) if filenames: par_dir = os.path.abspath(os.path.join(filenames[0], os.pardir)) else: par_dir = None return (filenames, par_dir) def selectDirDialog(): """ Open a dialog box to select a directory. :return: full path to selected directory """ root = tk.Tk() root.withdraw() return filedialog.askdirectory() def SaveFileDialog(filename, dirname=None, ext=None): ''' Open a dialog box to save file. :param filename: filename :param dirname: initial directory :param ext: default extension :return: full path to the chosen filename ''' root = tk.Tk() root.withdraw() filename_out = filedialog.asksaveasfilename( defaultextension=ext, initialdir=dirname, initialfile=filename) return filename_out def downsample(t_dense, y, nsparse): """ Decimate periodic signals to a specified number of samples.""" if(y.ndim) > 1: nsignals = y.shape[0] else: nsignals = 1 y = np.array([y]) # determine time step and period of input signal T = t_dense[-1] - t_dense[0] dt_dense = t_dense[1] - t_dense[0] # resample time vector linearly t_ds = np.linspace(t_dense[0], t_dense[-1], nsparse) # create MAV window nmav = int(0.03 * T / dt_dense) if nmav % 2 == 0: nmav += 1 mav = np.ones(nmav) / nmav # determine signals padding npad = int((nmav - 1) / 2) # determine indexes of sampling on convolved signals ids = np.round(np.linspace(0, t_dense.size - 1, nsparse)).astype(int) y_ds = np.empty((nsignals, nsparse)) # loop through signals for i in range(nsignals): # pad, convolve and resample pad_left = y[i, -(npad + 2):-2] pad_right = y[i, 1:npad + 1] y_ext = np.concatenate((pad_left, y[i, :], pad_right), axis=0) y_mav = np.convolve(y_ext, mav, mode='valid') y_ds[i, :] = y_mav[ids] if nsignals == 1: y_ds = y_ds[0, :] return (t_ds, y_ds) def rescale(x, lb=None, ub=None, lb_new=0, ub_new=1): ''' Rescale a value to a specific interval by linear transformation. ''' if lb is None: lb = x.min() if ub is None: ub = x.max() xnorm = (x - lb) / (ub - lb) return xnorm * (ub_new - lb_new) + lb_new def extractCompTimes(filenames): ''' Extract computation times from a list of simulation files. ''' tcomps = np.empty(len(filenames)) for i, fn in enumerate(filenames): logger.info('Loading data from "%s"', fn) with open(fn, 'rb') as fh: frame = pickle.load(fh) meta = frame['meta'] tcomps[i] = meta['tcomp'] return tcomps def getNeuronLookupsFile(mechname): return os.path.join( os.path.split(__file__)[0], 'neurons', '{}_lookups.pkl'.format(mechname)) def getLookups2D(mechname, a, Fdrive): ''' Retrieve appropriate 2D lookup tables and reference vectors for a given membrane mechanism, sonophore diameter and US frequency. :param mechname: name of membrane density mechanism :param a: sonophore diameter (m) :param Fdrive: US frequency (Hz) :return: 3-tuple with 1D numpy arrays of reference acoustic amplitudes and charge densities, and a dictionary of 2D lookup numpy arrays ''' # Check lookup file existence lookup_path = getNeuronLookupsFile(mechname) if not os.path.isfile(lookup_path): raise FileNotFoundError('Missing lookup file: "{}"'.format(lookup_path)) # Load lookups dictionary logger.debug('Loading lookup table') with open(lookup_path, 'rb') as fh: lookups4D = pickle.load(fh) # Retrieve 1D inputs from lookups dictionary aref = lookups4D.pop('a') Fref = lookups4D.pop('f') Aref = lookups4D.pop('A') Qref = lookups4D.pop('Q') # Check that sonophore diameter is within lookup range arange = (aref.min() - 1e-12, aref.max() + 1e-12) if a < arange[0] or a > arange[1]: raise ValueError('Invalid sonophore diameter: {}m (must be within {}m - {}m lookup interval)' .format(*si_format([a, *arange], precision=2, space=' '))) # Check that US frequency is within lookup range Frange = (Fref.min() - 1e-9, Fref.max() + 1e-9) if Fdrive < Frange[0] or Fdrive > Frange[1]: raise ValueError('Invalid frequency: {}Hz (must be within {}Hz - {}Hz lookup interval)' .format(*si_format([Fdrive, *Frange], precision=2, space=' '))) # Interpolate 4D lookups at sonophore diameter and then at US frequency logger.debug('Interpolating lookups at a = {}m'.format(si_format(a, space=' '))) lookups3D = {key: interp1d(aref, y4D, axis=0)(a) for key, y4D in lookups4D.items()} logger.debug('Interpolating lookups at f = {}Hz'.format(si_format(Fdrive, space=' '))) lookups2D = {key: interp1d(Fref, y3D, axis=0)(Fdrive) for key, y3D in lookups3D.items()} return Aref, Qref, lookups2D diff --git a/scripts/run_astim.py b/scripts/run_astim.py index d82fb6b..eaa7d72 100644 --- a/scripts/run_astim.py +++ b/scripts/run_astim.py @@ -1,148 +1,149 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-02-13 18:16:09 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-25 14:34:49 +# @Last Modified time: 2018-09-25 17:59:55 ''' Run A-STIM simulations of a specific point-neuron. ''' import os import logging import numpy as np +import matplotlib.pyplot as plt from argparse import ArgumentParser from PySONIC.core import NeuronalBilayerSonophore from PySONIC.utils import logger, selectDirDialog, Intensity2Pressure from PySONIC.neurons import getNeuronsDict from PySONIC.batches import createSimQueue, runBatch from PySONIC.plt import plotBatch # Default parameters defaults = dict( neuron='RS', diams=[32.0], # nm freqs=[500.0], # kHz amps=[100.0], # kPa durations=[100.0], # ms PRFs=[100.0], # Hz DCs=[100.0], # % offsets=[50.], # ms method='sonic' ) def runAStimBatch(outdir, nbls, stim_params, method, mpi=False): ''' Run batch A-STIM simulations of the system for various neuron types and stimulation parameters. :param outdir: full path to output directory :param stim_params: dictionary containing sweeps for all stimulation parameters :param method: numerical integration method ("classic", "hybrid" or "sonic") :param mpi: boolean statting wether or not to use multiprocessing :return: list of full paths to the output files ''' mandatory_params = ['freqs', 'durations', 'offsets', 'PRFs', 'DCs'] for mparam in mandatory_params: if mparam not in stim_params: raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam)) logger.info("Starting A-STIM simulation batch") # Generate queue nofreq_queue = createSimQueue(stim_params.get('amps', [None]), stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DCs']) # Repeat queue for each US frequency nofreq_queue = np.array(nofreq_queue) freqs = stim_params['freqs'] nf = len(freqs) nqueue = nofreq_queue.shape[0] queue = np.tile(nofreq_queue, (nf, 1)) freqs_col = np.vstack([np.ones(nqueue) * f for f in freqs]).reshape(nf * nqueue, 1) queue = np.hstack((freqs_col, queue)).tolist() # Add method to queue items for item in queue: item.append(method) # Run batch return runBatch(nbls, 'runAndSave', queue, extra_params=[outdir], mpi=mpi) def main(): ap = ArgumentParser() # Runtime options ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results') ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory') ap.add_argument('-t', '--titrate', default=False, action='store_true', help='Perform titration') ap.add_argument('-m', '--method', type=str, default=defaults['method'], help='Numerical integration method ("classic", "hybrid" or "sonic")') # Stimulation parameters ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'], help='Neuron name (string)') ap.add_argument('-a', '--diams', nargs='+', type=float, help='Sonophore diameter (nm)') ap.add_argument('-f', '--freqs', nargs='+', type=float, help='US frequency (kHz)') ap.add_argument('-A', '--amps', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)') ap.add_argument('-I', '--intensities', nargs='+', type=float, help='Acoustic intensity (W/cm2)') ap.add_argument('-d', '--durations', nargs='+', type=float, help='Stimulus duration (ms)') ap.add_argument('--offset', nargs='+', type=float, help='Offset duration (ms)') ap.add_argument('--PRF', nargs='+', type=float, help='PRF (Hz)') ap.add_argument('--DC', nargs='+', type=float, help='Duty cycle (%%)') # Parse arguments args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None} loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO logger.setLevel(loglevel) outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog() mpi = args['mpi'] plot = args['plot'] titrate = args['titrate'] method = args['method'] neuron_str = args['neuron'] diams = np.array(args.get('diams', defaults['diams'])) * 1e-9 # m if 'amps' in args: amps = np.array(args['amps']) * 1e3 # Pa elif 'intensities' in args: amps = Intensity2Pressure(np.array(args['intensities']) * 1e4) # Pa else: amps = np.array(defaults['amps']) * 1e3 # Pa stim_params = dict( freqs=np.array(args.get('freqs', defaults['freqs'])) * 1e3, # Hz amps=amps, # Pa durations=np.array(args.get('durations', defaults['durations'])) * 1e-3, # s PRFs=np.array(args.get('PRFs', defaults['PRFs'])), # Hz DCs=np.array(args.get('DCs', defaults['DCs'])) * 1e-2, # (-) offsets=np.array(args.get('offsets', defaults['offsets'])) * 1e-3 # s ) if titrate: stim_params['amps'] = [None] # Run A-STIM batch if neuron_str not in getNeuronsDict(): logger.error('Unknown neuron type: "%s"', neuron_str) return neuron = getNeuronsDict()[neuron_str]() pkl_filepaths = [] for a in diams: nbls = NeuronalBilayerSonophore(a, neuron) pkl_filepaths += runAStimBatch(outdir, nbls, stim_params, method, mpi=mpi) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles if plot: - yvars = {'Q_m': ['Qm']} - plotBatch(pkl_dir, pkl_filepaths, yvars) + plotBatch(pkl_filepaths, vars_dict={'Q_m': ['Qm']}) + plt.show() if __name__ == '__main__': main() diff --git a/scripts/run_estim.py b/scripts/run_estim.py index 818514b..333119a 100644 --- a/scripts/run_estim.py +++ b/scripts/run_estim.py @@ -1,108 +1,109 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-08-24 11:55:07 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-25 10:19:07 +# @Last Modified time: 2018-09-25 17:58:47 ''' Run E-STIM simulations of a specific point-neuron. ''' import os import logging import numpy as np +import matplotlib.pyplot as plt from argparse import ArgumentParser from PySONIC.utils import logger, selectDirDialog from PySONIC.neurons import * from PySONIC.batches import createSimQueue, runBatch from PySONIC.plt import plotBatch # Default parameters defaults = dict( neuron='RS', amps=[10.0], # mA/m2 durations=[100.0], # ms PRFs=[100.0], # Hz DCs=[100.0], # % offsets=[50.], # ms method='sonic' ) def runEStimBatch(outdir, neuron, stim_params, mpi=False): ''' Run batch E-STIM simulations of the system for various neuron types and stimulation parameters. :param outdir: full path to output directory :param stim_params: dictionary containing sweeps for all stimulation parameters :param mpi: boolean statting wether or not to use multiprocessing :return: list of full paths to the output files ''' mandatory_params = ['durations', 'offsets', 'PRFs', 'DCs'] for mparam in mandatory_params: if mparam not in stim_params: raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam)) logger.info("Starting E-STIM simulation batch") # Generate simulations queue queue = createSimQueue(stim_params.get('amps', [None]), stim_params['durations'], stim_params['offsets'], stim_params['PRFs'], stim_params['DCs']) # Run batch return runBatch(neuron, 'runAndSave', queue, extra_params=[outdir], mpi=mpi) def main(): ap = ArgumentParser() # Runtime options ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results') ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory') ap.add_argument('-t', '--titrate', default=False, action='store_true', help='Perform titration') # Stimulation parameters ap.add_argument('-n', '--neuron', type=str, default=defaults['neuron'], help='Neuron name (string)') ap.add_argument('-A', '--amps', nargs='+', type=float, help='Injected current density (mA/m2)') ap.add_argument('-d', '--durations', nargs='+', type=float, help='Stimulus duration (ms)') ap.add_argument('--offset', nargs='+', type=float, help='Offset duration (ms)') ap.add_argument('--PRF', nargs='+', type=float, help='PRF (Hz)') ap.add_argument('--DC', nargs='+', type=float, help='Duty cycle (%%)') # Parse arguments args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None} loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO logger.setLevel(loglevel) outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog() mpi = args['mpi'] plot = args['plot'] titrate = args['titrate'] neuron_str = args['neuron'] stim_params = dict( amps=np.array(args.get('amps', defaults['amps'])), # mA/m2 durations=np.array(args.get('durations', defaults['durations'])) * 1e-3, # s PRFs=np.array(args.get('PRFs', defaults['PRFs'])), # Hz DCs=np.array(args.get('DCs', defaults['DCs'])) * 1e-2, # (-) offsets=np.array(args.get('offsets', defaults['offsets'])) * 1e-3 # s ) if titrate: stim_params['amps'] = [None] # Run E-STIM batch if neuron_str not in getNeuronsDict(): logger.error('Unknown neuron type: "%s"', neuron_str) return neuron = getNeuronsDict()[neuron_str]() pkl_filepaths = runEStimBatch(outdir, neuron, stim_params, mpi=mpi) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles if plot: - yvars = {'V_m': ['Vm']} - plotBatch(pkl_dir, pkl_filepaths, yvars) + plotBatch(pkl_filepaths, vars_dict={'V_m': ['Vm']}) + plt.show() if __name__ == '__main__': main() diff --git a/scripts/run_mech.py b/scripts/run_mech.py index 9b3f7ef..2e89810 100644 --- a/scripts/run_mech.py +++ b/scripts/run_mech.py @@ -1,120 +1,122 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-11-21 10:46:56 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-25 14:36:20 +# @Last Modified time: 2018-09-25 18:01:00 ''' Run simulations of the NICE mechanical model. ''' import os import logging import numpy as np +import matplotlib.pyplot as plt from argparse import ArgumentParser from PySONIC.core import BilayerSonophore from PySONIC.utils import logger, selectDirDialog, Intensity2Pressure from PySONIC.neurons import CorticalRS from PySONIC.batches import createQueue, runBatch from PySONIC.plt import plotBatch # Default parameters defaults = dict( Cm0=CorticalRS().Cm0 * 1e2, # uF/m2 Qm0=CorticalRS().Vm0, # nC/m2 diams=[32.0], # nm embeddings=[0.], # um freqs=[500.0], # kHz amps=[100.0], # kPa charges=[0.] # nC/cm2 ) def runMechBatch(outdir, bls, stim_params, mpi=False): ''' Run batch simulations of the mechanical system with imposed values of charge density. :param outdir: full path to output directory :param bls: BilayerSonophore object :param stim_params: dictionary containing sweeps for all stimulation parameters :param mpi: boolean stating whether or not to use multiprocessing :return: list of full paths to the output files ''' # Checking validity of stimulation parameters mandatory_params = ['freqs', 'amps', 'charges'] for mparam in mandatory_params: if mparam not in stim_params: raise ValueError('Missing stimulation parameter field: "{}"'.format(mparam)) logger.info("Starting mechanical simulation batch") # Unpack stimulation parameters freqs = np.array(stim_params['freqs']) amps = np.array(stim_params['amps']) charges = np.array(stim_params['charges']) # Generate simulations queue and run batch queue = createQueue((freqs, amps, charges)) return runBatch(bls, 'runAndSave', queue, extra_params=[outdir], mpi=mpi) def main(): ap = ArgumentParser() # Runtime options ap.add_argument('--mpi', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-p', '--plot', default=False, action='store_true', help='Plot results') ap.add_argument('-o', '--outputdir', type=str, default=None, help='Output directory') # Stimulation parameters ap.add_argument('-a', '--diams', nargs='+', type=float, help='Sonophore diameter (nm)') ap.add_argument('--Cm0', type=float, default=defaults['Cm0'], help='Resting membrane capacitance (uF/cm2)') ap.add_argument('--Qm0', type=float, default=defaults['Qm0'], help='Resting membrane charge density (nC/cm2)') ap.add_argument('-d', '--embeddings', nargs='+', type=float, help='Embedding depth (um)') ap.add_argument('-f', '--freqs', nargs='+', type=float, help='US frequency (kHz)') ap.add_argument('-A', '--amps', nargs='+', type=float, help='Acoustic pressure amplitude (kPa)') ap.add_argument('-I', '--intensities', nargs='+', type=float, help='Acoustic intensity (W/cm2)') ap.add_argument('-Q', '--charges', nargs='+', type=float, help='Membrane charge density (nC/cm2)') # Parse arguments args = {key: value for key, value in vars(ap.parse_args()).items() if value is not None} loglevel = logging.DEBUG if args['verbose'] is True else logging.INFO logger.setLevel(loglevel) outdir = args['outputdir'] if 'outputdir' in args else selectDirDialog() mpi = args['mpi'] plot = args['plot'] Cm0 = args['Cm0'] * 1e-2 # F/m2 Qm0 = args['Qm0'] * 1e-5 # C/m2 diams = np.array(args.get('diams', defaults['diams'])) * 1e-9 # m embeddings = np.array(args.get('embeddings', defaults['embeddings'])) * 1e-6 # m if 'amps' in args: amps = np.array(args['amps']) * 1e3 # Pa elif 'intensities' in args: amps = Intensity2Pressure(np.array(args['intensities']) * 1e4) # Pa else: amps = np.array(defaults['amps']) * 1e3 # Pa stim_params = dict( freqs=np.array(args.get('freqs', defaults['freqs'])) * 1e3, # Hz amps=amps, # Pa charges=np.array(args.get('charges', defaults['charges'])) * 1e-5 # C/m2 ) # Run MECH batch pkl_filepaths = [] for a in diams: for d in embeddings: bls = BilayerSonophore(a, Cm0, Qm0, embedding_depth=d) pkl_filepaths += runMechBatch(outdir, bls, stim_params, mpi=mpi) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles if plot: - plotBatch(pkl_dir, pkl_filepaths) + plotBatch(pkl_filepaths, vars_dict={'Z': ['Z']}) + plt.show() if __name__ == '__main__': main()