diff --git a/PySONIC/lookups/BLS_lookups.json b/PySONIC/lookups/BLS_lookups.json index 7275c93..31ab9bb 100644 --- a/PySONIC/lookups/BLS_lookups.json +++ b/PySONIC/lookups/BLS_lookups.json @@ -1,436 +1,447 @@ { "32.0": { "-80.00": { "LJ_approx": { "x0": 1.7875580514692446e-09, "C": 14506.791031634148, "nrep": 3.911252335063797, "nattr": 0.9495868868453603 }, "Delta_eq": 1.2344323203763867e-09 }, "-71.40": { "LJ_approx": { "x0": 1.710159362626304e-09, "C": 16757.44053535206, "nrep": 3.9149844779001572, "nattr": 0.9876139143736086 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.588820457014353e-09, "C": 21124.28839722447, "nrep": 3.9219530533405726, "nattr": 1.0531179666960837 }, "Delta_eq": 1.302942739961778e-09 }, "-71.90": { "LJ_approx": { "x0": 1.7142977983903395e-09, "C": 16627.43538695451, "nrep": 3.9147721975981384, "nattr": 0.9855168537576823 }, "Delta_eq": 1.2553492695740507e-09 }, "-89.50": { "LJ_approx": { "x0": 1.8913883171160228e-09, "C": 12016.525797229067, "nrep": 3.9069373029335464, "nattr": 0.9021994595277029 }, "Delta_eq": 1.2107230911508513e-09 }, "-61.93": { "LJ_approx": { "x0": 1.6390264131559902e-09, "C": 19180.97634755811, "nrep": 3.9188840789597705, "nattr": 1.0250251620607604 }, "Delta_eq": 1.281743450351987e-09 }, "-53.00": { "LJ_approx": { "x0": 1.5830321174402216e-09, "C": 21361.655211483354, "nrep": 3.9223254792281588, "nattr": 1.0564645995714745 }, "Delta_eq": 1.305609024046854e-09 }, "-53.58": { "LJ_approx": { "x0": 1.5863754403940291e-09, "C": 21224.206759769622, "nrep": 3.922109852077156, "nattr": 1.0545313303221624 }, "Delta_eq": 1.3040630712174578e-09 }, "-48.87": { "LJ_approx": { "x0": 1.5603070731170595e-09, "C": 22321.280954434333, "nrep": 3.9238276518118833, "nattr": 1.0698008224359472 }, "Delta_eq": 1.3165739825437056e-09 }, "0.00": { "LJ_approx": { "x0": 1.429523524073023e-09, "C": 28748.036227122713, "nrep": 3.9338919786768276, "nattr": 1.1551044201542804 }, "Delta_eq": 1.4e-09 }, "-70.00": { "LJ_approx": { "x0": 1.698788510560293e-09, "C": 17120.631318195712, "nrep": 3.915575488491436, "nattr": 0.9934238714780391 }, "Delta_eq": 1.2603339470322538e-09 } }, "64.0": { "-80.00": { "LJ_approx": { "x0": 1.783357531675752e-09, "C": 14639.319598806138, "nrep": 3.9113027551187414, "nattr": 0.9404151935643594 }, "Delta_eq": 1.2344323203763867e-09 }, "-71.90": { "LJ_approx": { "x0": 1.7103254451796522e-09, "C": 16775.90747591089, "nrep": 3.9148582072320104, "nattr": 0.9747613204242506 }, "Delta_eq": 1.2553492695740507e-09 }, "-71.40": { "LJ_approx": { "x0": 1.7061996856525807e-09, "C": 16906.878806702443, "nrep": 3.9150725853841957, "nattr": 0.976778398349503 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.585264156646392e-09, "C": 21303.176047683613, "nrep": 3.92211004079812, "nattr": 1.0402641595550777 }, "Delta_eq": 1.302942739961778e-09 }, "-89.50": { "LJ_approx": { "x0": 1.886870747500225e-09, "C": 12129.260307725155, "nrep": 3.906945602966425, "nattr": 0.895598309376088 }, "Delta_eq": 1.2107230911508513e-09 }, "-61.93": { "LJ_approx": { "x0": 1.635297932833928e-09, "C": 19347.359224637305, "nrep": 3.9190113341505843, "nattr": 1.0129039638328117 }, "Delta_eq": 1.281743450351987e-09 } }, "50.0": { "-71.90": { "LJ_approx": { "x0": 1.7114794411958874e-09, "C": 16732.841575829825, "nrep": 3.914827883392826, "nattr": 0.9781401389550711 }, "Delta_eq": 1.2553492695740507e-09 } }, "10.0": { "0.00": { "LJ_approx": { "x0": 1.4403460578039628e-09, "C": 27932.27792195569, "nrep": 3.9334138654752686, "nattr": 1.19526523864855 }, "Delta_eq": 1.4e-09 } }, "100.0": { "0.00": { "LJ_approx": { "x0": 1.4254455131143225e-09, "C": 29048.417918044444, "nrep": 3.9342659189249254, "nattr": 1.1351227816904121 }, "Delta_eq": 1.4e-09 }, "-71.90": { "LJ_approx": { "x0": 1.7087681652667724e-09, "C": 16833.83962398515, "nrep": 3.914908533680663, "nattr": 0.9697102045586926 }, "Delta_eq": 1.2553492695740507e-09 }, "-71.40": { "LJ_approx": { "x0": 1.7046480280451768e-09, "C": 16965.15489682674, "nrep": 3.9151238997284845, "nattr": 0.9716928395857687 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.5838767580748841e-09, "C": 21372.412565003375, "nrep": 3.922191259312523, "nattr": 1.0344167918733638 }, "Delta_eq": 1.302942739961778e-09 }, "-89.50": { "LJ_approx": { "x0": 1.8850831984948754e-09, "C": 12173.857567383837, "nrep": 3.906959887367545, "nattr": 0.8923154523792497 }, "Delta_eq": 1.2107230911508513e-09 }, "-61.93": { "LJ_approx": { "x0": 1.6338405482612552e-09, "C": 19411.96069791924, "nrep": 3.9190798895919405, "nattr": 1.0073171165886772 }, "Delta_eq": 1.281743450351987e-09 } }, "500.0": { "0.00": { "LJ_approx": { "x0": 1.4236928207491738e-09, "C": 29174.423851140436, "nrep": 3.934489087598531, "nattr": 1.1244706663694597 }, "Delta_eq": 1.4e-09 } }, "15.0": { "-71.90": { "LJ_approx": { "x0": 1.722207527516432e-09, "C": 16331.124295102756, "nrep": 3.914721476976037, "nattr": 1.0021304453280393 }, "Delta_eq": 1.2553492695740507e-09 }, "-71.40": { "LJ_approx": { "x0": 1.7180439454408102e-09, "C": 16459.15520614834, "nrep": 3.9149304238974905, "nattr": 1.0043742068193204 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.5959361190370466e-09, "C": 20763.823187183425, "nrep": 3.921785737782814, "nattr": 1.0737976563473925 }, "Delta_eq": 1.302942739961778e-09 }, "-89.50": { "LJ_approx": { "x0": 1.9002701433896942e-09, "C": 11795.576706463997, "nrep": 3.9070059150621814, "nattr": 0.9114123498082551 }, "Delta_eq": 1.2107230911508513e-09 }, "-61.93": { "LJ_approx": { "x0": 1.646473044478369e-09, "C": 18846.9803104403, "nrep": 3.918766624746869, "nattr": 1.044220870098494 }, "Delta_eq": 1.281743450351987e-09 } }, "70.0": { "-61.93": { "LJ_approx": { "x0": 1.6349587829329923e-09, "C": 19362.426097728276, "nrep": 3.9190260877993097, "nattr": 1.0116609492531905 }, "Delta_eq": 1.281743450351987e-09 }, "-71.90": { "LJ_approx": { "x0": 1.7099628637265945e-09, "C": 16789.420291577666, "nrep": 3.9148688185890483, "nattr": 0.9736446481917082 }, "Delta_eq": 1.2553492695740507e-09 }, "-71.40": { "LJ_approx": { "x0": 1.7058388214243274e-09, "C": 16920.455606556396, "nrep": 3.9150834388621805, "nattr": 0.9756530742858303 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.584940743805642e-09, "C": 21319.35643207058, "nrep": 3.9221277053335264, "nattr": 1.0389567310289958 }, "Delta_eq": 1.302942739961778e-09 }, "-89.50": { "LJ_approx": { "x0": 1.886457143504779e-09, "C": 12139.584658299475, "nrep": 3.9069481854501382, "nattr": 0.8948872453823127 }, "Delta_eq": 1.2107230911508513e-09 } }, "150.0": { "-71.90": { "LJ_approx": { "x0": 1.707781542586275e-09, "C": 16870.340248462282, "nrep": 3.914948339378328, "nattr": 0.9660711186661354 }, "Delta_eq": 1.2553492695740507e-09 }, "-71.40": { "LJ_approx": { "x0": 1.7036651779798684e-09, "C": 17001.86272239105, "nrep": 3.9151643485118117, "nattr": 0.9680342060844059 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.5830041112065094e-09, "C": 21415.64649760579, "nrep": 3.9222512220871013, "nattr": 1.0303387653647968 }, "Delta_eq": 1.302942739961778e-09 }, "-89.50": { "LJ_approx": { "x0": 1.883936689519767e-09, "C": 12202.390459945682, "nrep": 3.906974649816042, "nattr": 0.8898102498996743 }, "Delta_eq": 1.2107230911508513e-09 }, "-61.93": { "LJ_approx": { "x0": 1.6329212539031057e-09, "C": 19452.44141782297, "nrep": 3.9191317215599946, "nattr": 1.0033715654432673 }, "Delta_eq": 1.281743450351987e-09 } }, "16.0": { "-71.90": { "LJ_approx": { "x0": 1.721337196662225e-09, "C": 16363.738563915058, "nrep": 3.9147202446411637, "nattr": 1.0005179992350384 }, "Delta_eq": 1.2553492695740507e-09 }, "-71.40": { "LJ_approx": { "x0": 1.717176984027311e-09, "C": 16491.97282711717, "nrep": 3.9149293522242252, "nattr": 1.0027563589061772 }, "Delta_eq": 1.2566584760815426e-09 }, "-54.00": { "LJ_approx": { "x0": 1.5951507107307484e-09, "C": 20803.713978702526, "nrep": 3.921796023871708, "nattr": 1.0717457519944718 }, "Delta_eq": 1.302942739961778e-09 }, "-89.50": { "LJ_approx": { "x0": 1.899300769652515e-09, "C": 11819.650350288994, "nrep": 3.906993085458305, "nattr": 0.9105930969008011 }, "Delta_eq": 1.2107230911508513e-09 }, "-61.93": { "LJ_approx": { "x0": 1.6456524054221337e-09, "C": 18883.851022856303, "nrep": 3.918771854603072, "nattr": 1.042336408142498 }, "Delta_eq": 1.281743450351987e-09 } + }, + "40.0": { + "-71.90": { + "LJ_approx": { + "x0": 1.7127556130633909e-09, + "C": 16685.139617894773, + "nrep": 3.9147997833590855, + "nattr": 0.9816110605284186 + }, + "Delta_eq": 1.2553492695740507e-09 + } } } \ No newline at end of file diff --git a/PySONIC/lookups/SONIC_RS.pkl b/PySONIC/lookups/SONIC_RS.pkl new file mode 100644 index 0000000..e869332 Binary files /dev/null and b/PySONIC/lookups/SONIC_RS.pkl differ diff --git a/PySONIC/utils.py b/PySONIC/utils.py index 8abcec0..78f418b 100644 --- a/PySONIC/utils.py +++ b/PySONIC/utils.py @@ -1,563 +1,561 @@ #!/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-19 15:43:48 +# @Last Modified time: 2018-09-19 20:07:27 """ Definition of generic utility functions used in other modules """ from enum import Enum import operator import os import pickle import tkinter as tk from tkinter import filedialog import inspect import json import yaml from openpyxl import load_workbook import numpy as np import colorlog from scipy.interpolate import interp1d from . import neurons 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() class InputError(Exception): ''' Exception raised for errors in the input. ''' pass class PmCompMethod(Enum): """ Enum: types of computation method for the intermolecular pressure """ direct = 1 predict = 2 def loadYAML(filepath): """ Load a dictionary of parameters from an external YAML file. :param filepath: full path to the YAML file :return: parameters dictionary """ _, filename = os.path.split(filepath) logger.debug('Loading parameters from "%s"', filename) with open(filepath, 'r') as f: stream = f.read() params = yaml.load(stream) return ParseNestedDict(params) def getLookupDir(): """ Return the location of the directory holding lookups files. :return: absolute path to the directory """ this_dir, _ = os.path.split(__file__) return os.path.join(this_dir, 'lookups') def ParseNestedDict(dict_in): """ Loop through a nested dictionary object and convert all string fields to floats. """ for key, value in dict_in.items(): if isinstance(value, dict): # If value itself is dictionary dict_in[key] = ParseNestedDict(value) elif isinstance(dict_in[key], str): dict_in[key] = float(dict_in[key]) return dict_in 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 ImportExcelCol(filename, sheetname, colstr, startrow): """ Load a specific column of an excel workbook as a numpy array. :param filename: absolute or relative path to the Excel workbook :param sheetname: name of the Excel spreadsheet to which data is appended :param colstr: string of the column to import :param startrow: index of the first row to consider :return: 1D numpy array with the column data """ wb = load_workbook(filename, read_only=True) ws = wb[sheetname] range_start_str = colstr + str(startrow) range_stop_str = colstr + str(ws.max_row) tmp = np.array([[i.value for i in j] for j in ws[range_start_str:range_stop_str]]) return tmp[:, 0] def ConstructMatrix(serialized_inputA, serialized_inputB, serialized_output): """ Construct a 2D output matrix from serialized input. :param serialized_inputA: serialized input variable A :param serialized_inputB: serialized input variable B :param serialized_output: serialized output variable :return: 4-tuple with vectors of unique values of A (m) and B (n), output variable 2D matrix (m,n) and number of holes in the matrix """ As = np.unique(serialized_inputA) Bs = np.unique(serialized_inputB) nA = As.size nB = Bs.size output = np.zeros((nA, nB)) output[:] = np.NAN nholes = 0 for i in range(nA): iA = np.where(serialized_inputA == As[i]) for j in range(nB): iB = np.where(serialized_inputB == Bs[j]) iMatch = np.intersect1d(iA, iB) if iMatch.size == 0: nholes += 1 elif iMatch.size > 1: logger.warning('Identical serialized inputs with values (%f, %f)', As[i], Bs[j]) else: iMatch = iMatch[0] output[i, j] = serialized_output[iMatch] return (As, Bs, output, nholes) 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 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 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 find_nearest(array, value): ''' Find nearest element in 1D array. ''' idx = (np.abs(array - value)).argmin() return (idx, array[idx]) 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 LennardJones(x, beta, alpha, C, m, n): """ Generic expression of a Lennard-Jones function, adapted for the context of symmetric deflection (distance = 2x). :param x: deflection (i.e. half-distance) :param beta: x-shifting factor :param alpha: x-scaling factor :param C: y-scaling factor :param m: exponent of the repulsion term :param n: exponent of the attraction term :return: Lennard-Jones potential at given distance (2x) """ return C * (np.power((alpha / (2 * x + beta)), m) - np.power((alpha / (2 * x + beta)), n)) def getNeuronsDict(): ''' Return dictionary of neurons classes that can be instantiated. ''' neurons_dict = {} for _, obj in inspect.getmembers(neurons): if inspect.isclass(obj) and isinstance(obj.name, str): neurons_dict[obj.name] = obj return neurons_dict def get_BLS_lookups(): lookup_path = getLookupDir() + '/BLS_lookups.json' try: with open(lookup_path) as fh: sample = json.load(fh) return sample except FileNotFoundError: return {} def save_BLS_lookups(lookups): """ Save BLS parameter into specific lookup file :return: absolute path to the directory """ lookup_path = getLookupDir() + '/BLS_lookups.json' with open(lookup_path, 'w') as fh: json.dump(lookups, fh, indent=2) 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 computeMeshEdges(x, scale='lin'): ''' Compute the appropriate edges of a mesh that quads a linear or logarihtmic distribution. :param x: the input vector :param scale: the type of distribution ('lin' for linear, 'log' for logarihtmic) :return: the edges vector ''' if scale is 'log': x = np.log10(x) dx = x[1] - x[0] if scale is 'lin': y = np.linspace(x[0] - dx / 2, x[-1] + dx / 2, x.size + 1) elif scale is 'log': y = np.logspace(x[0] - dx / 2, x[-1] + dx / 2, x.size + 1) return y 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 getCycleAverage(t, y, T): ''' Compute the cycle-averaged profile of a signal given a specific periodicity. :param t: time vector (s) :param y: signal vector :param T: period (s) :return: cycle-averaged signal vector ''' nsamples = y.size ncycles = int((t[-1] - t[0]) // T) npercycle = int(nsamples // ncycles) return np.mean(np.reshape(y[:ncycles * npercycle], (ncycles, npercycle)), axis=1) def itrpLookupsFreq(lookups3D, freqs, Fdrive): """ Interpolate a dictionary of 3D lookups at a given frequency. :param lookups3D: dictionary of 3D lookups :param freqs: array of lookup frequencies (Hz) :param Fdrive: acoustic drive frequency (Hz) :return: a dictionary of 2D lookups interpolated a the given frequency """ # Converting frequencies to integers for a better equality check int_freqs = list(map(int, freqs)) # If Fdrive in lookup frequencies, simply take (A, Q) slice at Fdrive index if int(Fdrive) in int_freqs: iFdrive = np.searchsorted(int_freqs, int(Fdrive)) # logger.debug('Using lookups directly at %.2f kHz', freqs[iFdrive] * 1e-3) lookups2D = {key: np.squeeze(lookups3D[key][iFdrive, :, :]) for key in lookups3D.keys()} # Otherwise, interpolate linearly between 2 (A, Q) slices at Fdrive bounding values indexes else: ilb = np.searchsorted(freqs, Fdrive) - 1 iub = ilb + 1 # logger.debug('Interpolating lookups between %.2f kHz and %.2f kHz', # freqs[ilb] * 1e-3, freqs[iub] * 1e-3) lookups2D_lb = {key: np.squeeze(lookups3D[key][ilb, :, :]) for key in lookups3D.keys()} lookups2D_ub = {key: np.squeeze(lookups3D[key][iub, :, :]) for key in lookups3D.keys()} Fratio = (Fdrive - freqs[ilb]) / (freqs[iub] - freqs[ilb]) lookups2D = {key: lookups2D_lb[key] + (lookups2D_ub[key] - lookups2D_lb[key]) * Fratio for key in lookups3D.keys()} return lookups2D 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_file = '{}_lookups_a{:.1f}nm.pkl'.format(mechname, a * 1e9) - # lookup_file = '{}_coeffs.pkl'.format(mechname) + lookup_file = 'SONIC_{}.pkl'.format(mechname) lookup_path = '{}/{}'.format(getLookupDir(), lookup_file) if not os.path.isfile(lookup_path): raise InputError('Missing lookup file: "{}"'.format(lookup_file)) # Load lookups dictionary + logger.debug('Loading lookup table') with open(lookup_path, 'rb') as fh: - # lookups4D = pickle.load(fh) - lookups3D = pickle.load(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') - Fref = lookups3D.pop('f') - Aref = lookups3D.pop('A') - Qref = lookups3D.pop('Q') + 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 InputError('Invalid sonophore diameter: {}m (must be within {}m - {}m lookup interval)' - # .format(*si_format([a, *arange], precision=2, space=' '))) + arange = (aref.min() - 1e-12, aref.max() + 1e-12) + if a < arange[0] or a > arange[1]: + raise InputError('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 InputError('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 - # lookups3D = {key: interp1d(aref, y4D, axis=0)(a) for key, y4D in lookups4D.items()} + 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 def nDindexes(dims, index): ''' Find index positions in a n-dimensional array. :param dims: dimensions of the n-dimensional array (tuple or list) :param index: index position in the flattened n-dimensional array :return: list of indexes along each array dimension ''' dims = list(dims) # Find size of each array dimension dims.reverse() dimsizes = [1] r = 1 for x in dims[:-1]: r *= x dimsizes.append(r) dims.reverse() dimsizes.reverse() # Find indexes indexes = [] remainder = index for dimsize in dimsizes[:-1]: idim, remainder = divmod(remainder, dimsize) indexes.append(idim) indexes.append(remainder) return indexes 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 checkNumBounds(values, bounds): ''' Check if a set of numbers is within predefined bounds. ''' # checking parameters against reference bounds for x, bound in zip(values, bounds): if x < bound[0] or x > bound[1]: raise ValueError('Input value {} out of [{}, {}] range'.format(x, bound[0], bound[1])) pass def getDefaultIndexes(params, defaults): ''' Return the indexes of default values found in lists of parameters. :param params: dictionary of parameter arrays :param defaults: dictionary of default values :return: dictionary of resolved default indexes ''' idefs = {} for key, default in defaults.items(): if default not in params[key]: raise Exception('default {} ({}) not found in parameter values'.format(key, default)) idefs[key] = np.where(params[key] == default)[0][0] return idefs diff --git a/sim/lookups_ASTIM.py b/sim/lookups_ASTIM.py index 1286b23..0eb9a1b 100644 --- a/sim/lookups_ASTIM.py +++ b/sim/lookups_ASTIM.py @@ -1,102 +1,102 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-02 17:50:10 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2018-09-19 15:54:52 +# @Last Modified time: 2018-09-19 19:59:55 """ Create lookup table for specific neuron. """ import os import sys import pickle import logging import numpy as np from argparse import ArgumentParser from PySONIC.solvers import computeAStimLookups, computeTestLookups from PySONIC.utils import logger, InputError, getNeuronsDict, getLookupDir # Default parameters default = { 'neuron': 'RS', 'diams': np.array([16.0, 32.0, 64.0]), # nm 'freqs': np.array([20., 100., 500., 1e3, 2e3, 3e3, 4e3]), # kHz 'amps': np.insert(np.logspace(np.log10(0.1), np.log10(600), num=50), 0, 0.0), # kPa } def main(): # Define argument parser ap = ArgumentParser() # Stimulation parameters ap.add_argument('-n', '--neuron', type=str, default=default['neuron'], help='Neuron name (string)') ap.add_argument('-a', '--diameters', nargs='+', type=float, default=None, help='Sonophore diameters (nm)') ap.add_argument('-f', '--frequencies', nargs='+', type=float, default=None, help='Acoustic drive frequencies (kHz)') ap.add_argument('-A', '--amplitudes', nargs='+', type=float, default=None, help='Acoustic pressure amplitudes (kPa)') # Boolean parameters ap.add_argument('-m', '--multiprocessing', default=False, action='store_true', help='Use multiprocessing') ap.add_argument('-v', '--verbose', default=False, action='store_true', help='Increase verbosity') ap.add_argument('-t', '--test', default=False, action='store_true', help='Run test batch') # Parse arguments args = ap.parse_args() neuron_str = args.neuron diams = (default['diams'] if args.diameters is None else np.array(args.diameters)) * 1e-9 # m freqs = (default['freqs'] if args.frequencies is None else np.array(args.frequencies)) * 1e3 # Hz amps = (default['amps'] if args.amplitudes is None else np.array(args.amplitudes)) * 1e3 # Pa if args.verbose: logger.setLevel(logging.DEBUG) else: logger.setLevel(logging.INFO) # Check neuron name validity if neuron_str not in getNeuronsDict(): raise InputError('Unknown neuron type: "{}"'.format(neuron_str)) neuron = getNeuronsDict()[neuron_str]() # Check if lookup file already exists - lookup_file = '{}_lookups.pkl'.format(neuron.name) + lookup_file = 'SONIC_{}.pkl'.format(neuron.name) lookup_filepath = '{0}/{1}'.format(getLookupDir(), lookup_file) if os.path.isfile(lookup_filepath): logger.warning('"%s" file already exists and will be overwritten. ' + 'Continue? (y/n)', lookup_file) user_str = input() if user_str not in ['y', 'Y']: logger.info('%s Lookup creation canceled', neuron.name) sys.exit(0) try: if args.test: # Compute test lookups lookup_dict = computeTestLookups(neuron, diams, freqs, amps, multiprocess=args.multiprocessing) else: # Compute real lookups lookup_dict = computeAStimLookups(neuron, diams, freqs, amps, multiprocess=args.multiprocessing) # Save dictionary in lookup file logger.info('Saving %s neuron lookup table in file: "%s"', neuron.name, lookup_file) with open(lookup_filepath, 'wb') as fh: pickle.dump(lookup_dict, fh) except InputError as err: logger.error(err) sys.exit(1) if __name__ == '__main__': main()