diff --git a/PySONIC/constants.py b/PySONIC/constants.py index a16dce7..d782e16 100644 --- a/PySONIC/constants.py +++ b/PySONIC/constants.py @@ -1,72 +1,72 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2016-11-04 13:23:31 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-08-26 14:52:37 +# @Last Modified time: 2019-08-26 15:14:47 ''' Numerical constants used in the package. ''' # -------------------------- Biophysical constants -------------------------- FARADAY = 9.64853e4 # Faraday constant (C/mol) Rg = 8.31342 # Universal gas constant (Pa.m^3.mol^-1.K^-1 or J.mol^-1.K^-1) Z_Ca = 2 # Calcium valence Z_Na = 1 # Sodium valence Z_K = 1 # Potassium valence CELSIUS_2_KELVIN = 273.15 # Celsius to Kelvin conversion constant # -------------------------- Intermolecular pressure fitting -------------------------- LJFIT_PM_MAX = 1e8 # Pm value at the deflection lower bound for LJ fitting (Pa) PNET_EQ_MAX = 1e-1 # Pnet error threshold at computed equilibrium position (Pa) -PMAVG_STD_ERR_MAX = 3000 # error threshold in intermolecular pressure nonlinear fit (Pa) +PMAVG_STD_ERR_MAX = 5e3 # error threshold in intermolecular pressure nonlinear fit (Pa) # -------------------------- Simulations -------------------------- MAX_RMSE_PTP_RATIO = 1e-4 # threshold RMSE / peak-to-peak ratio for periodic convergence Z_ERR_MAX = 1e-11 # periodic convergence threshold for deflection (m) NG_ERR_MAX = 1e-24 # periodic convergence threshold for gas content (mol) NCYCLES_MAX = 10 # max number of cycles in periodic simulations -CHARGE_RANGE = (-200e-5, 150e-5) # physiological charge range constraining the membrane (C/m2) +CHARGE_RANGE = (-300e-5, 150e-5) # physiological charge range constraining the membrane (C/m2) SOLVER_NSTEPS = 1000 # max number of steps during one ODE solver call CLASSIC_TARGET_DT = 1e-8 # target time step in output arrays of detailed simulations NPC_DENSE = 1000 # nb of samples per acoustic period in detailed simulations NPC_SPARSE = 40 # nb of samples per acoustic period in sparse simulations HYBRID_UPDATE_INTERVAL = 5e-4 # time interval between two hybrid integrations (s) DT_EFFECTIVE = 5e-5 # time step for effective integration (s) MIN_SAMPLES_PER_PULSE_INTERVAL = 1 # minimal number of time points per pulse interval (TON of TOFF) # -------------------------- Post-processing -------------------------- SPIKE_MIN_DT = 5e-4 # minimal time interval for spike detection on charge signal (s) SPIKE_MIN_QAMP = 5e-5 # threshold amplitude for spike detection on charge signal (C/m2) SPIKE_MIN_QPROM = 20e-5 # threshold prominence for spike detection on charge signal (C/m2) SPIKE_MIN_VAMP = 10.0 # threshold amplitude for spike detection on potential signal (mV) SPIKE_MIN_VPROM = 20.0 # threshold prominence for spike detection on potential signal (mV) MIN_NSPIKES_SPECTRUM = 3 # minimum number of spikes to compute firing rate spectrum # -------------------------- Titrations -------------------------- AMP_UPPER_BOUND_ESTIM = 50.0 # initial current density upper bound for titration (mA/m2) THRESHOLD_CONV_RANGE_ESTIM = 0.1 # current density search range threshold for titration (mA/m2) THRESHOLD_CONV_RANGE_ASTIM = 1e2 # acoustic pressure search range threshold for titration (Pa) # -------------------------- QSS stability analysis -------------------------- QSS_REL_OFFSET = .05 # relative state perturbation amplitude: s = s0 * (1 +/- x) QSS_HISTORY_INTERVAL = 30e-3 # recent history interval (s) QSS_INTEGRATION_INTERVAL = 1e-3 # iterative integration interval (s) QSS_MAX_INTEGRATION_DURATION = 1000e-3 # iterative integration interval (s) QSS_Q_CONV_THR = 1e-7 # max. charge deviation to infer convergence (C/m2) QSS_Q_DIV_THR = 1e-4 # min. charge deviation to infer divergence (C/m2) TMIN_STABILIZATION = 500e-3 # time window for stabilization analysis (s) def getConstantsDict(): cdict = {} for k, v in globals().items(): if not k.startswith('__') and k != 'getConstantsDict': cdict[k] = v return cdict diff --git a/PySONIC/core/bls.py b/PySONIC/core/bls.py index dc275f0..fbd8877 100644 --- a/PySONIC/core/bls.py +++ b/PySONIC/core/bls.py @@ -1,790 +1,790 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2016-09-29 16:16:19 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-08-14 17:40:22 +# @Last Modified time: 2019-08-26 15:18:05 from enum import Enum import os import json import numpy as np import pandas as pd import scipy.integrate as integrate from scipy.optimize import brentq, curve_fit from .model import Model from .simulators import PeriodicSimulator from ..utils import logger, si_format, debug from ..constants import * class PmCompMethod(Enum): ''' Enum: types of computation method for the intermolecular pressure ''' direct = 1 predict = 2 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 lookup(func): ''' Load parameters from lookup file, or compute them and sotre them in lookup file. ''' lookup_path = os.path.join(os.path.split(__file__)[0], 'bls_lookups.json') def wrapper(obj): # lookup_path = obj.getLookupsPath() akey = '{:.1f}'.format(obj.a * 1e9) Qkey = '{:.2f}'.format(obj.Qm0 * 1e5) # Open lookup files try: with open(lookup_path, 'r') as fh: lookups = json.load(fh) except FileNotFoundError: lookups = {} # If info not in lookups, compute parameters and add them if akey not in lookups or Qkey not in lookups[akey]: func(obj) if akey not in lookups: lookups[akey] = {Qkey: {'LJ_approx': obj.LJ_approx, 'Delta_eq': obj.Delta}} else: lookups[akey][Qkey] = {'LJ_approx': obj.LJ_approx, 'Delta_eq': obj.Delta} logger.debug('Saving BLS derived parameters to lookup file') with open(lookup_path, 'w') as fh: json.dump(lookups, fh, indent=2) # If lookup exists, load parameters from it else: logger.debug('Loading BLS derived parameters from lookup file') obj.LJ_approx = lookups[akey][Qkey]['LJ_approx'] obj.Delta = lookups[akey][Qkey]['Delta_eq'] return wrapper class BilayerSonophore(Model): ''' Definition of the Bilayer Sonophore Model - geometry - pressure terms - cavitation dynamics ''' # BIOMECHANICAL PARAMETERS T = 309.15 # Temperature (K) delta0 = 2.0e-9 # Thickness of the leaflet (m) Delta_ = 1.4e-9 # Initial gap between the two leaflets on a non-charged membrane at equil. (m) pDelta = 1.0e5 # Attraction/repulsion pressure coefficient (Pa) m = 5.0 # Exponent in the repulsion term (dimensionless) n = 3.3 # Exponent in the attraction term (dimensionless) rhoL = 1075.0 # Density of the surrounding fluid (kg/m^3) muL = 7.0e-4 # Dynamic viscosity of the surrounding fluid (Pa.s) muS = 0.035 # Dynamic viscosity of the leaflet (Pa.s) kA = 0.24 # Area compression modulus of the leaflet (N/m) alpha = 7.56 # Tissue shear loss modulus frequency coefficient (Pa.s) C0 = 0.62 # Initial gas molar concentration in the surrounding fluid (mol/m^3) kH = 1.613e5 # Henry's constant (Pa.m^3/mol) P0 = 1.0e5 # Static pressure in the surrounding fluid (Pa) Dgl = 3.68e-9 # Diffusion coefficient of gas in the fluid (m^2/s) xi = 0.5e-9 # Boundary layer thickness for gas transport across leaflet (m) c = 1515.0 # Speed of sound in medium (m/s) # BIOPHYSICAL PARAMETERS epsilon0 = 8.854e-12 # Vacuum permittivity (F/m) epsilonR = 1.0 # Relative permittivity of intramembrane cavity (dimensionless) tscale = 'us' # relevant temporal scale of the model simkey = 'MECH' # keyword used to characterize simulations made with this model def __init__(self, a, Cm0, Qm0, Fdrive=None, embedding_depth=0.0): ''' Constructor of the class. :param a: in-plane radius of the sonophore structure within the membrane (m) :param Cm0: membrane resting capacitance (F/m2) :param Qm0: membrane resting charge density (C/m2) :param Fdrive: frequency of acoustic perturbation (Hz) :param embedding_depth: depth of the embedding tissue around the membrane (m) ''' # Extract resting constants and geometry self.Cm0 = Cm0 self.Qm0 = Qm0 self.a = a self.d = embedding_depth self.S0 = np.pi * self.a**2 # Derive frequency-dependent tissue elastic modulus if Fdrive is not None: G_tissue = self.alpha * Fdrive # G'' (Pa) self.kA_tissue = 2 * G_tissue * self.d # kA of the tissue layer (N/m) else: self.kA_tissue = 0. # Compute Pm params self.computePMparams() # Compute initial volume and gas content self.V0 = np.pi * self.Delta * self.a**2 self.ng0 = self.gasPa2mol(self.P0, self.V0) def __repr__(self): s = '{}({:.1f} nm'.format(self.__class__.__name__, self.a * 1e9) if self.d > 0.: s += ', d={}m'.format(si_format(self.d, precision=1, space=' ')) return s + ')' @staticmethod def inputs(): return { 'a': { 'desc': 'sonophore radius', 'label': 'a', 'unit': 'nm', 'factor': 1e9, 'precision': 0 }, 'Fdrive': { 'desc': 'US frequency', 'label': 'f', 'unit': 'kHz', 'factor': 1e-3, 'precision': 0 }, 'Adrive': { 'desc': 'US pressure amplitude', 'label': 'A', 'unit': 'kPa', 'factor': 1e-3, 'precision': 2 }, 'Qm': { 'desc': 'membrane charge density', 'label': 'Q_m', 'unit': 'nC/cm^2', 'factor': 1e5, 'precision': 1 } } def filecodes(self, Fdrive, Adrive, Qm, PmCompMethod='predict'): return { 'simkey': self.simkey, 'a': '{:.0f}nm'.format(self.a * 1e9), 'Fdrive': '{:.0f}kHz'.format(Fdrive * 1e-3), 'Adrive': '{:.2f}kPa'.format(Adrive * 1e-3), 'Qm': '{:.1f}nCcm2'.format(Qm * 1e5) } @staticmethod def getPltVars(wrapleft='df["', wrapright='"]'): return { 'Pac': { 'desc': 'acoustic pressure', 'label': 'P_{AC}', 'unit': 'kPa', 'factor': 1e-3, 'func': 'Pacoustic({0}t{1}, meta["Adrive"] * {0}stimstate{1}, meta["Fdrive"])'.format( wrapleft, wrapright) }, 'Z': { 'desc': 'leaflets deflection', 'label': 'Z', 'unit': 'nm', 'factor': 1e9, 'bounds': (-1.0, 10.0) }, 'ng': { 'desc': 'gas content', 'label': 'n_g', 'unit': '10^{-22}\ mol', 'factor': 1e22, 'bounds': (1.0, 15.0) }, 'Pmavg': { 'desc': 'average intermolecular pressure', 'label': 'P_M', 'unit': 'kPa', 'factor': 1e-3, 'func': 'PMavgpred({0}Z{1})'.format(wrapleft, wrapright) }, 'Telastic': { 'desc': 'leaflet elastic tension', 'label': 'T_E', 'unit': 'mN/m', 'factor': 1e3, 'func': 'TEleaflet({0}Z{1})'.format(wrapleft, wrapright) }, 'Cm': { 'desc': 'membrane capacitance', 'label': 'C_m', 'unit': 'uF/cm^2', 'factor': 1e2, 'bounds': (0.0, 1.5), 'func': 'v_capacitance({0}Z{1})'.format(wrapleft, wrapright) } } @staticmethod def getPltScheme(): return { 'P_{AC}': ['Pac'], 'Z': ['Z'], 'n_g': ['ng'] } def curvrad(self, Z): ''' Leaflet curvature radius (signed variable) :param Z: leaflet apex deflection (m) :return: leaflet curvature radius (m) ''' if Z == 0.0: return np.inf else: return (self.a**2 + Z**2) / (2 * Z) def v_curvrad(self, Z): ''' Vectorized curvrad function ''' return np.array(list(map(self.curvrad, Z))) def surface(self, Z): ''' Surface area of the stretched leaflet (spherical cap formula) :param Z: leaflet apex deflection (m) :return: stretched leaflet surface (m^2) ''' return np.pi * (self.a**2 + Z**2) def volume(self, Z): ''' Volume of the inter-leaflet space (cylinder +/- 2 spherical caps) :param Z: leaflet apex deflection (m) :return: bilayer sonophore inner volume (m^3) ''' return np.pi * self.a**2 * self.Delta\ * (1 + (Z / (3 * self.Delta) * (3 + Z**2 / self.a**2))) def arealStrain(self, Z): ''' Areal strain of the stretched leaflet epsilon = (S - S0)/S0 = (Z/a)^2 :param Z: leaflet apex deflection (m) :return: areal strain (dimensionless) ''' return (Z / self.a)**2 def capacitance(self, Z): ''' Membrane capacitance (parallel-plate capacitor evaluated at average inter-layer distance) :param Z: leaflet apex deflection (m) :return: capacitance per unit area (F/m2) ''' if Z == 0.0: return self.Cm0 else: return ((self.Cm0 * self.Delta / self.a**2) * (Z + (self.a**2 - Z**2 - Z * self.Delta) / (2 * Z) * np.log((2 * Z + self.Delta) / self.Delta))) def v_capacitance(self, Z): ''' Vectorized capacitance function ''' return np.array(list(map(self.capacitance, Z))) def derCapacitance(self, Z, U): ''' Evolution of membrane capacitance :param Z: leaflet apex deflection (m) :param U: leaflet apex deflection velocity (m/s) :return: time derivative of capacitance per unit area (F/m2.s) ''' dCmdZ = ((self.Cm0 * self.Delta / self.a**2) * ((Z**2 + self.a**2) / (Z * (2 * Z + self.Delta)) - ((Z**2 + self.a**2) * np.log((2 * Z + self.Delta) / self.Delta)) / (2 * Z**2))) return dCmdZ * U @staticmethod def localDeflection(r, Z, R): ''' Local leaflet deflection at specific radial distance (signed) :param r: in-plane distance from center of the sonophore (m) :param Z: leaflet apex deflection (m) :param R: leaflet curvature radius (m) :return: local transverse leaflet deviation (m) ''' if np.abs(Z) == 0.0: return 0.0 else: return np.sign(Z) * (np.sqrt(R**2 - r**2) - np.abs(R) + np.abs(Z)) @staticmethod def Pacoustic(t, Adrive, Fdrive, phi=np.pi): ''' Time-varying acoustic pressure :param t: time (s) :param Adrive: acoustic drive amplitude (Pa) :param Fdrive: acoustic drive frequency (Hz) :param phi: acoustic drive phase (rad) ''' return Adrive * np.sin(2 * np.pi * Fdrive * t - phi) def PMlocal(self, r, Z, R): ''' Local intermolecular pressure :param r: in-plane distance from center of the sonophore (m) :param Z: leaflet apex deflection (m) :param R: leaflet curvature radius (m) :return: local intermolecular pressure (Pa) ''' z = self.localDeflection(r, Z, R) relgap = (2 * z + self.Delta) / self.Delta_ return self.pDelta * ((1 / relgap)**self.m - (1 / relgap)**self.n) def PMavg(self, Z, R, S): ''' Average intermolecular pressure across the leaflet (computed by quadratic integration) :param Z: leaflet apex outward deflection value (m) :param R: leaflet curvature radius (m) :param S: surface of the stretched leaflet (m^2) :return: averaged intermolecular resultant pressure (Pa) .. warning:: quadratic integration is computationally expensive. ''' # Integrate intermolecular force over an infinitely thin ring of radius r from 0 to a fTotal, _ = integrate.quad(lambda r, Z, R: 2 * np.pi * r * self.PMlocal(r, Z, R), 0, self.a, args=(Z, R)) return fTotal / S def v_PMavg(self, Z, R, S): ''' Vectorized PMavg function ''' return np.array(list(map(self.PMavg, Z, R, S))) def LJfitPMavg(self): ''' Determine optimal parameters of a Lennard-Jones expression approximating the average intermolecular pressure. These parameters are obtained by a nonlinear fit of the Lennard-Jones function for a range of deflection values between predetermined Zmin and Zmax. :return: 3-tuple with optimized LJ parameters for PmAvg prediction (Map) and the standard and max errors of the prediction in the fitting range (in Pascals) ''' # Determine lower bound of deflection range: when Pm = Pmmax PMmax = LJFIT_PM_MAX # Pa Zminlb = -0.49 * self.Delta Zminub = 0.0 Zmin = brentq(lambda Z, Pmmax: self.PMavg(Z, self.curvrad(Z), self.surface(Z)) - PMmax, Zminlb, Zminub, args=(PMmax), xtol=1e-16) # Create vectors for geometric variables Zmax = 2 * self.a Z = np.arange(Zmin, Zmax, 1e-11) Pmavg = self.v_PMavg(Z, self.v_curvrad(Z), self.surface(Z)) # Compute optimal nonlinear fit of custom LJ function with initial guess x0_guess = self.delta0 C_guess = 0.1 * self.pDelta nrep_guess = self.m nattr_guess = self.n pguess = (x0_guess, C_guess, nrep_guess, nattr_guess) popt, _ = curve_fit(lambda x, x0, C, nrep, nattr: LennardJones(x, self.Delta, x0, C, nrep, nattr), - Z, Pmavg, p0=pguess, maxfev=10000) + Z, Pmavg, p0=pguess, maxfev=100000) (x0_opt, C_opt, nrep_opt, nattr_opt) = popt Pmavg_fit = LennardJones(Z, self.Delta, x0_opt, C_opt, nrep_opt, nattr_opt) # Compute prediction error residuals = Pmavg - Pmavg_fit ss_res = np.sum(residuals**2) N = residuals.size std_err = np.sqrt(ss_res / N) max_err = max(np.abs(residuals)) logger.debug('LJ approx: x0 = %.2f nm, C = %.2f kPa, m = %.2f, n = %.2f', x0_opt * 1e9, C_opt * 1e-3, nrep_opt, nattr_opt) LJ_approx = {"x0": x0_opt, "C": C_opt, "nrep": nrep_opt, "nattr": nattr_opt} return (LJ_approx, std_err, max_err) @lookup def computePMparams(self): # Find Delta that cancels out Pm + Pec at Z = 0 (m) if self.Qm0 == 0.0: D_eq = self.Delta_ else: (D_eq, Pnet_eq) = self.findDeltaEq(self.Qm0) assert Pnet_eq < PNET_EQ_MAX, 'High Pnet at Z = 0 with ∆ = %.2f nm' % (D_eq * 1e9) self.Delta = D_eq # Find optimal Lennard-Jones parameters to approximate PMavg (self.LJ_approx, std_err, _) = self.LJfitPMavg() assert std_err < PMAVG_STD_ERR_MAX, 'High error in PmAvg nonlinear fit:'\ ' std_err = %.2f Pa' % std_err def PMavgpred(self, Z): ''' Approximated average intermolecular pressure (using nonlinearly fitted Lennard-Jones function) :param Z: leaflet apex deflection (m) :return: predicted average intermolecular pressure (Pa) ''' return LennardJones(Z, self.Delta, self.LJ_approx['x0'], self.LJ_approx['C'], self.LJ_approx['nrep'], self.LJ_approx['nattr']) def Pelec(self, Z, Qm): ''' Electrical pressure term :param Z: leaflet apex deflection (m) :param Qm: membrane charge density (C/m2) :return: electrical pressure (Pa) ''' relS = self.S0 / self.surface(Z) abs_perm = self.epsilon0 * self.epsilonR # F/m return - relS * Qm**2 / (2 * abs_perm) # Pa def findDeltaEq(self, Qm): ''' Compute the Delta that cancels out the (Pm + Pec) equation at Z = 0 for a given membrane charge density, using the Brent method to refine the pressure root iteratively. :param Qm: membrane charge density (C/m2) :return: equilibrium value (m) and associated pressure (Pa) ''' def dualPressure(Delta): x = (self.Delta_ / Delta) return (self.pDelta * (x**self.m - x**self.n) + self.Pelec(0.0, Qm)) Delta_eq = brentq(dualPressure, 0.1 * self.Delta_, 2.0 * self.Delta_, xtol=1e-16) logger.debug('∆eq = %.2f nm', Delta_eq * 1e9) return (Delta_eq, dualPressure(Delta_eq)) def gasFlux(self, Z, P): ''' Gas molar flux through the sonophore boundary layers :param Z: leaflet apex deflection (m) :param P: internal gas pressure (Pa) :return: gas molar flux (mol/s) ''' dC = self.C0 - P / self.kH return 2 * self.surface(Z) * self.Dgl * dC / self.xi @classmethod def gasmol2Pa(cls, ng, V): ''' Internal gas pressure for a given molar content :param ng: internal molar content (mol) :param V: sonophore inner volume (m^3) :return: internal gas pressure (Pa) ''' return ng * Rg * cls.T / V @classmethod def gasPa2mol(cls, P, V): ''' Internal gas molar content for a given pressure :param P: internal gas pressure (Pa) :param V: sonophore inner volume (m^3) :return: internal gas molar content (mol) ''' return P * V / (Rg * cls.T) def PtotQS(self, Z, ng, Qm, Pac, Pm_comp_method): ''' Net quasi-steady pressure for a given acoustic pressure (Ptot = Pm + Pg + Pec - P0 - Pac) :param Z: leaflet apex deflection (m) :param ng: internal molar content (mol) :param Qm: membrane charge density (C/m2) :param Pac: acoustic pressure (Pa) :param Pm_comp_method: computation method for average intermolecular pressure :return: total balance pressure (Pa) ''' if Pm_comp_method is PmCompMethod.direct: Pm = self.PMavg(Z, self.curvrad(Z), self.surface(Z)) elif Pm_comp_method is PmCompMethod.predict: Pm = self.PMavgpred(Z) return Pm + self.gasmol2Pa(ng, self.volume(Z)) - self.P0 - Pac + self.Pelec(Z, Qm) def balancedefQS(self, ng, Qm, Pac=0.0, Pm_comp_method=PmCompMethod.predict): ''' Quasi-steady equilibrium deflection for a given acoustic pressure (computed by approximating the root of quasi-steady pressure) :param ng: internal molar content (mol) :param Qm: membrane charge density (C/m2) :param Pac: external acoustic perturbation (Pa) :param Pm_comp_method: computation method for average intermolecular pressure :return: leaflet deflection canceling quasi-steady pressure (m) ''' lb, ub = -0.49 * self.Delta, self.a Plb, Pub = [self.PtotQS(x, ng, Qm, Pac, Pm_comp_method) for x in [lb, ub]] assert (Plb > 0 > Pub), '[{}, {}] is not a sign changing interval for PtotQS'.format(lb, ub) return brentq(self.PtotQS, lb, ub, args=(ng, Qm, Pac, Pm_comp_method), xtol=1e-16) def TEleaflet(self, Z): ''' Elastic tension in leaflet :param Z: leaflet apex deflection (m) :return: circumferential elastic tension (N/m) ''' return self.kA * self.arealStrain(Z) def TEtissue(self, Z): ''' Elastic tension in surrounding viscoelastic layer :param Z: leaflet apex deflection (m) :return: circumferential elastic tension (N/m) ''' return self.kA_tissue * self.arealStrain(Z) def TEtot(self, Z): ''' Total elastic tension (leaflet + surrounding viscoelastic layer) :param Z: leaflet apex deflection (m) :return: circumferential elastic tension (N/m) ''' return self.TEleaflet(Z) + self.TEtissue(Z) def PEtot(self, Z, R): ''' Total elastic tension pressure (leaflet + surrounding viscoelastic layer) :param Z: leaflet apex deflection (m) :param R: leaflet curvature radius (m) :return: elastic tension pressure (Pa) ''' return - self.TEtot(Z) / R @classmethod def PVleaflet(cls, U, R): ''' Viscous stress pressure in leaflet :param U: leaflet apex deflection velocity (m/s) :param R: leaflet curvature radius (m) :return: leaflet viscous stress pressure (Pa) ''' return - 12 * U * cls.delta0 * cls.muS / R**2 @classmethod def PVfluid(cls, U, R): ''' Viscous stress pressure in surrounding medium :param U: leaflet apex deflection velocity (m/s) :param R: leaflet curvature radius (m) :return: fluid viscous stress pressure (Pa) ''' return - 4 * U * cls.muL / np.abs(R) @classmethod def accP(cls, Ptot, R): ''' Leaflet transverse acceleration resulting from pressure imbalance :param Ptot: net pressure (Pa) :param R: leaflet curvature radius (m) :return: pressure-driven acceleration (m/s^2) ''' return Ptot / (cls.rhoL * np.abs(R)) @staticmethod def accNL(U, R): ''' Leaflet transverse nonlinear acceleration :param U: leaflet apex deflection velocity (m/s) :param R: leaflet curvature radius (m) :return: nonlinear acceleration term (m/s^2) .. note:: A simplified version of nonlinear acceleration (neglecting dR/dH) is used here. ''' # return - (3/2 - 2*R/H) * U**2 / R return -(3 * U**2) / (2 * R) @staticmethod def checkInputs(Fdrive, Adrive, Qm, phi): ''' Check validity of stimulation parameters :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param phi: acoustic drive phase (rad) :param Qm: imposed membrane charge density (C/m2) ''' if not all(isinstance(param, float) for param in [Fdrive, Adrive, Qm, phi]): raise TypeError('Invalid stimulation parameters (must be float typed)') if Fdrive <= 0: raise ValueError('Invalid US driving frequency: {} kHz (must be strictly positive)' .format(Fdrive * 1e-3)) if Adrive < 0: raise ValueError('Invalid US pressure amplitude: {} kPa (must be positive or null)' .format(Adrive * 1e-3)) if Qm < CHARGE_RANGE[0] or Qm > CHARGE_RANGE[1]: raise ValueError('Invalid applied charge: {} nC/cm2 (must be within [{}, {}] interval' .format(Qm * 1e5, CHARGE_RANGE[0] * 1e5, CHARGE_RANGE[1] * 1e5)) if phi < 0 or phi >= 2 * np.pi: raise ValueError('Invalid US pressure phase: {:.2f} rad (must be within [0, 2 PI[ rad' .format(phi)) def derivatives(self, t, y, Fdrive, Adrive, Qm, phi, Pm_comp_method=PmCompMethod.predict): ''' Evolution of the mechanical system :param t: time instant (s) :param y: vector of HH system variables at time t :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param Qm: membrane charge density (F/m2) :param phi: acoustic drive phase (rad) :param Pm_comp_method: computation method for average intermolecular pressure :return: vector of mechanical system derivatives at time t ''' # Split input vector explicitly U, Z, ng = y # Correct deflection value is below critical compression if Z < -0.5 * self.Delta: logger.warning('Deflection out of range: Z = %.2f nm', Z * 1e9) Z = -0.49 * self.Delta # Compute curvature radius R = self.curvrad(Z) # Compute total pressure Pg = self.gasmol2Pa(ng, self.volume(Z)) if Pm_comp_method is PmCompMethod.direct: Pm = self.PMavg(Z, self.curvrad(Z), self.surface(Z)) elif Pm_comp_method is PmCompMethod.predict: Pm = self.PMavgpred(Z) Ptot = (Pm + Pg - self.P0 - self.Pacoustic(t, Adrive, Fdrive, phi) + self.PEtot(Z, R) + self.PVleaflet(U, R) + self.PVfluid(U, R) + self.Pelec(Z, Qm)) # Compute derivatives dUdt = self.accP(Ptot, R) + self.accNL(U, R) dZdt = U dngdt = self.gasFlux(Z, Pg) # Return derivatives vector return [dUdt, dZdt, dngdt] def computeInitialDeflection(self, Adrive, Fdrive, phi, Qm, dt, Pm_comp_method=PmCompMethod.predict): ''' Compute non-zero deflection value for a small perturbation (solving quasi-steady equation). ''' Pac = self.Pacoustic(dt, Adrive, Fdrive, phi) return self.balancedefQS(self.ng0, Qm, Pac, Pm_comp_method) @Model.addMeta def simulate(self, Fdrive, Adrive, Qm, phi=np.pi, Pm_comp_method=PmCompMethod.predict): ''' Simulate until periodic stabilization for a specific set of ultrasound parameters, and return output data in a dataframe. :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param phi: acoustic drive phase (rad) :param Qm: imposed membrane charge density (C/m2) :param Pm_comp_method: type of method used to compute average intermolecular pressure :return: 2-tuple with the output dataframe and computation time. ''' logger.info('%s: simulation @ f = %sHz, A = %sPa, Q = %sC/cm2', self, *si_format([Fdrive, Adrive, Qm * 1e-4], 2, space=' ')) # Check validity of stimulation parameters self.checkInputs(Fdrive, Adrive, Qm, phi) # Determine time step dt = 1 / (NPC_DENSE * Fdrive) # Compute initial non-zero deflection Z = self.computeInitialDeflection(Adrive, Fdrive, phi, Qm, dt, Pm_comp_method=Pm_comp_method) # Set initial conditions y0 = np.array([0., 0., self.ng0]) y1 = np.array([0., Z, self.ng0]) # Initialize simulator and compute solution simulator = PeriodicSimulator( lambda t, y: self.derivatives(t, y, Fdrive, Adrive, Qm, phi, Pm_comp_method), ivars_to_check=[1, 2]) t, y, stim = simulator(y1, dt, 1. / Fdrive) # Prepend initial conditions (prior to stimulation) t, y, stim = simulator.prependSolution(t, y, stim, y0=y0) # Set last stimulation state to zero stim[-1] = 0 # Store output in dataframe and return return pd.DataFrame({ 't': t, 'stimstate': stim, 'Z': y[:, 1], 'ng': y[:, 2] }) def meta(self, Fdrive, Adrive, Qm, Pm_comp_method): return { 'simkey': self.simkey, 'a': self.a, 'd': self.d, 'Cm0': self.Cm0, 'Qm0': self.Qm0, 'Fdrive': Fdrive, 'Adrive': Adrive, 'Qm': Qm, 'Pm_comp_method': Pm_comp_method } def getCycleProfiles(self, Fdrive, Adrive, Qm): ''' Simulate mechanical system and compute pressures over the last acoustic cycle :param Fdrive: acoustic drive frequency (Hz) :param Adrive: acoustic drive amplitude (Pa) :param Qm: imposed membrane charge density (C/m2) :return: dataframe with the time, kinematic and pressure profiles over the last cycle. ''' # Run default simulation and retrieve last cycle solution logger.info('Running mechanical simulation (a = %sm, f = %sHz, A = %sPa)', si_format(self.a, 1), si_format(Fdrive, 1), si_format(Adrive, 1)) data = self.simulate( Fdrive, Adrive, Qm, Pm_comp_method=PmCompMethod.direct)[0].iloc[-NPC_DENSE:, :] # Extract relevant variables and de-offset time vector t, Z, ng = [data[key].values for key in ['t', 'Z', 'ng']] dt = (t[-1] - t[0]) / (NPC_DENSE - 1) t -= t[0] # Compute pressure cyclic profiles logger.info('Computing pressure cyclic profiles') R = self.v_curvrad(Z) U = np.diff(Z) / dt U = np.hstack((U, U[-1])) data = { 't': t, 'Z': Z, 'Cm': self.v_capacitance(Z), 'P_M': self.v_PMavg(Z, R, self.surface(Z)), 'P_Q': self.Pelec(Z, Qm), 'P_{VE}': self.PEtot(Z, R) + self.PVleaflet(U, R), 'P_V': self.PVfluid(U, R), 'P_G': self.gasmol2Pa(ng, self.volume(Z)), 'P_0': - np.ones(Z.size) * self.P0 } return pd.DataFrame(data, columns=data.keys()) diff --git a/PySONIC/core/bls_lookups.json b/PySONIC/core/bls_lookups.json index 20b442f..498bc37 100644 --- a/PySONIC/core/bls_lookups.json +++ b/PySONIC/core/bls_lookups.json @@ -1,658 +1,685 @@ { "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 }, "-58.00": { "LJ_approx": { "x0": 1.6131662659976035e-09, "C": 20156.47605325608, "nrep": 3.9204295233162925, "nattr": 1.0393049952285334 }, "Delta_eq": 1.2922508866011204e-09 }, "-140.00": { "LJ_approx": { "x0": 3.5396589580589484e-09, "C": 1255.8321160636908, "nrep": 3.879907809497444, "nattr": 0.4190657482583384 }, "Delta_eq": 1.1019265101358682e-09 + }, + "-200.00": { + "LJ_approx": { + "x0": 5.467498689940948e-09, + "C": 298.4575230665454, + "nrep": 3.8382806376855165, + "nattr": 0.0014805372073950717 + }, + "Delta_eq": 1.0050623818936587e-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 }, "-58.00": { "LJ_approx": { "x0": 1.6095250707781465e-09, "C": 20329.27848623273, "nrep": 3.92057191136993, "nattr": 1.0267862446034561 }, "Delta_eq": 1.2922508866011204e-09 }, "-140.00": { "LJ_approx": { "x0": 3.5398097121754107e-09, "C": 1255.8690004347025, "nrep": 3.8798490490747404, "nattr": 0.4604604439108636 }, "Delta_eq": 1.1019265101358682e-09 + }, + "-200.00": { + "LJ_approx": { + "x0": 6.6005456899992844e-09, + "C": 144.85407475420143, + "nrep": 3.838295127819493, + "nattr": 0.00992375961984237 + }, + "Delta_eq": 1.0050623818936587e-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 }, "-71.90": { "LJ_approx": { "x0": 1.7286986021591825e-09, "C": 16087.514816365254, "nrep": 3.9147885683678543, "nattr": 1.012616990226475 }, "Delta_eq": 1.2553492695740507e-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 }, "-58.00": { "LJ_approx": { "x0": 1.6196423302303851e-09, "C": 19847.346684716063, "nrep": 3.920294617611314, "nattr": 1.0573210567487794 }, "Delta_eq": 1.2922508866011204e-09 }, "-140.00": { "LJ_approx": { "x0": 3.536300349187651e-09, "C": 1259.968827699725, "nrep": 3.8800144975132485, "nattr": 0.35722933384459793 }, "Delta_eq": 1.1019265101358682e-09 + }, + "-200.00": { + "LJ_approx": { + "x0": 4.448077576688761e-09, + "C": 658.9502615105886, + "nrep": 3.8382596040618573, + "nattr": -0.0007412836249291593 + }, + "Delta_eq": 1.0050623818936587e-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 } }, "30.0": { "-71.90": { "LJ_approx": { "x0": 1.7147994934556977e-09, "C": 16608.65261608246, "nrep": 3.914764637335829, "nattr": 0.9867268079339511 }, "Delta_eq": 1.2553492695740507e-09 } }, "25.4": { "-71.90": { "LJ_approx": { "x0": 1.7162265501918752e-09, "C": 16555.217050637588, "nrep": 3.9147464050871856, "nattr": 0.9900398844251959 }, "Delta_eq": 1.2553492695740507e-09 }, "-61.93": { "LJ_approx": { "x0": 1.6408385667642563e-09, "C": 19099.831853142834, "nrep": 3.9188405082474103, "nattr": 1.0301853851264013 }, "Delta_eq": 1.281743450351987e-09 } }, "40.3": { "-71.90": { "LJ_approx": { "x0": 1.7127058661258177e-09, "C": 16686.9995037205, "nrep": 3.914800799246736, "nattr": 0.981479342025535 }, "Delta_eq": 1.2553492695740507e-09 }, "-61.93": { "LJ_approx": { "x0": 1.637531858311385e-09, "C": 19247.790954282773, "nrep": 3.918928365198691, "nattr": 1.020450016688262 }, "Delta_eq": 1.281743450351987e-09 } }, "45.0": { "-71.90": { "LJ_approx": { "x0": 1.712051746586677e-09, "C": 16711.456749369456, "nrep": 3.914814642254888, "nattr": 0.979726839958776 }, "Delta_eq": 1.2553492695740507e-09 } }, "20.0": { "-71.90": { "LJ_approx": { "x0": 1.7186388439609698e-09, "C": 16464.85454836168, "nrep": 3.9147266088229755, "nattr": 0.9952322612710219 }, "Delta_eq": 1.2553492695740507e-09 } }, "60.0": { "-71.90": { "LJ_approx": { "x0": 1.710603828012074e-09, "C": 16765.5278159216, "nrep": 3.9148503824940146, "nattr": 0.9756024890579372 }, "Delta_eq": 1.2553492695740507e-09 } }, "34.0": { "-71.90": { "LJ_approx": { "x0": 1.7138494069983225e-09, "C": 16644.21766668786, "nrep": 3.9147795488410644, "nattr": 0.9844103642751335 }, "Delta_eq": 1.2553492695740507e-09 } }, "22.6": { "-71.90": { "LJ_approx": { "x0": 1.717347083819261e-09, "C": 16513.244775760173, "nrep": 3.914735582752492, "nattr": 0.9925083846044295 }, "Delta_eq": 1.2553492695740507e-09 } }, "45.3": { "-71.90": { "LJ_approx": { "x0": 1.7120143529753677e-09, "C": 16712.8535932463, "nrep": 3.9148154954695285, "nattr": 0.9796234512533043 }, "Delta_eq": 1.2553492695740507e-09 } }, "31.0": { "0.00": { "LJ_approx": { "x0": 1.429704934686545e-09, "C": 28734.55271989346, "nrep": 3.9338783401809607, "nattr": 1.155904475115303 }, "Delta_eq": 1.4e-09 }, "-71.90": { "LJ_approx": { "x0": 1.714533452342314e-09, "C": 16618.61152256, "nrep": 3.9147686086169227, "nattr": 0.9860860917336786 }, "Delta_eq": 1.2553492695740507e-09 } }, "31.1": { "-71.90": { "LJ_approx": { "x0": 1.714515996985078e-09, "C": 16619.26565583038, "nrep": 3.914768856863524, "nattr": 0.986044873313721 }, "Delta_eq": 1.2553492695740507e-09 } } } \ No newline at end of file diff --git a/PySONIC/core/model.py b/PySONIC/core/model.py index 14c3e9d..b303216 100644 --- a/PySONIC/core/model.py +++ b/PySONIC/core/model.py @@ -1,244 +1,244 @@ # -*- coding: utf-8 -*- # @Author: Theo Lemaire # @Email: theo.lemaire@epfl.ch # @Date: 2017-08-03 11:53:04 # @Last Modified by: Theo Lemaire -# @Last Modified time: 2019-08-23 19:52:25 +# @Last Modified time: 2019-08-26 15:11:35 import os from functools import wraps from inspect import signature, getdoc import pickle import abc import inspect import numpy as np from .batches import Batch from ..utils import logger, loadData, timer, si_format, plural, debug class Model(metaclass=abc.ABCMeta): ''' Generic model interface. ''' @property @abc.abstractmethod def tscale(self): ''' Relevant temporal scale of the model. ''' raise NotImplementedError @property @abc.abstractmethod def simkey(self): ''' Keyword used to characterize simulations made with the model. ''' raise NotImplementedError @property @abc.abstractmethod def __repr__(self): ''' String representation. ''' raise NotImplementedError def params(self): ''' Return a dictionary of all model parameters (class and instance attributes) ''' def toAvoid(p): return (p.startswith('__') and p.endswith('__')) or p.startswith('_abc_') class_attrs = inspect.getmembers(self.__class__, lambda a: not(inspect.isroutine(a))) inst_attrs = inspect.getmembers(self, lambda a: not(inspect.isroutine(a))) class_attrs = [a for a in class_attrs if not toAvoid(a[0])] inst_attrs = [a for a in inst_attrs if not toAvoid(a[0]) and a not in class_attrs] params_dict = {a[0]: a[1] for a in class_attrs + inst_attrs} return params_dict @classmethod def description(cls): return getdoc(cls).split('\n', 1)[0].strip() @staticmethod @abc.abstractmethod def inputs(): ''' Return an informative dictionary on input variables used to simulate the model. ''' raise NotImplementedError @property @abc.abstractmethod def filecodes(self, *args): ''' Return a dictionary of string-encoded inputs used for file naming. ''' raise NotImplementedError def filecode(self, *args): ''' Generate file code given a specific combination of model input parameters. ''' # If meta dictionary was passed, generate inputs list from it if len(args) == 1 and isinstance(args[0], dict): meta = args[0] meta.pop('tcomp', None) nparams = len(signature(self.meta).parameters) args = list(meta.values())[-nparams:] # Create file code by joining string-encoded inputs with underscores codes = self.filecodes(*args).values() return '_'.join([x for x in codes if x is not None]) @classmethod @abc.abstractmethod def getPltVars(self, *args, **kwargs): ''' Return a dictionary with information about all plot variables related to the model. ''' raise NotImplementedError @classmethod @abc.abstractmethod def getPltScheme(self): ''' Return a dictionary model plot variables grouped by context. ''' raise NotImplementedError @staticmethod def checkOutputDir(queue, outputdir): ''' Check if an outputdir is provided in input arguments, and if so, add it as the first element of each item in the returned queue. ''' if outputdir is not None: for item in queue: item.insert(0, outputdir) else: if len(queue) > 5: logger.warning('Running more than 5 simulations without file saving') return queue @classmethod def simQueue(cls, *args, outputdir=None): return cls.checkOutputDir(Batch.createQueue(*args), outputdir) @staticmethod @abc.abstractmethod def checkInputs(self, *args): ''' Check the validity of simulation input parameters. ''' raise NotImplementedError @property @abc.abstractmethod def derivatives(self, *args, **kwargs): ''' Compute ODE derivatives for a specific set of ODE states and external parameters. ''' raise NotImplementedError @property @abc.abstractmethod def simulate(self, *args, **kwargs): ''' Simulate the model's differential system for specific input parameters and return output data in a dataframe. ''' raise NotImplementedError @classmethod @abc.abstractmethod def meta(self, *args): ''' Return an informative dictionary about model and simulation parameters. ''' raise NotImplementedError @staticmethod def addMeta(simfunc): ''' Add an informative dictionary about model and simulation parameters to simulation output ''' @wraps(simfunc) def wrapper(self, *args, **kwargs): data, tcomp = timer(simfunc)(self, *args, **kwargs) logger.debug('completed in %ss', si_format(tcomp, 1)) # Add keyword arguments from simfunc signature if not provided bound_args = signature(simfunc).bind(self, *args, **kwargs) bound_args.apply_defaults() target_args = dict(bound_args.arguments) # Try to retrieve meta information try: meta_params_names = list(signature(self.meta).parameters.keys()) meta_params = [target_args[k] for k in meta_params_names] meta = self.meta(*meta_params) except KeyError as err: - logger.error(f'Could not find {err} parameter in simulate function') + logger.error(f'Could not find {err} parameter in "{simfunc.__name__}" function') meta = {} # Add computation time to it meta['tcomp'] = tcomp # Return data with meta dict return data, meta return wrapper @staticmethod def logNSpikes(simfunc): ''' Log number of detected spikes on charge profile of simulation output. ''' @wraps(simfunc) def wrapper(self, *args, **kwargs): out = simfunc(self, *args, **kwargs) if out is None: return None data, meta = out nspikes = self.getNSpikes(data) logger.debug('{} spike{} detected'.format(nspikes, plural(nspikes))) return data, meta return wrapper @staticmethod def checkTitrate(argname): ''' If no None provided in the list of input parameters, perform a titration to find the threshold parameter and add it to the list. ''' def wrapper_with_args(simfunc): @wraps(simfunc) def wrapper(self, *args, **kwargs): # Get argument index from function signature func_args = list(signature(simfunc).parameters.keys())[1:] iarg = func_args.index(argname) # If argument is None if args[iarg] is None: # Generate new args list without argument args = list(args) new_args = args.copy() del new_args[iarg] # Perform titration to find threshold argument value xthr = self.titrate(*new_args) if np.isnan(xthr): logger.error(f'Could not find threshold {argname}') return None # Re-insert it into arguments list args[iarg] = xthr # Execute simulation function return simfunc(self, *args, **kwargs) return wrapper return wrapper_with_args def simAndSave(self, outdir, *args): ''' Simulate the model and save the results in a specific output directory. ''' out = self.simulate(*args) if out is None: return None data, meta = out if None in args: args = list(args) iNone = next(i for i, arg in enumerate(args) if arg is None) args[iNone] = meta['Adrive'] fpath = '{}/{}.pkl'.format(outdir, self.filecode(*args)) with open(fpath, 'wb') as fh: pickle.dump({'meta': meta, 'data': data}, fh) logger.debug('simulation data exported to "%s"', fpath) return fpath def getOutput(self, outdir, *args): ''' Get simulation output data for a specific parameters combination, by looking for an output file into a specific directory. If a corresponding output file is not found in the specified directory, the model is first run and results are saved in the output file. ''' fpath = '{}/{}.pkl'.format(outdir, self.filecode(*args)) if not os.path.isfile(fpath): self.simAndSave(outdir, *args) return loadData(fpath)