diff --git a/PointNICE/utils.py b/PointNICE/utils.py index ab5bc5e..4e0caae 100644 --- a/PointNICE/utils.py +++ b/PointNICE/utils.py @@ -1,254 +1,254 @@ #!/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: 2017-08-25 13:02:15 +# @Last Modified time: 2017-08-28 14:19:14 """ Definition of generic utility functions used in other modules """ from enum import Enum from functools import partial import os import logging import tkinter as tk from tkinter import filedialog from openpyxl import load_workbook import numpy as np import yaml # Get package logger logger = logging.getLogger('PointNICE') class PmCompMethod(Enum): """ Enum: types of computation method for the intermolecular pressure """ direct = 1 predict = 2 -def LoadParamsFile(filename): +def loadParamsFile(filename): """ Load a dictionary of parameters for the BLS model from an external yaml file. :param filename: name of the input file :return: parameters dictionary """ logger.info('Loading parameters from "%s"', filename) with open(filename, 'r') as f: stream = f.read() params = yaml.load(stream) return ParseNestedDict(params) -LoadParams = partial(LoadParamsFile, filename=os.path.split(__file__)[0] + '/params.yaml') +load_BLS_params = partial(loadParamsFile, filename=os.path.split(__file__)[0] + '/params.yaml') def getLookupDir(): """ Return the location of the directory holding lookups files. :return: absolute path to the directory """ this_dir, _ = os.path.split(__file__) return 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 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.get_sheet_by_name(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, c): """ 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, c): """ 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, ub, lb_new=0, ub_new=1): ''' Rescale a value to a specific interval by linear transformation. ''' xnorm = (x - lb) / (ub - lb) return xnorm * (ub_new - lb_new) + lb_new def printPct(pct, precision): print(('{:.' + str(precision) + 'f}%').format(pct), end='', flush=True) print('\r' * (precision + 3), end='') 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)) diff --git a/plot/plot_forces.py b/plot/plot_forces.py index 535b82f..c7bbce1 100644 --- a/plot/plot_forces.py +++ b/plot/plot_forces.py @@ -1,232 +1,232 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-10-07 10:22:24 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-18 15:53:50 +# @Last Modified time: 2017-08-28 14:17:56 """ Analysis of the system geometric variables and interplaying forces at stake in a static quasi-steady NICE system. """ import time import numpy as np import matplotlib.pyplot as plt import PointNICE -from PointNICE.utils import LoadParams, PmCompMethod +from PointNICE.utils import load_BLS_params, PmCompMethod plt_bool = 1 # Initialization: create a BLS instance -params = LoadParams() +params = load_BLS_params() a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} Fdrive = 0.0 # dummy stimulation frequency Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -71.9e-5 # membrane resting charge density (C/m2) bls = PointNICE.BilayerSonophore(geom, params, Fdrive, Cm0, Qm0) # Input 1: leaflet deflections ZMin = -0.45 * bls.Delta ZMax = 2 * bls.a nZ = 3000 Z = np.linspace(ZMin, ZMax, nZ) Zlb = -0.5 * bls.Delta Zub = bls.a # Input 2: acoustic perturbations PacMax = 9.5e4 nPac1 = 5 nPac2 = 100 Pac1 = np.linspace(-PacMax, PacMax, nPac1) Pac2 = np.linspace(-PacMax, PacMax, nPac2) # Input 3: membrane charge densities QmMin = bls.Qm0 QmMax = 50.0e-5 nQm = 7 Qm = np.linspace(QmMin, QmMax, nQm) # Outputs R = np.empty(nZ) Cm = np.empty(nZ) Pm_apex = np.empty(nZ) Pm_avg = np.empty(nZ) Pm_avg_predict = np.empty(nZ) Pg = np.empty(nZ) Pec = np.empty(nZ) Pel = np.empty(nZ) P0 = np.ones(nZ) * bls.P0 Pnet = np.empty(nZ) Pqs = np.empty((nZ, nPac1)) Pecdense = np.empty((nZ, nQm)) Pnetdense = np.empty((nZ, nQm)) Zeq = np.empty(nPac1) Zeq_dense = np.empty(nPac2) t0 = time.time() # Check net QS pressure at Z = 0 Peq0 = bls.PtotQS(0.0, bls.ng0, bls.Qm0, 0.0, PmCompMethod.direct) print('Net QS pressure at Z = 0.0 without perturbation: ' + '{:.2e}'.format(Peq0) + ' Pa') # Loop through the deflection vector for i in range(nZ): # 1-dimensional output vectors R[i] = bls.curvrad(Z[i]) Cm[i] = bls.Capct(Z[i]) Pm_apex[i] = bls.PMlocal(0.0, Z[i], R[i]) Pm_avg[i] = bls.PMavg(Z[i], R[i], bls.surface(Z[i])) Pm_avg_predict[i] = bls.PMavgpred(Z[i]) Pel[i] = bls.PEtot(Z[i], R[i]) Pg[i] = bls.gasmol2Pa(bls.ng0, bls.volume(Z[i])) Pec[i] = bls.Pelec(Z[i], bls.Qm0) Pnet[i] = bls.PtotQS(Z[i], bls.ng0, bls.Qm0, 0.0, PmCompMethod.direct) # loop through the acoustic perturbation vector an compute 2-dimensional # balance pressure output vector for j in range(nPac1): Pqs[i, j] = bls.PtotQS(Z[i], bls.ng0, bls.Qm0, Pac1[j], PmCompMethod.direct) for j in range(nQm): Pecdense[i, j] = bls.Pelec(Z[i], Qm[j]) Pnetdense[i, j] = bls.PtotQS(Z[i], bls.ng0, Qm[j], 0.0, PmCompMethod.direct) # Compute min local intermolecular pressure Pm_apex_min = np.amin(Pm_apex) iPm_apex_min = np.argmin(Pm_apex) print("min local intermolecular resultant pressure = %.2e Pa for z = %.2f nm" % (Pm_apex_min, Z[iPm_apex_min] * 1e9)) for j in range(nPac1): Zeq[j] = bls.balancedefQS(bls.ng0, bls.Qm0, Pac1[j], PmCompMethod.direct) for j in range(nPac2): Zeq_dense[j] = bls.balancedefQS(bls.ng0, bls.Qm0, Pac2[j], PmCompMethod.direct) t1 = time.time() print("computation completed in " + '{:.2f}'.format(t1 - t0) + " s") if plt_bool == 1: # 1: Intermolecular pressures fig1, ax = plt.subplots() fig1.canvas.set_window_title("1: integrated vs. predicted average intermolecular pressure") ax.set_xlabel('Z $(nm)$', fontsize=18) ax.set_ylabel('Pressures $(MPa)$', fontsize=18) ax.grid(True) ax.plot([Zlb * 1e9, Zlb * 1e9], [np.amin(Pm_avg) * 1e-6, np.amax(Pm_avg) * 1e-6], '--', color="blue", label="$-\Delta /2$") ax.plot([Zub * 1e9, Zub * 1e9], [np.amin(Pm_avg) * 1e-6, np.amax(Pm_avg) * 1e-6], '--', color="red", label="$a$") ax.plot(Z * 1e9, Pm_avg * 1e-6, '-', label="$P_{M, avg}$", color="green", linewidth=2.0) ax.plot(Z * 1e9, Pm_avg_predict * 1e-6, '-', label="$P_{M, avg-predict}$", color="red", linewidth=2.0) ax.set_xlim(ZMin * 1e9 - 5, ZMax * 1e9) ax.legend(fontsize=24) # 2: Capacitance and electric pressure fig2, ax = plt.subplots() fig2.canvas.set_window_title("2: Capacitance and electric equivalent pressure") ax.set_xlabel('Z $(nm)$', fontsize=18) ax.set_ylabel('$C_m \ (uF/cm^2)$', fontsize=18) ax.plot(Z * 1e9, Cm * 1e2, '-', label="$C_{m}$", color="black", linewidth=2.0) ax.set_xlim(ZMin * 1e9 - 5, ZMax * 1e9) ax2 = ax.twinx() ax2.set_ylabel('$P_{EC}\ (MPa)$', fontsize=18, color='magenta') ax2.plot(Z * 1e9, Pec * 1e-6, '-', label="$P_{EC}$", color="magenta", linewidth=2.0) # tmp: electric pressure for varying membrane charge densities figtmp, ax = plt.subplots() figtmp.canvas.set_window_title("electric pressure for varying membrane charges") ax.set_xlabel('Z $(nm)$', fontsize=18) ax.set_ylabel('$P_{EC} \ (MPa)$', fontsize=18) for j in range(nQm): lbl = "$Q_m$ = " + '{:.2f}'.format(Qm[j] * 1e5) + " nC/cm2" ax.plot(Z * 1e9, Pecdense[:, j] * 1e-6, '-', label=lbl, linewidth=2.0) ax.set_xlim(ZMin * 1e9 - 5, ZMax * 1e9) ax.legend() # tmp: net pressure for varying membrane potentials figtmp, ax = plt.subplots() figtmp.canvas.set_window_title("net pressure for varying membrane charges") ax.set_xlabel('Z $(nm)$', fontsize=18) ax.set_ylabel('$P_{net} \ (MPa)$', fontsize=18) for j in range(nQm): lbl = "$Q_m$ = " + '{:.2f}'.format(Qm[j] * 1e5) + " nC/cm2" ax.plot(Z * 1e9, Pnetdense[:, j] * 1e-6, '-', label=lbl, linewidth=2.0) ax.set_xlim(ZMin * 1e9 - 5, ZMax * 1e9) ax.legend() # 3: Net pressure without perturbation fig3, ax = plt.subplots() fig3.canvas.set_window_title("3: Net QS pressure without perturbation") ax.set_xlabel('Z $(nm)$', fontsize=18) ax.set_ylabel('Pressures $(kPa)$', fontsize=18) # ax.grid(True) # ax.plot([Zlb * 1e9, Zlb * 1e9], [np.amin(Pec) * 1e-3, np.amax(Pm_avg) * 1e-3], '--', # color="blue", label="$-\Delta / 2$") # ax.plot([Zub * 1e9, Zub * 1e9], [np.amin(Pec) * 1e-3, np.amax(Pm_avg) * 1e-3], '--', # color="red", label="$a$") ax.plot(Z * 1e9, Pg * 1e-3, '-', label="$P_{gas}$", linewidth=3.0, color='C0') ax.plot(Z * 1e9, -P0 * 1e-3, '-', label="$-P_{0}$", linewidth=3.0, color='C1') ax.plot(Z * 1e9, Pm_avg * 1e-3, '-', label="$P_{mol}$", linewidth=3.0, color='C2') ax.plot(Z * 1e9, Pec * 1e-3, '-', label="$P_{elec}$", linewidth=3.0, color='C3') ax.plot(Z * 1e9, Pel * 1e-3, '-', label="$P_{elastic}$", linewidth=3.0, color='C4') # ax.plot(Z * 1e9, (Pg - P0 + Pm_avg + Pec + Pel) * 1e-3, '--', label="$P_{net}$", linewidth=2.0, # color='black') # ax.plot(Z * 1e9, (Pg - P0 + Pm_avg + Pec - Pnet) * 1e-6, '--', label="$P_{net} diff$", # linewidth=2.0, color="blue") ax.set_xlim(ZMin * 1e9 - 5, 30) ax.set_ylim(-1500, 2000) ax.legend(fontsize=24) # ax.grid(True) # 4: QS pressure for different perturbations fig4, ax = plt.subplots() fig4.canvas.set_window_title("4: Net QS pressure for different acoustic perturbations") ax.set_xlabel('Z $(nm)$', fontsize=18) ax.set_ylabel('Pressures $(MPa)$', fontsize=18) ax.grid(True) ax.plot([Zlb * 1e9, Zlb * 1e9], [np.amin(Pqs[:, 0]) * 1e-6, np.amax(Pqs[:, nPac1 - 1]) * 1e-6], '--', color="blue", label="$-\Delta/2$") ax.plot([Zub * 1e9, Zub * 1e9], [np.amin(Pqs[:, 0]) * 1e-6, np.amax(Pqs[:, nPac1 - 1]) * 1e-6], '--', color="red", label="$a$") ax.set_xlim(ZMin * 1e9 - 5, ZMax * 1e9) for j in range(nPac1): lbl = "$P_{A}$ = %.2f MPa" % (Pac1[j] * 1e-6) ax.plot(Z * 1e9, Pqs[:, j] * 1e-6, '-', label=lbl, linewidth=2.0) ax.plot([Zeq[j] * 1e9, Zeq[j] * 1e9], [np.amin(Pqs[:, nPac1 - 1]) * 1e-6, np.amax(Pqs[:, 0]) * 1e-6], '--', color="black") ax.legend(fontsize=24) # 5: QS balance deflection for different acoustic perturbations fig5, ax = plt.subplots() fig5.canvas.set_window_title("5: QS balance deflection for different acoustic perturbations ") ax.set_xlabel('Perturbation $(MPa)$', fontsize=18) ax.set_ylabel('Z $(nm)$', fontsize=18) ax.plot([np.amin(Pac2) * 1e-6, np.amax(Pac2) * 1e-6], [Zlb * 1e9, Zlb * 1e9], '--', color="blue", label="$-\Delta / 2$") ax.plot([np.amin(Pac2) * 1e-6, np.amax(Pac2) * 1e-6], [Zub * 1e9, Zub * 1e9], '--', color="red", label="$a$") ax.plot([-bls.P0 * 1e-6, -bls.P0 * 1e-6], [np.amin(Zeq_dense) * 1e9, np.amax(Zeq_dense) * 1e9], '--', color="black", label="$-P_0$") ax.plot(Pac2 * 1e-6, Zeq_dense * 1e9, '-', label="$Z_{eq}$", linewidth=2.0) ax.set_xlim(-0.12, 0.12) ax.set_ylim(ZMin * 1e9 - 5, bls.a * 1e9 + 5) ax.legend(fontsize=24) plt.show() diff --git a/postpro/postpro_latency.py b/postpro/postpro_latency.py index caa879e..cdab810 100644 --- a/postpro/postpro_latency.py +++ b/postpro/postpro_latency.py @@ -1,65 +1,65 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-10-31 10:10:41 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-18 15:39:46 +# @Last Modified time: 2017-08-28 14:18:34 """ Test relationship between stimulus intensity and response latency. """ import numpy as np import matplotlib.pyplot as plt -from PointNICE.utils import ImportExcelCol, LoadParams +from PointNICE.utils import ImportExcelCol, load_BLS_params # Load NICE parameters -params = LoadParams() +params = load_BLS_params() biomech = params['biomech'] ac_imp = biomech['rhoL'] * biomech['c'] # Rayl # Define import settings xls_file = "C:/Users/admin/Desktop/Model output/NBLS spikes 0.35MHz/nbls_log_spikes_0.35MHz.xlsx" sheet = 'Data' # Import data f = ImportExcelCol(xls_file, sheet, 'E', 2) * 1e3 # Hz A = ImportExcelCol(xls_file, sheet, 'F', 2) * 1e3 # Pa T = ImportExcelCol(xls_file, sheet, 'G', 2) * 1e-3 # s N = ImportExcelCol(xls_file, sheet, 'Q', 2) L = ImportExcelCol(xls_file, sheet, 'R', 2) # ms # Retrieve unique values of latencies (for min. 2 spikes) and corresponding amplitudes iremove = np.where(N < 2)[0] A_true = np.delete(A, iremove) L_true = np.delete(L, iremove).astype(np.float) latencies, indices = np.unique(L_true, return_index=True) amplitudes = A_true[indices] # Convert amplitudes to intensities intensities = amplitudes**2 / (2 * ac_imp) * 1e-4 # W/cm2 # Plot latency vs. amplitude fig1, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$Amplitude \ (kPa)$", fontsize=28) ax.set_ylabel("$Latency \ (ms)$", fontsize=28) ax.scatter(amplitudes * 1e-3, latencies, color='black', s=100) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) # Plot latency vs. intensity fig2, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$Intensity \ (W/cm^2)$", fontsize=28) ax.set_ylabel("$Latency \ (ms)$", fontsize=28) ax.scatter(intensities, latencies, color='black', s=100) ax.set_xticks([0, 0.2, 0.4, 0.6, 0.8]) ax.set_yticks([25, 35, 55, 65]) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) plt.show() diff --git a/postpro/postpro_spikerate.py b/postpro/postpro_spikerate.py index ce2446d..a2ba918 100644 --- a/postpro/postpro_spikerate.py +++ b/postpro/postpro_spikerate.py @@ -1,83 +1,83 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-10-31 11:27:34 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-18 15:46:49 +# @Last Modified time: 2017-08-28 14:18:32 """ Test relationship between stimulus intensity spike rate. """ import numpy as np from scipy.optimize import curve_fit import matplotlib.pyplot as plt -from PointNICE.utils import ImportExcelCol, LoadParams +from PointNICE.utils import ImportExcelCol, load_BLS_params def fitfunc(x, a, b): """ Fitting function """ return a * np.power(x, b) # Load NICE parameters -params = LoadParams() +params = load_BLS_params() biomech = params['biomech'] ac_imp = biomech['rhoL'] * biomech['c'] # Rayl # Import data xls_file = "C:/Users/admin/Desktop/Model output/NBLS spikes 0.35MHz/nbls_log_spikes_0.35MHz.xlsx" sheet = 'Data' f = ImportExcelCol(xls_file, sheet, 'E', 2) * 1e3 # Hz A = ImportExcelCol(xls_file, sheet, 'F', 2) * 1e3 # Pa T = ImportExcelCol(xls_file, sheet, 'G', 2) * 1e-3 # s N = ImportExcelCol(xls_file, sheet, 'Q', 2) FR = ImportExcelCol(xls_file, sheet, 'S', 2) # ms # Retrieve available spike rates values (for min. 3 spikes) and corresponding amplitudes iremove = np.where(N < 15)[0] A_true = np.delete(A, iremove) spikerates = np.delete(FR, iremove).astype(np.float) amplitudes = np.delete(A, iremove) # Convert amplitudes to intensities intensities = amplitudes**2 / (2 * ac_imp) * 1e-4 # W/cm2 # Power law least square fitting popt, pcov = curve_fit(fitfunc, intensities, spikerates) print('power product fit: FR = %.2f I^%.2f' % (popt[0], popt[1])) # Compute predicted data and associated error spikerates_predicted = fitfunc(intensities, popt[0], popt[1]) residuals = spikerates - spikerates_predicted ss_res = np.sum(residuals**2) ss_tot = np.sum((spikerates - np.mean(spikerates))**2) r_squared = 1 - (ss_res / ss_tot) print("R-squared = " + '{:.5f}'.format(r_squared)) # Plot latency vs. amplitude fig1, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$Amplitude \ (kPa)$", fontsize=28) ax.set_ylabel("$Spike\ Rate \ (spikes/ms)$", fontsize=28) ax.scatter(amplitudes * 1e-3, spikerates, color='black') ax.set_ylim(0, 1.1 * np.amax(spikerates)) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) # Plot latency vs. intensity fig2, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$Intensity \ (W/cm^2)$", fontsize=28) ax.set_ylabel("$Spike\ Rate \ (spikes/ms)$", fontsize=28) ax.scatter(intensities, spikerates, color='black', label='$data$') ax.plot(intensities, spikerates_predicted, color='blue', label='$%.2f\ I^{%.2f}$' % (popt[0], popt[1])) ax.set_ylim(0, 1.1 * np.amax(spikerates)) ax.legend(fontsize=28) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) plt.show() diff --git a/postpro/postpro_spikes.py b/postpro/postpro_spikes.py index 78fad20..b5aedab 100644 --- a/postpro/postpro_spikes.py +++ b/postpro/postpro_spikes.py @@ -1,126 +1,126 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-10-27 09:50:55 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-18 15:48:58 +# @Last Modified time: 2017-08-28 14:18:29 """ Test influence of acoustic intensity and duration on number of spikes. """ import numpy as np import matplotlib.pyplot as plt import matplotlib.cm as cm from mpl_toolkits.mplot3d import Axes3D -from PointNICE.utils import ImportExcelCol, ConstructMatrix, LoadParams +from PointNICE.utils import ImportExcelCol, ConstructMatrix, load_BLS_params # Define options plot2d_bool = 0 plot3d_show = 1 plot3d_save = 0 plt_root = "../Output/effective spikes 2D/" plt_save_ext = '.png' # Load NICE parameters -params = LoadParams() +params = load_BLS_params() biomech = params['biomech'] ac_imp = biomech['rhoL'] * biomech['c'] # Rayl # Import data xls_file = "../../Output/effective spikes 2D/nbls_log_spikes.xlsx" sheet = 'Data' f_all = ImportExcelCol(xls_file, sheet, 'E', 2) * 1e3 # Hz A_all = ImportExcelCol(xls_file, sheet, 'F', 2) * 1e3 # Pa T_all = ImportExcelCol(xls_file, sheet, 'G', 2) * 1e-3 # s N_all = ImportExcelCol(xls_file, sheet, 'Q', 2) # number of spikes freqs = np.unique(f_all) for Fdrive in freqs: # Select data A = A_all[f_all == Fdrive] T = T_all[f_all == Fdrive] N = N_all[f_all == Fdrive] # Reshape serialized data into 2 dimensions (durations, amps, nspikes, nholes) = ConstructMatrix(T, A, N) nspikes2 = nspikes.conj().T # conjugate tranpose of nspikes matrix (for surface plot) # Convert to appropriate units intensities = amps**2 / (2 * ac_imp) * 1e-4 # W/cm2 durations = durations * 1e3 # ms nDurations = durations.size nIntensities = intensities.size Tmax = np.amax(durations) Tmin = np.amin(durations) Imax = np.amax(intensities) Imin = np.amin(intensities) print(str(nholes) + " hole(s) in reconstructed matrix") mymap = cm.get_cmap('jet') if plot2d_bool == 1: # Plot spikes vs. intensity (with duration color code) fig, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$I \ (W/cm^2)$", fontsize=28) ax.set_ylabel("$\#\ spikes$", fontsize=28) for i in range(nIntensities): ax.plot(intensities, nspikes[i, :], c=mymap((durations[i] - Tmin) / (Tmax - Tmin)), label='t = ' + str(durations[i]) + ' ms') sm_duration = plt.cm.ScalarMappable(cmap=mymap, norm=plt.Normalize(Tmin, Tmax)) sm_duration._A = [] cbar = plt.colorbar(sm_duration) cbar.ax.set_ylabel('$duration \ (ms)$', fontsize=28) # Plot spikes vs. duration (with intensity color code) fig, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$duration \ (ms)$", fontsize=28) ax.set_ylabel("$\#\ spikes$", fontsize=28) for j in range(nDurations): ax.plot(durations, nspikes[:, j], c=mymap((intensities[j] - Imin) / (Imax - Imin)), label='I = ' + str(intensities[j]) + ' W/cm2') sm_int = plt.cm.ScalarMappable(cmap=mymap, norm=plt.Normalize(Imin, Imax)) sm_int._A = [] cbar = plt.colorbar(sm_int) cbar.ax.set_ylabel("$I \ (W/cm^2)$", fontsize=28) if plot3d_show == 1 and nholes == 0: # 3D surface plot: nspikes = f(duration, intensity) X, Y = np.meshgrid(durations, intensities) fig = plt.figure(figsize=(12, 9)) ax = fig.gca(projection=Axes3D.name) ax.plot_surface(X, Y, nspikes2, rstride=1, cstride=1, cmap=mymap, linewidth=0, antialiased=False) ax.set_xlabel("$duration \ (ms)$", fontsize=24, labelpad=20) ax.set_ylabel("$intensity \ (W/cm^2)$", fontsize=24, labelpad=20) ax.set_zlabel("$\#\ spikes$", fontsize=24, labelpad=20) csetx = ax.contour(X, Y, nspikes2, zdir='x', offset=150, cmap=cm.coolwarm) csety = ax.contour(X, Y, nspikes2, zdir='y', offset=0.8, cmap=cm.coolwarm) ax.view_init(33, -126) ax.set_xticks([0, 50, 100, 150]) ax.set_yticks([0, 0.2, 0.4, 0.6, 0.8]) ax.set_zticks([0, 20, 40, 60, 80]) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) for item in ax.get_zticklabels(): item.set_fontsize(24) # Save figure if needed if plot3d_save == 1: plt_filename = '{}spikes_{:.0f}KHz{}'.format(plt_root, Fdrive * 1e-3, plt_save_ext) plt.savefig(plt_filename) print('Saving figure to "' + plt_root + '"') plt.close() plt.show() diff --git a/postpro/postpro_threshold_duration.py b/postpro/postpro_threshold_duration.py index 8c51987..9e1de31 100644 --- a/postpro/postpro_threshold_duration.py +++ b/postpro/postpro_threshold_duration.py @@ -1,72 +1,72 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-10-30 21:48:45 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-18 15:49:11 +# @Last Modified time: 2017-08-28 14:18:27 """ Test relationship between stimulus duration and minimum acoustic amplitude / intensity / energy for AP generation. """ import numpy as np import matplotlib.pyplot as plt -from PointNICE.utils import ImportExcelCol, LoadParams +from PointNICE.utils import ImportExcelCol, load_BLS_params # Load NICE parameters -params = LoadParams() +params = load_BLS_params() biomech = params['biomech'] ac_imp = biomech['rhoL'] * biomech['c'] # Rayl # Import data xls_file = "C:/Users/admin/Desktop/Model output/NBLS titration duration 0.35MHz/nbls_log_titration_duration_0.35MHz.xlsx" sheet = 'Data' f = ImportExcelCol(xls_file, sheet, 'E', 2) * 1e3 # Hz A = ImportExcelCol(xls_file, sheet, 'F', 2) * 1e3 # Pa T = ImportExcelCol(xls_file, sheet, 'G', 2) * 1e-3 # s N = ImportExcelCol(xls_file, sheet, 'Q', 2) # Convert to appropriate units durations = T * 1e3 # ms Trange = np.amax(durations) - np.amin(durations) amplitudes = A * 1e-3 # kPa intensities = A**2 / (2 * ac_imp) * 1e-4 # W/cm2 energies = intensities * durations # mJ/cm2 # Plot threshold amplitude vs. duration fig1, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$duration \ (ms)$", fontsize=28) ax.set_ylabel("$Amplitude \ (kPa)$", fontsize=28) ax.scatter(durations, amplitudes, color='black', s=100) ax.set_xlim(np.amin(durations) - 0.1 * Trange, np.amax(durations) + 0.1 * Trange) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) # Plot threshold intensity vs. duration fig2, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$duration \ (ms)$", fontsize=28) ax.set_ylabel("$Intensity \ (W/cm^2)$", fontsize=28) ax.scatter(durations, intensities, color='black', s=100) ax.set_xlim(np.amin(durations) - 0.1 * Trange, np.amax(durations) + 0.1 * Trange) ax.set_yticks([np.floor(np.amin(intensities) * 1e2) / 1e2, np.ceil(np.amax(intensities) * 1e2) / 1e2]) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) # Plot threshold energy vs. duration fig3, ax = plt.subplots(figsize=(12, 9)) ax.set_xlabel("$duration \ (ms)$", fontsize=28) ax.set_ylabel("$Energy \ (mJ/cm^2)$", fontsize=28) ax.scatter(durations, energies, color='black', s=100) ax.set_xlim(np.amin(durations) - 0.1 * Trange, np.amax(durations) + 0.1 * Trange) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) plt.show() diff --git a/postpro/postpro_threshold_frequency.py b/postpro/postpro_threshold_frequency.py index dc5e348..603c0ba 100644 --- a/postpro/postpro_threshold_frequency.py +++ b/postpro/postpro_threshold_frequency.py @@ -1,69 +1,69 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2016-10-30 21:48:45 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-07-18 15:49:20 +# @Last Modified time: 2017-08-28 14:18:23 """ Test relationship between stimulus frequency and minimum acoustic intensity for AP generation. """ import numpy as np import matplotlib.pyplot as plt import matplotlib.ticker as ticker -from PointNICE.utils import ImportExcelCol, LoadParams +from PointNICE.utils import ImportExcelCol, load_BLS_params # Load NICE parameters -params = LoadParams() +params = load_BLS_params() biomech = params['biomech'] ac_imp = biomech['rhoL'] * biomech['c'] # Rayl # Import data xls_file = "C:/Users/admin/Desktop/Model output/NBLS titration frequency 30ms/nbls_log_titration_frequency_30ms.xlsx" sheet = 'Data' f = ImportExcelCol(xls_file, sheet, 'E', 2) * 1e3 # Hz A = ImportExcelCol(xls_file, sheet, 'F', 2) * 1e3 # Pa T = ImportExcelCol(xls_file, sheet, 'G', 2) * 1e-3 # s N = ImportExcelCol(xls_file, sheet, 'Q', 2) # Convert to appropriate units frequencies = f * 1e-6 # MHz amplitudes = A * 1e-3 # kPa intensities = A**2 / (2 * ac_imp) * 1e-4 # W/cm2 # Plot threshold amplitude vs. duration fig1, ax = plt.subplots(figsize=(12, 9)) ax.set_xscale('log') ax.set_xlabel("$Frequency \ (MHz)$", fontsize=28) ax.set_ylabel("$Amplitude \ (kPa)$", fontsize=28) ax.scatter(frequencies, amplitudes, color='black', s=100) ax.set_xlim(1.5e-1, 5e0) ax.set_xscale('log') ax.set_xticks([0.2, 1, 4]) ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter()) ax.set_yticks([np.floor(np.amin(amplitudes)), np.ceil(np.amax(amplitudes))]) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) # Plot threshold intensity vs. duration fig2, ax = plt.subplots(figsize=(12, 9)) ax.set_xscale('log') ax.set_xlabel("$Frequency \ (MHz)$", fontsize=28) ax.set_ylabel("$Intensity \ (W/cm^2)$", fontsize=28) ax.scatter(frequencies, intensities, color='black', s=100) ax.set_xlim(1.5e-1, 5e0) ax.set_xscale('log') ax.set_xticks([0.2, 1, 4]) ax.get_xaxis().set_major_formatter(ticker.ScalarFormatter()) ax.set_yticks([np.floor(np.amin(intensities) * 1e2) / 1e2, np.ceil(np.amax(intensities) * 1e2) / 1e2]) for item in ax.get_yticklabels(): item.set_fontsize(24) for item in ax.get_xticklabels(): item.set_fontsize(24) plt.show() diff --git a/sim/ASTIM_batch.py b/sim/ASTIM_batch.py index d62728b..5b6b98a 100644 --- a/sim/ASTIM_batch.py +++ b/sim/ASTIM_batch.py @@ -1,57 +1,57 @@ #!/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: 2017-08-25 18:23:59 +# @Last Modified time: 2017-08-28 14:18:21 """ Run batch acoustic simulations of specific "point-neuron" models. """ import os import logging import numpy as np from PointNICE.solvers import checkBatchLog, runAStimBatch from PointNICE.channels import * -from PointNICE.utils import LoadParams +from PointNICE.utils import load_BLS_params from PointNICE.plt import plotBatch # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # BLS parameters -bls_params = LoadParams() +bls_params = load_BLS_params() # Geometry of BLS structure a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} # Channels mechanisms neurons = [CorticalLTS()] # Stimulation parameters stim_params = { 'freqs': [3.5e5], # Hz 'amps': [150e3], # Pa 'durations': [50e-3], # s 'PRFs': [1e2], # Hz 'DFs': [0.5] } stim_params['offsets'] = [30e-3] * len(stim_params['durations']) # s # Select output directory try: (batch_dir, log_filepath) = checkBatchLog('A-STIM') except AssertionError as err: logger.error(err) quit() # Run A-STIM batch pkl_filepaths = runAStimBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles plotBatch(pkl_dir, pkl_filepaths) diff --git a/sim/ASTIM_lookups.py b/sim/ASTIM_lookups.py index 9004833..ba0f0ae 100644 --- a/sim/ASTIM_lookups.py +++ b/sim/ASTIM_lookups.py @@ -1,50 +1,50 @@ #!/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: 2017-08-22 14:30:33 +# @Last Modified time: 2017-08-28 14:18:17 """ Create lookup tables for different acoustic frequencies. """ import logging import numpy as np import PointNICE -from PointNICE.utils import LoadParams +from PointNICE.utils import load_BLS_params from PointNICE.channels import * # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # BLS parameters -params = LoadParams() +params = load_BLS_params() # Geometry of NBLS structure a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} # Channel mechanisms neurons = [ThalamoCortical()] # Stimulation parameters freqs = [690e3] # Hz amps = np.logspace(np.log10(0.1), np.log10(600), num=50) * 1e3 # Pa amps = np.insert(amps, 0, 0.0) # adding amplitude 0 logger.info('Starting batch lookup creation') for ch_mech in neurons: # Create a SolverUS instance (with dummy frequency parameter) solver = PointNICE.SolverUS(geom, params, ch_mech, 0.0) charges = np.arange(np.round(ch_mech.Vm0 - 10.0), 50.0 + 1.0, 1.0) * 1e-5 # C/m2 # Create lookups for each frequency for Fdrive in freqs: solver.createLookup(ch_mech, Fdrive, amps, charges) logger.info('Lookup tables successfully created') diff --git a/sim/ASTIM_titration_batch.py b/sim/ASTIM_titration_batch.py index 7440003..91a8604 100644 --- a/sim/ASTIM_titration_batch.py +++ b/sim/ASTIM_titration_batch.py @@ -1,56 +1,56 @@ #!/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: 2017-08-25 18:25:39 +# @Last Modified time: 2017-08-28 14:18:52 """ Run batch acoustic titrations of specific "point-neuron" models. """ import os import logging import numpy as np from PointNICE.solvers import checkBatchLog, titrateAStimBatch from PointNICE.channels import * -from PointNICE.utils import LoadParams +from PointNICE.utils import load_BLS_params from PointNICE.plt import plotBatch # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # BLS parameters -bls_params = LoadParams() +bls_params = load_BLS_params() # Geometry of NBLS structure a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} # Channels mechanisms neurons = [CorticalRS()] # Stimulation parameters stim_params = { 'freqs': [3.5e5], # Hz 'amps': [100e3], # Pa 'durations': [50e-3], # s 'PRFs': [1e2], # Hz # 'DFs': [1.0] } # Select output directory try: (batch_dir, log_filepath) = checkBatchLog('A-STIM') except AssertionError as err: logger.error(err) quit() # Run titration batch pkl_filepaths = titrateAStimBatch(batch_dir, log_filepath, neurons, bls_params, geom, stim_params) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles plotBatch({'Q_m': ['Qm']}, pkl_dir, pkl_filepaths) diff --git a/sim/MECH_batch.py b/sim/MECH_batch.py index 8b3809e..1bc585d 100644 --- a/sim/MECH_batch.py +++ b/sim/MECH_batch.py @@ -1,55 +1,55 @@ #!/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: 2017-08-28 10:23:04 +# @Last Modified time: 2017-08-28 14:18:54 """ Run batch simulations of the NICE mechanical model with imposed charge densities """ import os import logging import numpy as np -from PointNICE.utils import LoadParams +from PointNICE.utils import load_BLS_params from PointNICE.solvers import checkBatchLog, runMechBatch from PointNICE.plt import plotBatch # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # BLS parameters -bls_params = LoadParams() +bls_params = load_BLS_params() # Geometry of BLS structure a = 32e-9 # in-plane radius (m) d = 0.0e-6 # embedding tissue thickness (m) geom = {"a": a, "d": d} # Electrical properties of the membrane Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -80e-5 # membrane resting charge density (C/m2) # Stimulation parameters stim_params = { 'freqs': [3.5e5], # Hz 'amps': [100e3], # Pa 'charges': [80e-5] # C/m2 } # Select output directory try: (batch_dir, log_filepath) = checkBatchLog('mech') except AssertionError as err: logger.error(err) quit() # Run mechanical batch pkl_filepaths = runMechBatch(batch_dir, log_filepath, bls_params, geom, Cm0, Qm0, stim_params) pkl_dir, _ = os.path.split(pkl_filepaths[0]) # Plot resulting profiles plotBatch(pkl_dir, pkl_filepaths) diff --git a/tests/test_basic.py b/tests/test_basic.py index 94f9754..a65623d 100644 --- a/tests/test_basic.py +++ b/tests/test_basic.py @@ -1,271 +1,271 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-14 18:37:45 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-28 14:07:37 +# @Last Modified time: 2017-08-28 14:18:56 ''' Test the basic functionalities of the package. ''' import logging import inspect import numpy as np -from PointNICE.utils import LoadParams +from PointNICE.utils import load_BLS_params from PointNICE import BilayerSonophore, channels, SolverElec, SolverUS from PointNICE.solvers import detectSpikes, titrateEStim, titrateAStim from PointNICE.constants import * # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # List of implemented neurons neurons = [] for _, obj in inspect.getmembers(channels): if inspect.isclass(obj) and isinstance(obj.name, str): if obj.name != 'LeechT': neurons.append(obj) def test_MECH(): ''' Maximal negative and positive deflections of the BLS structure for a specific sonophore size, resting membrane properties and stimulation parameters. ''' logger.info('Starting test: Mechanical simulation') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} - params = LoadParams() + params = load_BLS_params() # Create BLS instance Fdrive = 350e3 # Hz Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -80e-5 # membrane resting charge density (C/m2) bls = BilayerSonophore(geom, params, Fdrive, Cm0, Qm0) # Run mechanical simulation Adrive = 100e3 # Pa Qm = 50e-5 # C/m2 bls.runMech(Fdrive, Adrive, Qm) _, y, _ = bls.runMech(Fdrive, Adrive, Qm) # Check validity of deflection extrema Z, _ = y Zlast = Z[-NPC_FULL:] Zmin, Zmax = (Zlast.min(), Zlast.max()) logger.info('Zmin = %.2f nm, Zmax = %.2f nm', Zmin * 1e9, Zmax * 1e9) Zmin_ref, Zmax_ref = (-0.116e-9, 5.741e-9) assert np.abs(Zmin - Zmin_ref) < 1e-12, 'Unexpected BLS compression amplitude' assert np.abs(Zmax - Zmax_ref) < 1e-12, 'Unexpected BLS expansion amplitude' logger.info('Passed test: Mechanical simulation') def test_resting_potential(): ''' Neurons membrane potential in free conditions should stabilize to their specified resting potential value. ''' conv_err_msg = ('{} neuron membrane potential in free conditions does not converge to ' 'stable value (gap after 20s: {:.2e} mV)') value_err_msg = ('{} neuron steady-state membrane potential in free conditions differs ' 'significantly from specified resting potential (gap = {:.2f} mV)') logger.info('Starting test: neurons in free conditions') # Initialize solver solver = SolverElec() for neuron_class in neurons: # Simulate each neuron in free conditions neuron = neuron_class() _, y, _ = solver.run(neuron, Astim=0.0, tstim=20.0, toffset=0.0) Vm_free, *_ = y # Check membrane potential convergence Vm_free_last, Vm_free_beforelast = (Vm_free[-1], Vm_free[-2]) Vm_free_conv = Vm_free_last - Vm_free_beforelast assert np.abs(Vm_free_conv) < 1e-5, conv_err_msg.format(neuron.name, Vm_free_conv) # Check membrane potential convergence to resting potential Vm_free_diff = Vm_free_last - neuron.Vm0 assert np.abs(Vm_free_diff) < 0.1, value_err_msg.format(neuron.name, Vm_free_diff) logger.info('Passed test: neurons in free conditions') def test_ESTIM(): ''' Threshold E-STIM amplitude and needed to obtain an action potential and response latency should match reference values. ''' Athr_err_msg = ('{} neuron threshold amplitude for excitation does not match reference value' '(gap = {:.2f} mA/m2)') latency_err_msg = ('{} neuron latency for excitation at threshold amplitude does not match ' 'reference value (gap = {:.2f} ms)') logger.info('Starting test: E-STIM titration') # Initialize solver solver = SolverElec() Arange = (0.0, 2 * TITRATION_ESTIM_A_MAX) # mA/m2 # Stimulation parameters tstim = 100e-3 # s toffset = 50e-3 # s # Reference values Athr_refs = {'FS': 6.91, 'LTS': 1.54, 'RS': 5.03, 'LeechT': 5.54, 'RE': 3.61, 'TC': 4.05} latency_refs = {'FS': 101.00e-3, 'LTS': 128.56e-3, 'RS': 103.81e-3, 'LeechT': 20.22e-3, 'RE': 148.50e-3, 'TC': 63.46e-3} for neuron_class in neurons: # Perform titration for each neuron neuron = neuron_class() print(neuron.name, 'neuron titration') (Athr, t, y, _, latency) = titrateEStim(solver, neuron, Arange, tstim, toffset, PRF=None, DF=1.0) Vm = y[0, :] # Check that final number of spikes is 1 n_spikes, _, _ = detectSpikes(t, Vm, SPIKE_MIN_VAMP, SPIKE_MIN_DT) assert n_spikes == 1, 'Number of spikes after titration should be exactly 1' # Check threshold amplitude Athr_diff = Athr - Athr_refs[neuron.name] assert np.abs(Athr_diff) < 0.1, Athr_err_msg.format(neuron.name, Athr_diff) # Check response latency lat_diff = (latency - latency_refs[neuron.name]) * 1e3 assert np.abs(lat_diff) < 1.0, latency_err_msg.format(neuron.name, lat_diff) logger.info('Passed test: E-STIM titration') def test_ASTIM(): ''' Threshold A-STIM amplitude and needed to obtain an action potential and response latency should match reference values. ''' Athr_err_msg = ('{} neuron threshold amplitude for excitation does not match reference value' '(gap = {:.2f} kPa)') latency_err_msg = ('{} neuron latency for excitation at threshold amplitude does not match ' 'reference value (gap = {:.2f} ms)') logger.info('Starting test: A-STIM titration') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} - params = LoadParams() + params = load_BLS_params() # Stimulation parameters Fdrive = 350e3 # Hz tstim = 50e-3 # s toffset = 30e-3 # s Arange = (0.0, 2 * TITRATION_ASTIM_A_MAX) # Pa # Reference values Athr_refs = {'FS': 38.67e3, 'LTS': 24.80e3, 'RS': 51.17e3, 'RE': 46.36e3, 'TC': 23.24e3, 'LeechT': None} latency_refs = {'FS': 63.72e-3, 'LTS': 61.92e-3, 'RS': 62.52e-3, 'RE': 79.20e-3, 'TC': 68.53e-3, 'LeechT': None} # Titration for each neuron for neuron_class in neurons: # Initialize neuron neuron = neuron_class() print(neuron.name, 'neuron titration') # Initialize solver solver = SolverUS(geom, params, neuron, Fdrive) # Perform titration (Athr, t, y, _, latency) = titrateAStim(solver, neuron, Fdrive, Arange, tstim, toffset, PRF=None, DF=1.0) Qm = y[2] # Check that final number of spikes is 1 n_spikes, _, _ = detectSpikes(t, Qm, SPIKE_MIN_QAMP, SPIKE_MIN_DT) assert n_spikes == 1, 'Number of spikes after titration should be exactly 1' # Check threshold amplitude Athr_diff = (Athr - Athr_refs[neuron.name]) * 1e-3 assert np.abs(Athr_diff) < 0.1, Athr_err_msg.format(neuron.name, Athr_diff) # Check response latency lat_diff = (latency - latency_refs[neuron.name]) * 1e3 assert np.abs(lat_diff) < 1.0, latency_err_msg.format(neuron.name, lat_diff) logger.info('Passed test: A-STIM titration') def test_ASTIM_classic(): ''' Short classic A-STIM simulation should not raise any errors. ''' logger.info('Starting test: A-STIM classic simulation') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} - params = LoadParams() + params = load_BLS_params() # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa tstim = 1e-3 # s toffset = 1e-3 # s # Initialize RS neuron rs_neuron = channels.CorticalRS() # Initialize solver solver = SolverUS(geom, params, rs_neuron, Fdrive) # Run short classic simulation of the system solver.run(rs_neuron, Fdrive, Adrive, tstim, toffset, sim_type='classic') # If no error is raised, test is passed successfully logger.info('Passed test: A-STIM classic simulation') def test_ASTIM_hybrid(): ''' Short hybrid A-STIM simulation should not raise any errors. ''' logger.info('Starting test: A-STIM hybrid simulation') # BLS geometry and parameters geom = {"a": 32e-9, "d": 0.0e-6} - params = LoadParams() + params = load_BLS_params() # Stimulation parameters Fdrive = 350e3 # Hz Adrive = 100e3 # Pa tstim = 10e-3 # s toffset = 1e-3 # s # Initialize RS neuron rs_neuron = channels.CorticalRS() # Initialize solver solver = SolverUS(geom, params, rs_neuron, Fdrive) # Run short classic simulation of the system solver.run(rs_neuron, Fdrive, Adrive, tstim, toffset, sim_type='hybrid') # If no error is raised, test is passed successfully logger.info('Passed test: A-STIM hybrid simulation') if __name__ == '__main__': logger.info('Starting tests') test_MECH() test_resting_potential() test_ESTIM() test_ASTIM() test_ASTIM_classic() test_ASTIM_hybrid() logger.info('All tests successfully passed') diff --git a/tests/test_graphs.py b/tests/test_graphs.py index 067eb66..05f0236 100644 --- a/tests/test_graphs.py +++ b/tests/test_graphs.py @@ -1,91 +1,91 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Date: 2017-06-14 18:37:45 # @Email: theo.lemaire@epfl.ch # @Last Modified by: Theo Lemaire -# @Last Modified time: 2017-08-25 10:43:11 +# @Last Modified time: 2017-08-28 14:19:39 ''' Test the basic functionalities of the package and output graphs of the call flows. ''' import logging from pycallgraph import PyCallGraph from pycallgraph.output import GraphvizOutput import PointNICE -from PointNICE.utils import LoadParams +from PointNICE.utils import load_BLS_params from PointNICE.channels import CorticalRS # Set logging options logging.basicConfig(format='%(asctime)s %(message)s', datefmt='%d/%m/%Y %H:%M:%S:') logger = logging.getLogger('PointNICE') logger.setLevel(logging.DEBUG) # Create Graphviz output object graphviz = GraphvizOutput() # Set geometry of NBLS structure geom = {"a": 32e-9, "d": 0.0e-6} # Loading model parameters -params = LoadParams() +params = load_BLS_params() # Defining general stimulation parameters Fdrive = 3.5e5 # Hz Adrive = 1e5 # Pa PRF = 1.5e3 # Hz DF = 1 logger.info('Graph 1: BLS initialization') Cm0 = 1e-2 # membrane resting capacitance (F/m2) Qm0 = -89e-5 # membrane resting charge density (C/m2) graphviz.output_file = 'graphs/bls_init.png' with PyCallGraph(output=graphviz): bls = PointNICE.BilayerSonophore(geom, params, Fdrive, Cm0, Qm0) logger.info('Graph 2: Channels mechanism initialization') graphviz.output_file = 'graphs/sim_mech.png' with PyCallGraph(output=graphviz): bls.runMech(Fdrive, 2e4, Qm0) logger.info('Graph 3: Channels mechanism initialization') graphviz.output_file = 'graphs/channel_init.png' with PyCallGraph(output=graphviz): rs_mech = CorticalRS() logger.info('Graph 4: SolverUS initialization') graphviz.output_file = 'graphs/solver_init.png' with PyCallGraph(output=graphviz): solver = PointNICE.SolverUS(geom, params, rs_mech, Fdrive) logger.info('Graph 5: classic simulation') tstim = 1e-3 # s toffset = 1e-3 # s graphviz.output_file = 'graphs/sim_classic.png' with PyCallGraph(output=graphviz): solver.run(rs_mech, Fdrive, Adrive, tstim, toffset, PRF, DF, 'classic') logger.info('Graph 6: effective simulation') tstim = 30e-3 # s toffset = 10e-3 # s graphviz.output_file = 'graphs/sim_effective.png' with PyCallGraph(output=graphviz): solver.run(rs_mech, Fdrive, Adrive, tstim, toffset, PRF, DF, 'effective') logger.info('Graph 7: hybrid simulation') tstim = 10e-3 # s toffset = 1e-3 # s graphviz.output_file = 'graphs/sim_hybrid.png' with PyCallGraph(output=graphviz): solver.run(rs_mech, Fdrive, Adrive, tstim, toffset, PRF, DF, 'hybrid')