diff --git a/.drone.yml b/.drone.yml index 07596f0d..24b26ac4 100644 --- a/.drone.yml +++ b/.drone.yml @@ -1,17 +1,18 @@ build: - image: bjodah/bjodahimg16dev:latest + image: bjodah/bjodahimg16dev:v1.3 environment: - - CPLUS_INCLUDE_PATH=/opt/boost_1_63_0/include - - LIBRARY_PATH=/opt/boost_1_63_0/lib - - LD_LIBRARY_PATH=/opt/boost_1_63_0/lib + - CPLUS_INCLUDE_PATH=/opt/boost_1_65_0/include + - LIBRARY_PATH=/opt/boost_1_65_0/lib + - LD_LIBRARY_PATH=/opt/boost_1_65_0/lib commands: - - bash -c '[[ $(python setup.py --version) =~ ^[0-9]+.* ]]' - - export CPLUS_INCLUDE_PATH=/opt/boost_1_63_0 + - python2 -m pip install --ignore-installed --no-deps quantities==0.12.1 # https://github.com/python-quantities/python-quantities/issues/122 + - python3 -m pip install --ignore-installed --no-deps quantities==0.12.1 - PYCVODES_LAPACK=lapack,blas ./scripts/ci.sh chempy - (cd examples/; python -m pip install --force bokeh==0.11.1; for f in bokeh_*.py; do python -m bokeh html $f; done) - ./scripts/prepare_deploy.sh - PATH=/opt/miniconda3/bin:$PATH conda config --add channels bjodah # sym, pyodesys, pyneqsys - PATH=/opt/miniconda3/bin:$PATH conda build conda-recipe + - bash -c '[[ $(python setup.py --version) =~ ^[0-9]+.* ]]' - if grep "DO-NOT-MERGE!" -R . --exclude ".drone.yml"; then exit 1; fi deploy: diff --git a/chempy/__init__.py b/chempy/__init__.py index b325998b..89392aaa 100644 --- a/chempy/__init__.py +++ b/chempy/__init__.py @@ -11,5 +11,4 @@ from .reactionsystem import ReactionSystem from .henry import Henry from .util.periodic import atomic_number - -from .kinetics import * +from .kinetics import EyringParam, EyringHS, MassAction diff --git a/chempy/chemistry.py b/chempy/chemistry.py index d247af94..9706a7e8 100644 --- a/chempy/chemistry.py +++ b/chempy/chemistry.py @@ -1128,7 +1128,7 @@ def balance_stoichiometry(reactants, products, substances=None, balanced products : dict """ - from sympy import Matrix + from sympy import Matrix, gcd _intersect = set.intersection(*map(set, (reactants, products))) if _intersect: @@ -1142,7 +1142,6 @@ def balance_stoichiometry(reactants, products, substances=None, subst_keys = list(substances.keys()) cks = Substance.composition_keys(substances.values()) - nsubs = len(substances) # ?C2H2 + ?O2 -> ?CO + ?H2O # Ax = 0 @@ -1157,61 +1156,10 @@ def _get(sk, ck): return substances[sk].composition.get(ck, 0) * (-1 if sk in reactants else 1) A = Matrix([[_get(sk, ck) for sk in subst_keys] for ck in cks]) - - # A2 x = b - # - # A2: x: b: - # - # O2 CO H2O C2H2 - # C 0 1 0 x0 2 - # H 0 0 2 x1 2 - # O -2 1 1 x2 0 - - A_aug, pivot = A.rref() - if len(pivot) < nsubs-1: - raise ValueError("Unsatisfiable system of equations") - x_aug = Matrix(A_aug[:len(pivot), 1:]).LUsolve(Matrix(-A_aug[:len(pivot), 0])) - - # Reorder to original indices - x = [1] - for si in range(1, nsubs): - ai = si - 1 # augmented index - if ai in pivot: - x.append(x_aug[pivot.index(ai)]) - else: - x.append(None) - - # Now solve for the redundant x:s - for si in range(1, nsubs): - elem = x[si] - if elem is None: - # solve - col = A[:, si] - for ri, cell in enumerate(col): - if cell == 0: - continue - others = 0 - for ci, comp in enumerate(A[ri, :]): - if ci == si: - continue - if x[ci] is None: - raise NotImplementedError("Need a second LU solve") - others += comp*x[ci] - x[si] = -others/cell - break - - x = Matrix(x) - while True: - for idx in range(nsubs): - elem = x[idx] - if not elem.is_integer: - numer, denom = elem.as_numer_denom() - x *= denom - break - else: - break - if 0 in x: - raise ValueError("Unable to balance stoichiometry (did you forget a product?)") + nullsp, = A.nullspace() + if 0 in nullsp: + raise ValueError("Superfluous species given.") + x = nullsp/reduce(gcd, nullsp.T) def _x(k): return x[subst_keys.index(k)] diff --git a/chempy/kinetics/__init__.py b/chempy/kinetics/__init__.py index fe32fba1..fac5e66a 100644 --- a/chempy/kinetics/__init__.py +++ b/chempy/kinetics/__init__.py @@ -3,4 +3,6 @@ """ from __future__ import (absolute_import, division, print_function) -from .rates import MassAction +from .rates import MassAction, EyringHS +from .eyring import EyringParam +from .arrhenius import ArrheniusParam diff --git a/chempy/kinetics/eyring.py b/chempy/kinetics/eyring.py index 334152ab..668a64c4 100644 --- a/chempy/kinetics/eyring.py +++ b/chempy/kinetics/eyring.py @@ -160,7 +160,8 @@ def equation_as_string(self, precision, tex=False): r" \exp \left(-\frac{{{}}}{{RT}} \right)").format( str_A, str_Ea + ' ' + str_Ea_unit), str_A_unit else: - return "kB*T/h*exp({}/R)*exp(-{}/(R*T))".format(str_A, str_Ea + ' ' + str_Ea_unit), str_A_unit + return "kB*T/h*exp({}/R)*exp(-{}/(R*T))".format( + str_A, str_Ea + ' ' + str_Ea_unit), str_A_unit def __str__(self): return self.equation_as_string('%.5g') @@ -175,7 +176,9 @@ def __call__(self, state, constants=default_constants, units=default_units, return super(EyringParamWithUnits, self).__call__( state, constants, units, backend) - def as_RateExpr(self, unique_keys=None, constants=default_constants, units=default_units, backend=None): + def as_RateExpr(self, unique_keys=None, constants=default_constants, + units=default_units, backend=None): if backend is None: backend = Backend() - return super(EyringParamWithUnits, self).as_RateExpr(unique_keys, constants, units, backend) + return super(EyringParamWithUnits, self).as_RateExpr( + unique_keys, constants, units, backend) diff --git a/chempy/kinetics/ode.py b/chempy/kinetics/ode.py index 460495b7..c3861f13 100644 --- a/chempy/kinetics/ode.py +++ b/chempy/kinetics/ode.py @@ -17,7 +17,9 @@ except ImportError: np = None -from ..units import to_unitless, get_derived_unit, rescale +from ..units import ( + to_unitless, get_derived_unit, rescale, magnitude, unitless_in_registry +) from ..util._expr import Expr from .rates import RateExpr, MassAction @@ -93,7 +95,7 @@ def dCdt_list(rsys, rates): def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, unit_registry=None, - output_conc_unit=None, output_time_unit=None, cstr=False, **kwargs): + output_conc_unit=None, output_time_unit=None, cstr=False, constants=None, **kwargs): """ Creates a :class:`pyneqsys.SymbolicSys` from a :class:`ReactionSystem` The parameters passed to RateExpr will contain the key ``'time'`` corresponding to the @@ -120,6 +122,9 @@ def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, output_time_unit : unit (Optional) cstr : bool Generate expressions for continuously stirred tank reactor. + constants : module + e.g. ``chempy.units.default_constants``, parameter keys not found in + substitutions will be looked for as an attribute of ``constants`` when provided. \*\*kwargs : Keyword arguemnts passed on to `SymbolicSys`. @@ -152,7 +157,6 @@ def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, from pyodesys.symbolic import SymbolicSys r_exprs = [rxn.rate_expr() for rxn in rsys.rxns] - _ori_pk = set.union(*(ratex.all_parameter_keys() for ratex in r_exprs)) _ori_uk = set.union(*(ratex.all_unique_keys() for ratex in r_exprs)) _subst_pk = set() @@ -161,7 +165,7 @@ def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, substitutions = substitutions or {} unique = OrderedDict() - unique_units = {} # OrderedDict() + unique_units = {} cstr_fr_fc = ( 'feedratio', @@ -176,15 +180,17 @@ def get_odesys(rsys, include_params=True, substitutions=None, SymbolicSys=None, def _reg_unique_unit(k, arg_dim, idx): if unit_registry is None: return - unique_units[k] = reduce(mul, [unit_registry[dim]**v for dim, v in arg_dim[idx].items()]) + unique_units[k] = reduce(mul, [1]+[unit_registry[dim]**v for dim, v in arg_dim[idx].items()]) def _reg_unique(expr, rxn=None): if not isinstance(expr, Expr): raise NotImplementedError("Currently only Expr sub classes are supported.") + if unit_registry is None: arg_dim = None else: arg_dim = expr.args_dimensionality(reaction=rxn) + if expr.args is None: for idx, k in enumerate(expr.unique_keys): if k not in substitutions: @@ -193,7 +199,7 @@ def _reg_unique(expr, rxn=None): else: for idx, arg in enumerate(expr.args): if isinstance(arg, Expr): - _reg_unique(arg) + _reg_unique(arg, rxn=rxn) elif expr.unique_keys is not None and idx < len(expr.unique_keys): uk = expr.unique_keys[idx] if uk not in substitutions: @@ -210,7 +216,20 @@ def _reg_unique(expr, rxn=None): _reg_unique(sv) else: _passive_subst[sk] = sv - all_pk = list(filter(lambda x: x not in substitutions and x != 'time', _ori_pk.union(_subst_pk))) + + all_pk = [] + for pk in filter(lambda x: x not in substitutions and x != 'time', + _ori_pk.union(_subst_pk)): + if hasattr(constants, pk): + const = getattr(constants, pk) + if unit_registry is None: + const = magnitude(const) + else: + const = unitless_in_registry(const, unit_registry) + + _passive_subst[pk] = const + else: + all_pk.append(pk) if not include_params: for rxn, ratex in zip(rsys.rxns, r_exprs): @@ -231,10 +250,7 @@ def _reg_unique(expr, rxn=None): p_units = pk_units if include_params else (pk_units + [unique_units[k] for k in unique]) new_r_exprs = [] for ratex in r_exprs: - if ratex.args is None: - _new_ratex = ratex - else: - _pu, _new_ratex = ratex.dedimensionalisation(unit_registry) + _pu, _new_ratex = ratex.dedimensionalisation(unit_registry) new_r_exprs.append(_new_ratex) r_exprs = new_r_exprs @@ -264,7 +280,7 @@ def dydt(t, y, p, backend=math): raise ValueError("Key 'time' is reserved.") variables['time'] = t for k, act in _active_subst.items(): - if unit_registry is not None: + if unit_registry is not None and act.args: _, act = act.dedimensionalisation(unit_registry) variables[k] = act(variables, backend=backend) variables.update(_passive_subst) @@ -276,7 +292,7 @@ def reaction_rates(t, y, p, backend=math): raise ValueError("Key 'time' is reserved.") variables['time'] = t for k, act in _active_subst.items(): - if unit_registry is not None: + if unit_registry is not None and act.args: _, act = act.dedimensionalisation(unit_registry) variables[k] = act(variables, backend=backend) variables.update(_passive_subst) diff --git a/chempy/kinetics/rates.py b/chempy/kinetics/rates.py index 50259c5c..32a60ef6 100644 --- a/chempy/kinetics/rates.py +++ b/chempy/kinetics/rates.py @@ -13,12 +13,15 @@ import math from operator import add -from ..units import get_derived_unit +from ..units import get_derived_unit, default_units, energy, concentration from ..util._dimensionality import dimension_codes, base_registry from ..util.pyutil import memoize, deprecated from ..util._expr import Expr +_molar = getattr(default_units, 'molar', 1) # makes module importable. + + class RateExpr(Expr): """ Baseclass for rate expressions, see source code of e.g. MassAction & Radiolytic. """ @@ -157,7 +160,7 @@ def rate_coeff(self, variables, backend=math, **kwargs): return rat_c def __call__(self, variables, backend=math, reaction=None, **kwargs): - return self.rate_coeff(variables, backend=backend)*self.active_conc_prod( + return self.rate_coeff(variables, backend=backend, reaction=reaction)*self.active_conc_prod( variables, backend=backend, reaction=reaction, **kwargs) @classmethod @@ -171,7 +174,7 @@ def string(self, *args, **kwargs): return super(MassAction, self).string(*args, **kwargs) @classmethod - @deprecated(use_instead=Expr.from_callback) + @deprecated(use_instead='MassAction.from_callback') def subclass_from_callback(cls, cb, cls_attrs=None): """ Override MassAction.__call__ """ _RateExpr = super(MassAction, cls).subclass_from_callback(cb, cls_attrs=cls_attrs) @@ -202,20 +205,61 @@ class Arrhenius(Expr): argument_names = ('A', 'Ea_over_R') parameter_keys = ('temperature',) + def args_dimensionality(self, reaction): + order = reaction.order() + return ( + {'time': -1, + 'amount': 1-order, 'length': 3*(order - 1)}, + {'temperature': 1}, + ) + def __call__(self, variables, backend=math, **kwargs): A, Ea_over_R = self.all_args(variables, backend=backend, **kwargs) return A*backend.exp(-Ea_over_R/variables['temperature']) class Eyring(Expr): - """ Rate expression for Eyring eq: c0*T*exp(-c1/T) """ + """ Rate expression for Eyring eq: c0*T*exp(-c1/T) - argument_names = ('kB_h_times_exp_dS_R', 'dH_over_R') + Note that choice of standard state (c^0) will matter for order > 1. + """ + + argument_names = ('kB_h_times_exp_dS_R', 'dH_over_R', 'conc0') + argument_defaults = (1*_molar,) + parameter_keys = ('temperature',) + + def args_dimensionality(self, reaction): + order = reaction.order() + return ( + {'time': -1, 'temperature': -1, + 'amount': 1-order, 'length': 3*(order - 1)}, + {'temperature': 1}, + concentration + ) def __call__(self, variables, backend=math, **kwargs): - c0, c1 = self.all_args(variables, backend=backend, **kwargs) + c0, c1, conc0 = self.all_args(variables, backend=backend, **kwargs) T = variables['temperature'] - return c0*T*backend.exp(-c1/T) + return c0*T*backend.exp(-c1/T)*conc0**(1-kwargs['reaction'].order()) + + +class EyringHS(Expr): + argument_names = ('dH', 'dS', 'c0') + argument_defaults = (1*_molar,) + parameter_keys = ('temperature', 'molar_gas_constant', + 'Boltzmann_constant', 'Planck_constant') + + def args_dimensionality(self, **kwargs): + return ( + energy + {'amount': -1}, + energy + {'amount': -1, 'temperature': -1}, + concentration + ) + + def __call__(self, variables, backend=math, reaction=None, **kwargs): + dH, dS, c0 = self.all_args(variables, backend=backend, **kwargs) + T, R, kB, h = [variables[k] for k in self.parameter_keys] + return kB/h*T*backend.exp(-(dH-T*dS)/(R*T))*c0**(1-reaction.order()) class RampedTemp(Expr): @@ -223,6 +267,21 @@ class RampedTemp(Expr): argument_names = ('T0', 'dTdt') parameter_keys = ('time',) # consider e.g. a parameter such as 'init_time' + def args_dimensionality(self, **kwargs): + return ({'temperature': 1}, {'temperature': 1, 'time': -1}) + def __call__(self, variables, backend=None, **kwargs): T0, dTdt = self.all_args(variables, backend=backend, **kwargs) return T0 + dTdt*variables['time'] + + +class SinTemp(Expr): + argument_names = ('Tbase', 'Tamp', 'angvel', 'phase') + parameter_keys = ('time',) + + def args_dimensionality(self, **kwargs): + return ({'temperature': 1}, {'temperature': 1}, {'time': -1}, {}) + + def __call__(self, variables, backend=math, **kwargs): + Tbase, Tamp, angvel, phase = self.all_args(variables, backend=backend, **kwargs) + return Tbase + Tamp*backend.sin(angvel*variables['time'] + phase) diff --git a/chempy/kinetics/tests/test_ode.py b/chempy/kinetics/tests/test_ode.py index c61791c8..a3ebde63 100644 --- a/chempy/kinetics/tests/test_ode.py +++ b/chempy/kinetics/tests/test_ode.py @@ -16,7 +16,7 @@ from chempy.thermodynamics.expressions import MassActionEq from chempy.units import ( SI_base_registry, get_derived_unit, allclose, units_library, - to_unitless, default_units as u + to_unitless, default_constants as const, default_units as u ) from chempy.util._expr import Expr from chempy.util.testing import requires @@ -322,7 +322,7 @@ def test_get_odesys__time_dep_rate(): class RampedRate(Expr): argument_names = ('rate_constant', 'ramping_rate') - def __call__(self, variables, backend=math): + def __call__(self, variables, reaction, backend=math): rate_constant, ramping_rate = self.all_args(variables, backend=backend) return rate_constant * ramping_rate * variables['time'] @@ -670,3 +670,243 @@ def test_get_odesys_rsys_with_units__named_params(): with pytest.raises(Exception): odesys.integrate(tend, c0, p, integrator='odeint') + + +@requires('pycvodes', 'sym', units_library) +def test_get_odesys__Eyring(): + R = 8.314472 + T_K = 300 + dH = 80e3 + dS = 10 + rsys1 = ReactionSystem.from_string(""" + NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) + """.format(dH=dH, dS=dS), substance_keys='NOBr NO Br'.split()) + kref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K)) + NOBr0_M = 0.7 + init_cond = dict( + NOBr=NOBr0_M*u.M, + NO=0*u.M, + Br=0*u.M + ) + t = 5*u.second + params = dict( + temperature=T_K*u.K + ) + + def check(rsys): + odesys, extra = get_odesys(rsys, unit_registry=SI_base_registry, constants=const) + res = odesys.integrate(t, init_cond, params, integrator='cvode') + NOBr_ref = NOBr0_M*np.exp(-kref*to_unitless(res.xout, u.second)) + ref = np.array([NOBr_ref] + [NOBr0_M - NOBr_ref]*2).T + cmp = to_unitless(res.yout, u.M) + assert np.allclose(cmp, ref) + check(rsys1) + rsys2 = ReactionSystem.from_string(""" + NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) + """.format(dH=dH, dS=dS), substance_keys='NOBr NO Br'.split()) + check(rsys2) + + +@requires('pycvodes', 'sym', units_library) +def test_get_odesys__Eyring_2nd_order(): + R = 8.314472 + T_K = 300 + dH = 80e3 + dS = 10 + rsys1b = ReactionSystem.from_string(""" + NO + Br -> NOBr; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) + """.format(dH=dH, dS=dS)) + c0 = 1 # mol/dm3 === 1000 mol/m3 + kbref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))/c0 + NO0_M = 1.5 + Br0_M = 0.7 + init_cond = dict( + NOBr=0*u.M, + NO=NO0_M*u.M, + Br=Br0_M*u.M + ) + t = 5*u.second + params = dict( + temperature=T_K*u.K + ) + + def analytic_b(t): + U, V = NO0_M, Br0_M + d = U - V + return (U*(1 - np.exp(-kbref*t*d)))/(U/V - np.exp(-kbref*t*d)) + + def check(rsys): + odesys, extra = get_odesys(rsys, unit_registry=SI_base_registry, constants=const) + res = odesys.integrate(t, init_cond, params, integrator='cvode') + t_sec = to_unitless(res.xout, u.second) + NOBr_ref = analytic_b(t_sec) + cmp = to_unitless(res.yout, u.M) + ref = np.empty_like(cmp) + ref[:, odesys.names.index('NOBr')] = NOBr_ref + ref[:, odesys.names.index('Br')] = Br0_M - NOBr_ref + ref[:, odesys.names.index('NO')] = NO0_M - NOBr_ref + assert np.allclose(cmp, ref) + + check(rsys1b) + rsys2b = ReactionSystem.from_string(""" + NO + Br -> NOBr; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) + """.format(dH=dH, dS=dS)) + check(rsys2b) + + +@requires('pycvodes', 'sym', 'scipy', units_library) +def test_get_odesys__Eyring_1st_order_linearly_ramped_temperautre(): + from scipy.special import expi + + def analytic_unit0(t, T0, dH, dS): + R = 8.314472 + kB = 1.3806504e-23 + h = 6.62606896e-34 + A = kB/h*np.exp(dS/R) + B = dH/R + return np.exp(A*( + (-B**2*np.exp(B/T0)*expi(-B/T0) - T0*(B - T0))*np.exp(-B/T0) + + (B**2*np.exp(B/(t + T0))*expi(-B/(t + T0)) - (t + T0)*(-B + t + T0))*np.exp(-B/(t + T0)) + )/2) + + T_K = 290 + dH = 80e3 + dS = 10 + rsys1 = ReactionSystem.from_string(""" + NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) + """.format(dH=dH, dS=dS)) + + NOBr0_M = 0.7 + init_cond = dict( + NOBr=NOBr0_M*u.M, + NO=0*u.M, + Br=0*u.M + ) + t = 20*u.second + + def check(rsys): + odes, extra = get_odesys(rsys, unit_registry=SI_base_registry, constants=const, substitutions={ + 'temperature': RampedTemp([T_K*u.K, 1*u.K/u.s])}) + for odesys in [odes, odes.as_autonomous()]: + res = odesys.integrate(t, init_cond, integrator='cvode') + t_sec = to_unitless(res.xout, u.second) + NOBr_ref = NOBr0_M*analytic_unit0(t_sec, T_K, dH, dS) + cmp = to_unitless(res.yout, u.M) + ref = np.empty_like(cmp) + ref[:, odesys.names.index('NOBr')] = NOBr_ref + ref[:, odesys.names.index('Br')] = NOBr0_M - NOBr_ref + ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref + assert np.allclose(cmp, ref) + + check(rsys1) + rsys2 = ReactionSystem.from_string(""" + NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) + """.format(dH=dH, dS=dS)) + check(rsys2) + + +@requires('pycvodes', 'sym', 'scipy', units_library) +def test_get_odesys__Eyring_2nd_order_linearly_ramped_temperautre(): + from scipy.special import expi + + def analytic_unit0(t, k, m, dH, dS): + R = 8.314472 + kB = 1.3806504e-23 + h = 6.62606896e-34 + A = kB/h*np.exp(dS/R) + B = dH/R + return k*np.exp(B*(k*t + 2*m)/(m*(k*t + m)))/( + A*(-B**2*np.exp(B/(k*t + m))*expi(-B/(k*t + m)) - B*k*t - B*m + k**2*t**2 + 2*k*m*t + m**2)*np.exp(B/m) + + (A*B**2*np.exp(B/m)*expi(-B/m) - A*m*(-B + m) + k*np.exp(B/m))*np.exp(B/(k*t + m)) + ) + + T_K = 290 + dTdt_Ks = 3 + dH = 80e3 + dS = 10 + rsys1 = ReactionSystem.from_string(""" + 2 NO2 -> N2O4; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol) + """.format(dH=dH, dS=dS)) + + NO2_M = 1.0 + init_cond = dict( + NO2=NO2_M*u.M, + N2O4=0*u.M + ) + t = 20*u.second + + def check(rsys): + odes, extra = get_odesys(rsys, unit_registry=SI_base_registry, constants=const, substitutions={ + 'temperature': RampedTemp([T_K*u.K, dTdt_Ks*u.K/u.s])}) + for odesys in [odes, odes.as_autonomous()]: + res = odesys.integrate(t, init_cond, integrator='cvode') + t_sec = to_unitless(res.xout, u.second) + NO2_ref = analytic_unit0(t_sec, dTdt_Ks, T_K, dH, dS) + cmp = to_unitless(res.yout, u.M) + ref = np.empty_like(cmp) + ref[:, odesys.names.index('NO2')] = NO2_ref + ref[:, odesys.names.index('N2O4')] = (NO2_M - NO2_ref)/2 + assert np.allclose(cmp, ref) + + check(rsys1) + rsys2 = ReactionSystem.from_string(""" + 2 NO2 -> N2O4; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol])) + """.format(dH=dH, dS=dS)) + check(rsys2) + + +@requires('pycvodes', 'sym', units_library) +def test_get_odesys__Eyring_2nd_order_reversible(): + R = 8.314472 + T_K = 273.15 + 20 # 20 degree celsius + kB = 1.3806504e-23 + h = 6.62606896e-34 + + dHf = 74e3 + dSf = R*np.log(h/kB/T_K*1e16) + + dHb = 79e3 + dSb = dSf - 23 + + rsys1 = ReactionSystem.from_string(""" + Fe+3 + SCN- -> FeSCN+2; EyringParam(dH={dHf}*J/mol, dS={dSf}*J/K/mol) + FeSCN+2 -> Fe+3 + SCN-; EyringParam(dH={dHb}*J/mol, dS={dSb}*J/K/mol) + """.format(dHf=dHf, dSf=dSf, dHb=dHb, dSb=dSb)) + kf_ref = 20836643994.118652*T_K*np.exp(-(dHf - T_K*dSf)/(R*T_K)) + kb_ref = 20836643994.118652*T_K*np.exp(-(dHb - T_K*dSb)/(R*T_K)) + Fe0 = 6e-3 + SCN0 = 2e-3 + init_cond = { + 'Fe+3': Fe0*u.M, + 'SCN-': SCN0*u.M, + 'FeSCN+2': 0*u.M + } + t = 3*u.second + + def check(rsys, params): + odes, extra = get_odesys(rsys, include_params=False, + unit_registry=SI_base_registry, constants=const) + for odesys in [odes, odes.as_autonomous()]: + res = odesys.integrate(t, init_cond, params, integrator='cvode') + t_sec = to_unitless(res.xout, u.second) + FeSCN_ref = binary_rev(t_sec, kf_ref, kb_ref, 0, Fe0, SCN0) + cmp = to_unitless(res.yout, u.M) + ref = np.empty_like(cmp) + ref[:, odesys.names.index('FeSCN+2')] = FeSCN_ref + ref[:, odesys.names.index('Fe+3')] = Fe0 - FeSCN_ref + ref[:, odesys.names.index('SCN-')] = SCN0 - FeSCN_ref + assert np.allclose(cmp, ref) + + check(rsys1, {'temperature': T_K*u.K}) + rsys2 = ReactionSystem.from_string(""" + Fe+3 + SCN- -> FeSCN+2; MassAction(EyringHS([{dHf}*J/mol, {dSf}*J/K/mol])) + FeSCN+2 -> Fe+3 + SCN-; MassAction(EyringHS([{dHb}*J/mol, {dSb}*J/K/mol])) + """.format(dHf=dHf, dSf=dSf, dHb=dHb, dSb=dSb)) + check(rsys2, {'temperature': T_K*u.K}) + rsys3 = ReactionSystem.from_string(""" + Fe+3 + SCN- -> FeSCN+2; MassAction(EyringHS.fk('dHf', 'dSf')) + FeSCN+2 -> Fe+3 + SCN-; MassAction(EyringHS.fk('dHb', 'dSb')) + """) + check(rsys3, dict(temperature=T_K*u.K, + dHf=dHf*u.J/u.mol, dSf=dSf*u.J/u.mol/u.K, + dHb=dHb*u.J/u.mol, dSb=dSb*u.J/u.mol/u.K)) diff --git a/chempy/kinetics/tests/test_rates.py b/chempy/kinetics/tests/test_rates.py index ad80f3ab..31997ba6 100644 --- a/chempy/kinetics/tests/test_rates.py +++ b/chempy/kinetics/tests/test_rates.py @@ -82,7 +82,7 @@ def test_MassAction(): def test_Expr__from_callback(): - def rate_coeff(args, T, backend=math): + def rate_coeff(args, T, reaction, backend=math): return args[0]*backend.exp(args[1]/T) RateExpr = Expr.from_callback(rate_coeff, parameter_keys=('Tem',), nargs=2) k1 = RateExpr([2.1e10, -5132.2]) @@ -94,7 +94,7 @@ def rate_coeff(args, T, backend=math): def test_MassAction__subclass_from_callback(): - def rate_coeff(variables, all_args, backend): + def rate_coeff(variables, all_args, reaction, backend): return all_args[0]*backend.exp(all_args[1]/variables['temperature']) CustomMassAction = MassAction.subclass_from_callback( rate_coeff, cls_attrs=dict(parameter_keys=('temperature',), nargs=2)) @@ -108,7 +108,7 @@ def rate_coeff(variables, all_args, backend): @requires(units_library) def test_MassAction__subclass_from_callback__units(): - def rate_coeff(variables, all_args, backend): + def rate_coeff(variables, all_args, backend, **kwargs): return all_args[0]*backend.exp(all_args[1]/variables['temperature']) CustomMassAction = MassAction.subclass_from_callback( rate_coeff, cls_attrs=dict(parameter_keys=('temperature',), nargs=2)) @@ -142,7 +142,7 @@ def ref(v): assert abs((ma(var, reaction=r) - ref_var)/ref_var) < 1e-14 with pytest.raises(ValueError): - Arrhenius([A, Ea_over_R, A]) + Arrhenius([A, Ea_over_R, 1, A]) # assert ama.as_mass_action({T_: 273.15}).args[0] == A*math.exp(-Ea_over_R/273.15) @@ -255,14 +255,14 @@ def test_mk_Radiolytic(): class Eyring4(Expr): nargs = 4 - def __call__(self, variables, backend=math): + def __call__(self, variables, reaction, backend=math): Sact_fact, Hact_exp, Sref_fact, Href_exp = self.all_args(variables, backend=backend) T = variables['temperature'] return T * Sact_fact / Sref_fact * backend.exp((Href_exp-Hact_exp)/T) def test_EyringMassAction(): - args = kB_h_times_exp_dS_R, dH_over_R = 1.2e11/273.15, 40e3/8 + args = kB_h_times_exp_dS_R, dH_over_R, c0 = 1.2e11/273.15, 40e3/8, 1 ama = MassAction(Eyring(args, ('Sfreq', 'Hact'))) rxn1 = Reaction({'A': 2, 'B': 1}, {'C': 1}, ama, {'B': 1}) T_ = 'temperature' diff --git a/chempy/tests/test_reactionsystem.py b/chempy/tests/test_reactionsystem.py index 1c65035f..7574bb86 100644 --- a/chempy/tests/test_reactionsystem.py +++ b/chempy/tests/test_reactionsystem.py @@ -211,6 +211,7 @@ def test_ReactionSystem__upper_conc_bounds__a_substance_no_composition(): assert res == ref +@requires(parsing_library) def test_ReactionSystem__identify_equilibria(): rsys = ReactionSystem.from_string(""" 2 H2 + O2 -> 2 H2O ; 1e-3 diff --git a/chempy/units.py b/chempy/units.py index de7f3d9d..c12d9e74 100644 --- a/chempy/units.py +++ b/chempy/units.py @@ -21,6 +21,7 @@ from operator import mul from ._util import NameSpace +from .util.arithmeticdict import ArithmeticDict units_library = 'quantities' # info used for selective testing. @@ -108,6 +109,10 @@ def magnitude(value): except AttributeError: return value +energy = ArithmeticDict(int, {'mass': 1, 'length': 2, 'time': -2}) +volume = ArithmeticDict(int, {'length': 3}) +concentration = {'amount': 1} - volume + def get_derived_unit(registry, key): """ Get the unit of a physcial quantity in a provided unit system. diff --git a/chempy/util/_dimensionality.py b/chempy/util/_dimensionality.py index 5ece41c7..e83041a8 100644 --- a/chempy/util/_dimensionality.py +++ b/chempy/util/_dimensionality.py @@ -5,8 +5,8 @@ from .pyutil import defaultnamedtuple dimension_codes = OrderedDict(zip( - 'length mass time current temperature luminous_intensity amount'.split(), - 'L M T I Θ N J'.split() + 'length mass time current temperature amount'.split(), # not considering luminous_intensity + 'L M T I Θ N'.split() )) diff --git a/chempy/util/_expr.py b/chempy/util/_expr.py index fa284a2c..2f48dfd0 100644 --- a/chempy/util/_expr.py +++ b/chempy/util/_expr.py @@ -128,6 +128,11 @@ def __init__(self, args=None, unique_keys=None): self.args = args + @classmethod + def fk(cls, *args): + """ Alternative constructor "from keys", \\*args is used as ``unique_keys``. """ + return cls(unique_keys=args) + @classmethod def from_callback(cls, callback, attr='__call__', **kwargs): """ Factory of subclasses @@ -159,13 +164,14 @@ def from_callback(cls, callback, attr='__call__', **kwargs): True """ - def body(self, variables, backend=math, **kwargs): + def body(self, variables, backend=math, **kw): args = self.all_args(variables, backend=backend) params = self.all_params(variables, backend=backend) - return callback(args, *params, backend=backend, **kwargs) + return callback(args, *params, backend=backend, **kw) class Wrapper(cls): pass + setattr(Wrapper, attr, body) Wrapper.__name__ = callback.__name__ for k, v in kwargs.items(): @@ -305,20 +311,29 @@ def dedimensionalisation(self, unit_registry, variables={}, backend=math): self.__class__ instance: with dedimensioanlised arguments """ - from ..units import default_unit_in_registry, to_unitless - units = [None if isinstance(arg, Expr) else default_unit_in_registry(arg, unit_registry) for arg - in self.all_args(variables, backend=backend, evaluate=False)] - new_units, unitless_args = [], [] - for arg, unit in zip(self.all_args(variables, backend=backend, evaluate=False), units): - if isinstance(arg, Expr): - if unit is not None: - raise ValueError() - _unit, _dedim = arg.dedimensionalisation(unit_registry, variables, backend=backend) - else: - _unit, _dedim = unit, to_unitless(arg, unit) - new_units.append(_unit) - unitless_args.append(_dedim) - return new_units, self.__class__(unitless_args, self.unique_keys) + from ..units import default_unit_in_registry, to_unitless, unitless_in_registry + new_units = [] + if self.args is None: + unitless_args = None + else: + unitless_args = [] + units = [None if isinstance(arg, Expr) else default_unit_in_registry(arg, unit_registry) for arg + in self.all_args(variables, backend=backend, evaluate=False)] + for arg, unit in zip(self.all_args(variables, backend=backend, evaluate=False), units): + if isinstance(arg, Expr): + if unit is not None: + raise ValueError() + _unit, _dedim = arg.dedimensionalisation(unit_registry, variables, backend=backend) + else: + _unit, _dedim = unit, to_unitless(arg, unit) + new_units.append(_unit) + unitless_args.append(_dedim) + instance = self.__class__(unitless_args, self.unique_keys) + if self.argument_defaults is not None: + instance.argument_defaults = tuple(unitless_in_registry(arg, unit_registry) + for arg in self.argument_defaults) + + return new_units, instance def _sympy_format(self, method, variables, backend, default): variables = variables or {} @@ -499,6 +514,10 @@ class Log10(UnaryFunction): _func_name = 'log10' +class Exp(UnaryFunction): + _func_name = 'exp' + + def create_Piecewise(parameter_name, nan_fallback=False): """ Examples diff --git a/chempy/util/arithmeticdict.py b/chempy/util/arithmeticdict.py index 0f693435..ab5fd7b1 100644 --- a/chempy/util/arithmeticdict.py +++ b/chempy/util/arithmeticdict.py @@ -184,10 +184,10 @@ def all_non_negative(self): return False return True - def __hash__(self): - default = self.default_factory() - l = [self.default_factory] - for k, v in self.items(): - if v != default: - l.append((k, v)) - return hash(tuple(l)) + # def __hash__(self): + # default = self.default_factory() + # l = [self.default_factory] + # for k, v in self.items(): + # if v != default: + # l.append((k, v)) + # return hash(tuple(l)) diff --git a/chempy/util/parsing.py b/chempy/util/parsing.py index 52713bc9..402b120b 100644 --- a/chempy/util/parsing.py +++ b/chempy/util/parsing.py @@ -334,6 +334,8 @@ def to_reaction(line, substance_keys, token, Cls, globals_=None, **kwargs): globals_ = {k: getattr(rates, k) for k in dir(rates)} globals_.update({'chempy': chempy, 'default_units': default_units}) if default_units is not None: + globals_.update({k: v for k, v in chempy.__dict__.items() + if not k.startswith('_')}) globals_.update(default_units.as_dict()) try: stoich, param, kw = map(str.strip, line.rstrip('\n').split(';')) diff --git a/examples/Ammonia.ipynb b/examples/Ammonia.ipynb index 0b1668d5..a9c317a4 100644 --- a/examples/Ammonia.ipynb +++ b/examples/Ammonia.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from __future__ import division, print_function\n", @@ -27,11 +23,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "substance_names = ['H+', 'OH-', 'NH4+', 'NH3', 'H2O']\n", @@ -49,12 +41,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "rs = EqSystem(equilibria, subst)\n", @@ -65,12 +52,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "logx, logsol, sane = rs.root(init_conc, x0=x, NumSys=(NumSysLog,))\n", @@ -80,11 +62,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x - logx" @@ -93,11 +71,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "ny = len(substance_names)\n", @@ -113,11 +87,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "numsys_log = NumSysLog(rs, backend=sp)\n", @@ -128,11 +98,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "numsys_lin = NumSysLin(rs, backend=sp)\n", @@ -142,11 +108,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "A, ks = rs.stoichs_constants(False, backend=sp)\n", @@ -156,11 +118,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyneqsys.symbolic import SymbolicSys\n", @@ -173,11 +131,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "j = sp.Matrix(1, len(f), lambda _, q: f[q]).jacobian(y)\n", @@ -192,11 +146,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "j.simplify()\n", @@ -206,11 +156,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rs.composition_balance_vectors()" @@ -219,11 +165,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "numsys_rref_log = NumSysLog(rs, True, True, backend=sp)\n", @@ -233,11 +175,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "np.set_printoptions(4, linewidth=120)\n", @@ -249,12 +187,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "x, res, sane = rs.root(init_conc, rref_equil=True, rref_preserv=True)\n", @@ -264,11 +197,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x, res, sane = rs.root(init_conc, x0=rs.as_per_substance_array(\n", @@ -279,11 +208,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x, res, sane = rs.root(init_conc, x0=rs.as_per_substance_array(\n", @@ -294,11 +219,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "init_conc" @@ -307,11 +228,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "nc=60\n", @@ -332,11 +249,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_lin = plot_rref(method='lm')\n", @@ -346,12 +259,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "for col_id in range(len(substance_names)):\n", @@ -365,11 +273,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rs.plot_errors(res_lin[0][0], init_conc, Hp_0, 'H+')" @@ -378,11 +282,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "init_conc, rs.ns" @@ -391,11 +291,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_log = plot_rref(NumSys=NumSysLog)" @@ -404,11 +300,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rs.plot_errors(res_log[0][0], init_conc, Hp_0, 'H+')" @@ -417,11 +309,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_log_lin = plot_rref(NumSys=(NumSysLog, NumSysLin))\n", @@ -431,11 +319,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.equilibria import NumSysSquare\n", @@ -446,11 +330,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_sq = plot_rref(NumSys=(NumSysSquare,), method='lm')\n", @@ -460,11 +340,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x, res, sane = rs.root(x0, NumSys=NumSysLog, rref_equil=True, rref_preserv=True)\n", @@ -474,11 +350,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x, res, sane = rs.root(x, NumSys=NumSysLin, rref_equil=True, rref_preserv=True)\n", @@ -488,11 +360,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# x, res, sane = rs.root(x, NumSys=(NumSysLinTanh,), rref_equil=False, rref_preserv=False)\n", @@ -502,11 +370,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# res_tanh = plot_rref(NumSys=(NumSysLog, NumSysLinTanh,))\n", diff --git a/examples/NaCl_precipitation.ipynb b/examples/NaCl_precipitation.ipynb index 051f31a5..80d0f8ab 100644 --- a/examples/NaCl_precipitation.ipynb +++ b/examples/NaCl_precipitation.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.chemistry import Species, Equilibrium\n", @@ -17,11 +13,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "Na_p = Species('Na+', 1, composition={11: 1})\n", @@ -34,11 +26,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys.root({'Na+': 5, 'Cl-': 5, 'NaCl': 0})" @@ -47,11 +35,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys.root({'Na+': 1, 'Cl-': 1, 'NaCl': 0})" @@ -60,11 +44,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys.root({'Na+': 2, 'Cl-': 2, 'NaCl': 0})" @@ -73,11 +53,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys.root({'Na+': 0, 'Cl-': 0, 'NaCl': 2})" @@ -86,11 +62,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys.root({'Na+': 0, 'Cl-': 0, 'NaCl': 5.0})" @@ -99,12 +71,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys.root({'Na+': 0, 'Cl-': 0, 'NaCl': 5.0}, rref_preserv=True)" @@ -113,11 +80,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.equilibria import NumSysLog\n", @@ -127,11 +90,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.equilibria import NumSysLin\n", @@ -141,11 +100,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] } diff --git a/examples/Robertson_kinetics.ipynb b/examples/Robertson_kinetics.ipynb index fb6b643b..bbd5fb9e 100644 --- a/examples/Robertson_kinetics.ipynb +++ b/examples/Robertson_kinetics.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from __future__ import absolute_import, division, print_function\n", @@ -26,11 +22,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "A, B, C, D = map(Substance, 'ABCD')\n", @@ -47,11 +39,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from IPython.display import HTML\n", @@ -61,11 +49,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys2graph(rsys, 'robertson.png', save='.')\n", @@ -75,11 +59,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys, extra = get_odesys(rsys, include_params=True)\n", @@ -89,11 +69,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys.get_jac()" @@ -102,11 +78,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "c0 = defaultdict(float, {'A': 1})\n", @@ -117,12 +89,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "extra['rate_exprs_cb'](result.xout, result.yout)" @@ -131,12 +98,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "result.plot(xscale='log', yscale='log')" @@ -145,12 +107,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 2, figsize=(14, 6))\n", @@ -173,11 +130,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "str_massaction = \"\"\"\n", @@ -190,11 +143,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys3 = ReactionSystem.from_string(str_massaction, substance_factory=lambda formula: Substance(formula))" @@ -203,11 +152,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys3.substance_names()" @@ -216,12 +161,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys3, extra3 = get_odesys(rsys3, include_params=False, lower_bounds=[0, 0, 0])\n", @@ -231,11 +171,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys3.exprs, odesys3.params, odesys3.names, odesys3.param_names" @@ -244,11 +180,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot(A0=1.0, B0=0.0, C0=0.0, lg_k1=-2, lg_k2=4, lg_k3=7, lg_tend=9):\n", @@ -269,11 +201,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# We could also have used SymPy to construct symbolic rates:\n", @@ -283,7 +211,7 @@ "B + C -> A + C; sp.Symbol('k2')\n", "2 B -> B + C; sp.Symbol('k3')\n", "\"\"\", rxn_parse_kwargs=dict(globals_={'sp': sympy}), substance_factory=lambda formula: Substance(formula))\n", - "odesys_sym, _ = get_odesys(rsys_sym)\n", + "odesys_sym, _ = get_odesys(rsys_sym, params=True)\n", "for attr in 'exprs params names param_names'.split():\n", " print(getattr(odesys_sym, attr))" ] @@ -291,12 +219,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "# For larger systems it is easy to loose track of what substances are actually playing a part, here the html table can help:\n", @@ -305,6 +228,13 @@ "rsys.substances.pop('D')\n", "HTML(_html)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -327,7 +257,7 @@ }, "widgets": { "state": { - "75650d88adf74eae995032910a0f4866": { + "6296ca2fc06a4f0c9d9390c4dd51a493": { "views": [ { "cell_index": 16 diff --git a/examples/_4state_model.ipynb b/examples/_4state_model.ipynb index 61c43f59..ead3c966 100644 --- a/examples/_4state_model.ipynb +++ b/examples/_4state_model.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from collections import OrderedDict, defaultdict\n", @@ -29,11 +25,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "substances = OrderedDict([\n", @@ -48,11 +40,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def _gibbs(args, T, R, backend, **kwargs):\n", @@ -69,11 +57,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "Gibbs = EqExpr.from_callback(_gibbs, parameter_keys=('temperature', 'R'), argument_names=('H', 'S', 'Cp', 'Tref'))\n", @@ -83,11 +67,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "thermo_dis = Gibbs(unique_keys=('He_dis', 'Se_dis', 'Cp_dis', 'Tref_dis'))\n", @@ -100,11 +80,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "eq_dis = Equilibrium({'NL'}, {'N', 'L'}, thermo_dis, name='ligand-protein dissociation')\n", @@ -115,11 +91,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys = ReactionSystem(\n", @@ -132,11 +104,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "vecs, comp = rsys.composition_balance_vectors()\n", @@ -147,11 +115,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys2graph(rsys, '4state.png', save='.', include_inactive=False)\n", @@ -161,12 +125,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "from IPython.display import HTML\n", @@ -177,11 +136,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def pretty_replace(s, subs=None):\n", @@ -218,11 +173,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "ratexs = [autosymbols['r(\\mathrm{%s})' % rnames[rxn.name]] for rxn in rsys.rxns]\n", @@ -234,11 +185,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "default_c0 = defaultdict(float, {'N': 1e-9, 'L': 1e-8})\n", @@ -267,10 +214,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "We have the melting temperature $T_m$ as a free parameter, however, the model is expressed in terms of $\\Delta_u S ^\\circ$ so will need to derive the latter from the former:\n", "$$\n", @@ -298,11 +242,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def Se0_from_Tm(Tm, token):\n", @@ -315,11 +255,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "params_c0 = default_c0.copy()\n", @@ -331,11 +267,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.symbolic import ScaledSys" @@ -344,11 +276,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys, extra = get_odesys(rsys, include_params=False, SymbolicSys=ScaledSys, dep_scaling=1e9)\n", @@ -359,12 +287,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot(system, c0=None, first_step=None, t0=0, stiffness=False, nsteps=9000, **kwargs):\n", @@ -404,11 +327,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.native.cvode import NativeCvodeSys\n", @@ -418,12 +337,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "_, _, info_native = integrate_and_plot(native)\n", @@ -433,11 +347,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "info['time_wall']/info_native['time_wall']" @@ -446,11 +356,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.kinetics._native import get_native" @@ -459,11 +365,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "native2 = get_native(rsys, odesys, 'cvode')" @@ -472,11 +374,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "_, _, info_native2 = integrate_and_plot(native2, first_step=0.0)\n", @@ -485,10 +383,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "We have one complication, due to linear dependencies in our formulation of the system of ODEs our jacobian is singular:" ] @@ -496,11 +391,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "cses, (jac_in_cse,) = odesys.be.cse(odesys.get_jac())\n", @@ -510,11 +401,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys.jacobian_singular()" @@ -523,9 +410,7 @@ { "cell_type": "markdown", "metadata": { - "collapsed": true, - "deletable": true, - "editable": true + "collapsed": true }, "source": [ "Since implicit methods (which are required for stiff cases often encountered in kinetic modelling) uses the jacboian in the modified Newton's method we may get failures during integration. What we can do is to identify linear dependencies based on composition of the materials and exploit the invariants to reduce the dimensionality of the system of ODEs:" @@ -534,11 +419,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "A, comp_names = rsys.composition_balance_vectors()\n", @@ -548,11 +429,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "y0 = {odesys[k]: sympy.Symbol(k+'0') for k in rsys.substances.keys()}\n", @@ -565,11 +442,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.symbolic import PartiallySolvedSystem\n", @@ -582,11 +455,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psysLA = PartiallySolvedSystem(odesys, extra['linear_dependencies'](['L', 'A']), **no_invar)\n", @@ -597,12 +466,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(12,4))\n", @@ -621,11 +485,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.symbolic import SymbolicSys\n", @@ -635,11 +495,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psys_root.roots" @@ -648,11 +504,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psysLN['N']" @@ -661,11 +513,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psysLN.analytic_exprs" @@ -674,11 +522,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psysLN.names" @@ -687,11 +531,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psysLN.dep" @@ -700,11 +540,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "tout1, Cout1, info_root = integrate_and_plot(psys_root, first_step=0.0, return_on_root=True)" @@ -713,11 +549,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "xout2, yout2, info_LA = integrate_and_plot(psysLA, first_step=0.0, t0=tout1[-1], c0=dict(zip(odesys.names, Cout1[-1, :])))" @@ -726,11 +558,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "print('\\troot\\tLA\\troot+LA\\tLN')\n", @@ -741,11 +569,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.symbolic import symmetricsys\n", @@ -766,11 +590,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "_, _, info_t = integrate_and_plot(tsys, first_step=0.0)\n", @@ -780,12 +600,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "_, _, info_tLN = integrate_and_plot(tsysLN, first_step=0.0, nsteps=18000)\n", @@ -795,11 +610,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "_, _, info_tLA = integrate_and_plot(tsysLA, first_step=0.0)\n", @@ -809,11 +620,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "print(open(next(filter(lambda s: s.endswith('.cpp'), native2._native._written_files))).read())" @@ -836,7 +643,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.1" } }, "nbformat": 4, diff --git a/examples/_eyring_kinetics.ipynb b/examples/_eyring_kinetics.ipynb new file mode 100644 index 00000000..aa91cfbe --- /dev/null +++ b/examples/_eyring_kinetics.ipynb @@ -0,0 +1,126 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from chempy import ReactionSystem\n", + "from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const\n", + "from chempy.kinetics.ode import get_odesys\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = 8.314472\n", + "T_K = 300\n", + "dH=80e3\n", + "dS=10\n", + "rsys1 = ReactionSystem.from_string(\"\"\"\n", + "NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n", + "\"\"\".format(dH=dH, dS=dS))\n", + "kref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))\n", + "kref" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NOBr0_M = 0.7\n", + "init_cond = dict(\n", + " NOBr=NOBr0_M*u.M,\n", + " NO=0*u.M,\n", + " Br=0*u.M\n", + ")\n", + "t = 5*u.second\n", + "params = dict(\n", + " temperature=T_K*u.K\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(rsys):\n", + " odesys, extra = get_odesys(rsys, unit_registry=si, constants=const)\n", + " fig, axes = plt.subplots(1, 3, figsize=(14, 6))\n", + " res = odesys.integrate(t, init_cond, params, integrator='cvode')\n", + " t_sec = to_unitless(res.xout, u.second)\n", + " NOBr_ref = NOBr0_M*np.exp(-kref*t_sec)\n", + " cmp = to_unitless(res.yout, u.M)\n", + " ref = np.empty_like(cmp)\n", + " ref[:, odesys.names.index('NOBr')] = NOBr_ref\n", + " ref[:, odesys.names.index('Br')] = NOBr0_M - NOBr_ref\n", + " ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref\n", + " assert np.allclose(cmp, ref)\n", + " axes[0].plot(t_sec, cmp)\n", + " axes[1].plot(t_sec, cmp - ref)\n", + " res.plot_invariant_violations(ax=axes[2])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys2 = ReactionSystem.from_string(\"\"\"\n", + "NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))\n", + "\"\"\".format(dH=dH, dS=dS))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys2)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_eyring_kinetics_2nd_order.ipynb b/examples/_eyring_kinetics_2nd_order.ipynb new file mode 100644 index 00000000..98d66984 --- /dev/null +++ b/examples/_eyring_kinetics_2nd_order.ipynb @@ -0,0 +1,149 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from chempy import ReactionSystem\n", + "from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const\n", + "from chempy.kinetics.ode import get_odesys\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = 8.314472\n", + "T_K = 300\n", + "dH=80e3\n", + "dS=10\n", + "rsys1b = ReactionSystem.from_string(\"\"\"\n", + "NO + Br -> NOBr; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n", + "\"\"\".format(dH=dH, dS=dS))\n", + "c0 = 1 # mol/dm3 === 1000 mol/m3\n", + "kbref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))/c0\n", + "kbref" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NO0_M = 1.5\n", + "Br0_M = 0.7\n", + "init_cond = dict(\n", + " NOBr=0*u.M,\n", + " NO=NO0_M*u.M,\n", + " Br=Br0_M*u.M\n", + ")\n", + "t = 5*u.second\n", + "params = dict(\n", + " temperature=T_K*u.K\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def analytic_b(t):\n", + " U, V = NO0_M, Br0_M\n", + " d = U - V\n", + " return (U*(1 - np.exp(-kbref*t*d)))/(U/V - np.exp(-kbref*t*d))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(rsys):\n", + " odesys, extra = get_odesys(rsys, unit_registry=si, constants=const)\n", + " fig, axes = plt.subplots(1, 4, figsize=(14, 6))\n", + " res = odesys.integrate(t, init_cond, params, integrator='cvode')\n", + " t_sec = to_unitless(res.xout, u.second)\n", + " NOBr_ref = analytic_b(t_sec)\n", + " cmp = to_unitless(res.yout, u.M)\n", + " ref = np.empty_like(cmp)\n", + " ref[:, odesys.names.index('NOBr')] = NOBr_ref\n", + " ref[:, odesys.names.index('Br')] = Br0_M - NOBr_ref\n", + " ref[:, odesys.names.index('NO')] = NO0_M - NOBr_ref\n", + " axes[0].plot(t_sec, cmp)\n", + " axes[1].plot(t_sec, ref)\n", + " axes[2].plot(t_sec, cmp - ref)\n", + " res.plot_invariant_violations(ax=axes[3])\n", + " assert np.allclose(cmp, ref)\n", + " print({k: v for k, v in res.info.items() if not k.startswith('internal')}) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys1b)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys2b = ReactionSystem.from_string(\"\"\"\n", + "NO + Br -> NOBr; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))\n", + "\"\"\".format(dH=dH, dS=dS))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys2b)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_eyring_kinetics_2nd_order_reversible.ipynb b/examples/_eyring_kinetics_2nd_order_reversible.ipynb new file mode 100644 index 00000000..5ff237cf --- /dev/null +++ b/examples/_eyring_kinetics_2nd_order_reversible.ipynb @@ -0,0 +1,162 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sm\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from chempy import ReactionSystem\n", + "from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.integrated import binary_rev\n", + "sm.init_printing()\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = 8.314472\n", + "T_K = 273.15 + 20 # 20 degree celsius\n", + "kB = 1.3806504e-23\n", + "h = 6.62606896e-34\n", + "\n", + "dHf = 74e3\n", + "dSf = R*np.log(h/kB/T_K*1e16)\n", + "\n", + "dHb = 79e3\n", + "dSb = dSf - 23\n", + "\n", + "rsys1 = ReactionSystem.from_string(\"\"\"\n", + "Fe+3 + SCN- -> FeSCN+2; EyringParam(dH={dHf}*J/mol, dS={dSf}*J/K/mol)\n", + "FeSCN+2 -> Fe+3 + SCN-; EyringParam(dH={dHb}*J/mol, dS={dSb}*J/K/mol)\n", + "\"\"\".format(dHf=dHf, dSf=dSf, dHb=dHb, dSb=dSb))\n", + "kf_ref = 20836643994.118652*T_K*np.exp(-(dHf - T_K*dSf)/(R*T_K))\n", + "kb_ref = 20836643994.118652*T_K*np.exp(-(dHb - T_K*dSb)/(R*T_K))\n", + "kf_ref, kb_ref" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Fe0 = 6e-3\n", + "SCN0 = 2e-3\n", + "init_cond = {\n", + " 'Fe+3': Fe0*u.M,\n", + " 'SCN-': SCN0*u.M,\n", + " 'FeSCN+2': 0*u.M\n", + "}\n", + "t = 3*u.second" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(rsys, params):\n", + " odes, extra = get_odesys(rsys, include_params=False, unit_registry=si, constants=const)\n", + " fig, all_axes = plt.subplots(2, 3, figsize=(14, 6))\n", + " for axes, odesys in zip(all_axes, [odes, odes.as_autonomous()]):\n", + " res = odesys.integrate(t, init_cond, params, integrator='cvode')\n", + " t_sec = to_unitless(res.xout, u.second)\n", + " FeSCN_ref = binary_rev(t_sec, kf_ref, kb_ref, 0, Fe0, SCN0)\n", + " cmp = to_unitless(res.yout, u.M)\n", + " ref = np.empty_like(cmp)\n", + " ref[:, odesys.names.index('FeSCN+2')] = FeSCN_ref\n", + " ref[:, odesys.names.index('Fe+3')] = Fe0 - FeSCN_ref\n", + " ref[:, odesys.names.index('SCN-')] = SCN0 - FeSCN_ref\n", + " axes[0].plot(t_sec, cmp)\n", + " axes[1].plot(t_sec, cmp - ref)\n", + " res.plot_invariant_violations(ax=axes[2])\n", + " assert np.allclose(cmp, ref)\n", + " print({k: v for k, v in res.info.items() if not k.startswith('internal')}) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys1, {'temperature': T_K*u.K})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys2 = ReactionSystem.from_string(\"\"\"\n", + "Fe+3 + SCN- -> FeSCN+2; MassAction(EyringHS([{dHf}*J/mol, {dSf}*J/K/mol]))\n", + "FeSCN+2 -> Fe+3 + SCN-; MassAction(EyringHS([{dHb}*J/mol, {dSb}*J/K/mol]))\n", + "\"\"\".format(dHf=dHf, dSf=dSf, dHb=dHb, dSb=dSb))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys2, {'temperature': T_K*u.K})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys3 = ReactionSystem.from_string(\"\"\"\n", + "Fe+3 + SCN- -> FeSCN+2; MassAction(EyringHS.fk('dHf', 'dSf'))\n", + "FeSCN+2 -> Fe+3 + SCN-; MassAction(EyringHS.fk('dHb', 'dSb'))\n", + "\"\"\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys3, dict(temperature=T_K*u.K,\n", + " dHf=dHf*u.J/u.mol, dSf=dSf*u.J/u.mol/u.K,\n", + " dHb=dHb*u.J/u.mol, dSb=dSb*u.J/u.mol/u.K))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_eyring_kinetics_ramped_temperature.ipynb b/examples/_eyring_kinetics_ramped_temperature.ipynb new file mode 100644 index 00000000..b5de45dc --- /dev/null +++ b/examples/_eyring_kinetics_ramped_temperature.ipynb @@ -0,0 +1,215 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sm\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from chempy import ReactionSystem\n", + "from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.rates import RampedTemp\n", + "sm.init_printing()\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "t, t0, A, B, C1 = sm.symbols('t t0 A B C1')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "y = C1/sm.E**((A*(((t + t0)*(-B + t + t0))/sm.E**(B/(t + t0)) - B**2*sm.Ei(-(B/(t + t0)))))/2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "y" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(y.diff(t)/y).simplify().expand().simplify().factor().powsimp(force=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "y.subs(t, 0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "yunit0 = y.subs(C1, C1/y.subs(t, 0)).simplify()\n", + "yunit0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.special import expi\n", + "f = sm.lambdify([t, t0, A, B], yunit0, modules=['numpy', {'Ei': expi}])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = 8.314472\n", + "T_K = 290\n", + "kB = 1.3806504e-23\n", + "h = 6.62606896e-34\n", + "dH = 80e3\n", + "dS = 10\n", + "rsys1 = ReactionSystem.from_string(\"\"\"\n", + "NOBr -> NO + Br; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n", + "\"\"\".format(dH=dH, dS=dS))\n", + "kref = 20836643994.118652*T_K*np.exp(-(dH - T_K*dS)/(R*T_K))\n", + "kref" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_A = kB/h*np.exp(dS/R)\n", + "_B = dH/R" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f(np.array([0, 1, 5, 20]), 290, _A, _B)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NOBr0_M = 0.7\n", + "init_cond = dict(\n", + " NOBr=NOBr0_M*u.M,\n", + " NO=0*u.M,\n", + " Br=0*u.M\n", + ")\n", + "t = 20*u.second" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(rsys):\n", + " odes, extra = get_odesys(rsys, unit_registry=si, constants=const, substitutions={\n", + " 'temperature': RampedTemp([T_K*u.K, 1*u.K/u.s])})\n", + " fig, all_axes = plt.subplots(2, 3, figsize=(14, 6))\n", + " for axes, odesys in zip(all_axes, [odes, odes.as_autonomous()]):\n", + " res = odesys.integrate(t, init_cond, integrator='cvode')\n", + " t_sec = to_unitless(res.xout, u.second)\n", + " NOBr_ref = NOBr0_M*f(t_sec, T_K, _A, _B)\n", + " cmp = to_unitless(res.yout, u.M)\n", + " ref = np.empty_like(cmp)\n", + " ref[:, odesys.names.index('NOBr')] = NOBr_ref\n", + " ref[:, odesys.names.index('Br')] = NOBr0_M - NOBr_ref\n", + " ref[:, odesys.names.index('NO')] = NOBr0_M - NOBr_ref\n", + " axes[0].plot(t_sec, cmp)\n", + " axes[1].plot(t_sec, cmp - ref)\n", + " res.plot_invariant_violations(ax=axes[2])\n", + " assert np.allclose(cmp, ref)\n", + " print({k: v for k, v in res.info.items() if not k.startswith('internal')}) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys2 = ReactionSystem.from_string(\"\"\"\n", + "NOBr -> NO + Br; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))\n", + "\"\"\".format(dH=dH, dS=dS))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys2)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_eyring_kinetics_ramped_temperature_2nd_order.ipynb b/examples/_eyring_kinetics_ramped_temperature_2nd_order.ipynb new file mode 100644 index 00000000..4fb37604 --- /dev/null +++ b/examples/_eyring_kinetics_ramped_temperature_2nd_order.ipynb @@ -0,0 +1,223 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sm\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "from chempy import ReactionSystem\n", + "from chempy.units import to_unitless, SI_base_registry as si, default_units as u, default_constants as const\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.rates import RampedTemp\n", + "sm.init_printing()\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "symbs = t, k, m, A, B, C1 = sm.symbols('t k m A B C1')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "y = -sm.E**(B/(m + k*t))*k/(\n", + " A*B*m - A*m**2 + A*B*k*t - 2*A*k*m*t - A*k**2*t**2 +\n", + " sm.E**(B/(m + k*t))*k*C1 +\n", + " A*B**2*sm.E**(B/(m + k*t))*sm.Ei(-(B/(m + k*t)))\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "(y.diff(t)/y).simplify().expand().simplify().factor().powsimp(force=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_C1, = sm.solve(y.subs(t, 0) - 1, C1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "yunit0 = y.subs(C1, _C1).simplify()\n", + "yunit0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(sm.python(yunit0))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from scipy.special import expi\n", + "f = sm.lambdify(symbs[:-1], yunit0, modules=['numpy', {'Ei': expi}])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "R = 8.314472\n", + "T_K = 290\n", + "dTdt_Ks = 3\n", + "kB = 1.3806504e-23\n", + "h = 6.62606896e-34\n", + "dH = 80e3\n", + "dS = 10\n", + "rsys1 = ReactionSystem.from_string(\"\"\"\n", + "2 NO2 -> N2O4; EyringParam(dH={dH}*J/mol, dS={dS}*J/K/mol)\n", + "\"\"\".format(dH=dH, dS=dS))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "_A = kB/h*np.exp(dS/R)\n", + "_B = dH/R" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "f(np.array([0, 1, 5, 20]), dTdt_Ks, T_K, _A, _B)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NO2_M = 1.0\n", + "init_cond = dict(\n", + " NO2=NO2_M*u.M,\n", + " N2O4=0*u.M\n", + ")\n", + "t = 20*u.second" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(rsys):\n", + " odes, extra = get_odesys(rsys, unit_registry=si, constants=const, substitutions={\n", + " 'temperature': RampedTemp([T_K*u.K, dTdt_Ks*u.K/u.s])})\n", + " fig, all_axes = plt.subplots(2, 3, figsize=(14, 6))\n", + " for axes, odesys in zip(all_axes, [odes, odes.as_autonomous()]):\n", + " res = odesys.integrate(t, init_cond, integrator='cvode')\n", + " t_sec = to_unitless(res.xout, u.second)\n", + " NO2_ref = f(t_sec, dTdt_Ks, T_K, _A, _B)\n", + " cmp = to_unitless(res.yout, u.M)\n", + " ref = np.empty_like(cmp)\n", + " ref[:, odesys.names.index('NO2')] = NO2_ref\n", + " ref[:, odesys.names.index('N2O4')] = (NO2_M - NO2_ref)/2\n", + " axes[0].plot(t_sec, cmp)\n", + " axes[1].plot(t_sec, cmp - ref)\n", + " res.plot_invariant_violations(ax=axes[2])\n", + " assert np.allclose(cmp, ref)\n", + " print({k: v for k, v in res.info.items() if not k.startswith('internal')}) " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys2 = ReactionSystem.from_string(\"\"\"\n", + "2 NO2 -> N2O4; MassAction(EyringHS([{dH}*J/mol, {dS}*J/K/mol]))\n", + "\"\"\".format(dH=dH, dS=dS))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(rsys2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_interactive_kinetic_modelling.ipynb b/examples/_interactive_kinetic_modelling.ipynb index 6375637b..bdfe9cde 100644 --- a/examples/_interactive_kinetic_modelling.ipynb +++ b/examples/_interactive_kinetic_modelling.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from collections import defaultdict\n", @@ -23,11 +19,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys = ReactionSystem.from_string(\"\"\"\n", @@ -40,11 +32,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys, extra = get_odesys(rsys, include_params=False, unit_registry=SI_base_registry,\n", @@ -54,11 +42,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "native = get_native(rsys, odesys, 'cvode')" @@ -67,11 +51,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot(\n", @@ -89,11 +69,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "interact(integrate_and_plot)" @@ -102,11 +78,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] } @@ -127,7 +99,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.1" }, "widgets": { "state": { diff --git a/examples/_kinetic_model_fitting.ipynb b/examples/_kinetic_model_fitting.ipynb index 236d2bb2..a1d1d1d5 100644 --- a/examples/_kinetic_model_fitting.ipynb +++ b/examples/_kinetic_model_fitting.ipynb @@ -2,7 +2,10 @@ "cells": [ { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "# Fitting kinetic parameters to experimental data\n", "Author: Björn Dahlgren.\n", @@ -19,9 +22,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "import bz2, codecs, collections, functools, itertools, json\n", @@ -42,7 +43,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "Experimental conditions, the two solutions which were mixed in 1:1 volume ratio in a stopped flow apparatus:" ] @@ -50,9 +54,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "sol1 = QuantityDict(u.mM, {'SCN-': 3*u.mM, 'K+': 3*u.mM, 'Na+': 33*u.mM, 'H+': 50*u.mM, 'ClO4-': (33+50)*u.mM})\n", @@ -62,9 +64,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "sol = (sol1 + sol2)/2 # 1:1 volume ratio at mixing\n", @@ -77,9 +77,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "ionic_strength_keys = 'abcd'\n", @@ -93,7 +91,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "We will read the data from a preprocessed file:" ] @@ -101,9 +102,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "transform = np.array([[1e-3, 0], [0, 1e-4]]) # converts 1st col: ms -> s and 2nd col to absorbance\n", @@ -118,7 +117,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "Let's plot the data:" ] @@ -126,9 +128,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "def mk_subplots(nrows=1, subplots_adjust=True, **kwargs):\n", @@ -149,9 +149,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def plot_series(series):\n", @@ -171,7 +169,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "We see that one data series is off: 16.5 ℃ and 0.0862 molal. Let's ignore that for now and perform the fitting, let's start with a pseudo-first order assumption (poor but simple):" ] @@ -179,9 +180,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def fit_pseudo1(serie, ax=None):\n", @@ -200,10 +199,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "def fit_all(series, fit_cb=fit_pseudo1, plot=False):\n", @@ -231,9 +227,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "def pseudo_to_k2(v):\n", @@ -249,9 +243,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "k2_unit = 1/u.M/u.s\n", @@ -273,7 +265,10 @@ }, { "cell_type": "markdown", - "metadata": {}, + "metadata": { + "deletable": true, + "editable": true + }, "source": [ "We define the reaction studied:" ] @@ -281,10 +276,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "eqsys = chempy.equilibria.EqSystem.from_string('Fe+3 + SCN- = FeSCN+2')\n", @@ -294,9 +286,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.kinetics.integrated import binary_rev" @@ -305,9 +295,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true - }, + "metadata": {}, "outputs": [], "source": [ "# TODO: show how to fit data to closed form expression, maybe also to non-linear kinetic model with complex-formation." diff --git a/examples/_kinetics_cstr.ipynb b/examples/_kinetics_cstr.ipynb index ee384a23..ceed66c6 100644 --- a/examples/_kinetics_cstr.ipynb +++ b/examples/_kinetics_cstr.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "%load_ext autoreload\n", @@ -17,11 +13,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from collections import defaultdict\n", @@ -39,11 +31,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys = ReactionSystem.from_string(\"A -> B; 'k'\", substance_factory=lambda k: Substance(k))\n", @@ -53,11 +41,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys, extra = get_odesys(rsys, include_params=False)\n", @@ -66,10 +50,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "We can change the expressions of the ODE system manually to account for source and sink terms from the flux:" ] @@ -77,11 +58,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "t, c1, c2, IA, IB, f, k, fc_A, fc_B = map(odesys.be.Symbol, 't c1 c2 I_A I_B f k phi_A phi_B'.split())\n", @@ -105,11 +82,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "cstr.param_names" @@ -118,11 +91,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "init_c, pars = {'A': .15, 'B': .1}, {'k': 0.8, 'f': 0.3, 'phi_A': .7, 'phi_B': .1}\n", @@ -132,11 +101,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res.plot()" @@ -145,9 +110,7 @@ { "cell_type": "markdown", "metadata": { - "collapsed": true, - "deletable": true, - "editable": true + "collapsed": true }, "source": [ "We can derive an analytic solution to the ODE system:" @@ -156,11 +119,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "k = cstr.params[0]\n", @@ -176,11 +135,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "exprs0 = [expr.subs(t, 0) for expr in exprs]\n", @@ -190,11 +145,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sol = cstr.be.solve([expr - c0 for expr, c0 in zip(exprs0, (IA, IB))], (c1, c2))\n", @@ -204,11 +155,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "exprs2 = [expr.subs(sol) for expr in exprs]\n", @@ -218,11 +165,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "IA" @@ -231,11 +174,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import sympy as sp\n", @@ -248,11 +187,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "exprs2_0 = [expr.subs(t, 0).simplify() for expr in exprs2]\n", @@ -262,11 +197,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "_cb = cstr.be.lambdify([t, IA, IB, k, f, fc_A, fc_B], exprs2)\n", @@ -277,11 +208,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def get_ref(result, parameters=None):\n", @@ -295,11 +222,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "yref = get_ref(res)\n", @@ -309,11 +232,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n", @@ -324,10 +243,7 @@ }, { "cell_type": "markdown", - "metadata": { - "deletable": true, - "editable": true - }, + "metadata": {}, "source": [ "Implementing CSTR expressions directly ``ReactionSystem.rates`` & ``get_odesys`` simplifies the code considerably:" ] @@ -335,11 +251,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "cstr2, extra2 = get_odesys(rsys, include_params=False, cstr=True)\n", @@ -349,11 +261,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "cstr2.param_names" @@ -362,11 +270,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "renamed_pars = {'k': pars['k'], 'fc_A': pars['phi_A'], 'fc_B': pars['phi_B'], 'feedratio': pars['f']}\n", @@ -378,11 +282,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.kinetics.integrated import unary_irrev_cstr\n", @@ -392,11 +292,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fr, fc = extra2['cstr_fr_fc']" @@ -405,11 +301,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def get_ref2(result):\n", @@ -427,11 +319,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res2.plot(y=res2.yout - get_ref2(res2))" @@ -440,11 +328,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "assert np.allclose(res2.yout, get_ref2(res2))" @@ -467,7 +351,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.5.2" + "version": "3.6.1" } }, "nbformat": 4, diff --git a/examples/_kinetics_ode_time_dependent.ipynb b/examples/_kinetics_ode_time_dependent.ipynb new file mode 100644 index 00000000..594e09b8 --- /dev/null +++ b/examples/_kinetics_ode_time_dependent.ipynb @@ -0,0 +1,142 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "from collections import defaultdict\n", + "from chempy import ReactionSystem\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.rates import SinTemp\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys = ReactionSystem.from_string(\"\"\"\n", + "2 HNO2 -> H2O + NO + NO2; EyringParam(dH=85e3, dS=10)\n", + "2 NO2 -> N2O4; EyringParam(dH=70e3, dS=20)\n", + "\"\"\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "st = SinTemp(unique_keys='Tbase Tamp Tangvel Tphase'.split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys, extra = get_odesys(rsys, include_params=False, substitutions={'temperature': st})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "init_conc = defaultdict(lambda: 0, HNO2=1, H2O=55)\n", + "params = dict(\n", + " Tbase=300,\n", + " Tamp=10,\n", + " Tangvel=2*math.pi/10,\n", + " Tphase=-math.pi/2\n", + ")\n", + "duration = 60" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(system):\n", + " result = system.integrate(duration, init_conc, params, integrator='cvode', nsteps=2000)\n", + " result.plot(names='NO HNO2 N2O4 NO2'.split())\n", + " print({k: v for k, v in sorted(result.info.items()) if not k.startswith('internal')})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(odesys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "autsys = odesys.as_autonomous()\n", + "assert len(autsys.dep) == len(odesys.dep) + 1" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "autsys.exprs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(autsys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_kinetics_ode_time_dependent_autonomous.ipynb b/examples/_kinetics_ode_time_dependent_autonomous.ipynb new file mode 100644 index 00000000..10e20ac5 --- /dev/null +++ b/examples/_kinetics_ode_time_dependent_autonomous.ipynb @@ -0,0 +1,274 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "from collections import defaultdict\n", + "import matplotlib.pyplot as plt\n", + "from chempy import ReactionSystem\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.rates import SinTemp\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys = ReactionSystem.from_string(\"\"\"\n", + "2 HNO2 -> H2O + NO + NO2; MassAction(EyringHS.fk('dH1', 'dS1'))\n", + "2 NO2 -> N2O4; MassAction(EyringHS.fk('dH2', 'dS2'))\n", + "\"\"\") # fictitious thermodynamic parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "st = SinTemp(unique_keys='Tbase Tamp Tangvel Tphase'.split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "class Constants:\n", + " Planck_constant = 6.63e-34\n", + " Boltzmann_constant = 1.38e-23\n", + " molar_gas_constant = 8.314\n", + "\n", + "odesys, extra = get_odesys(rsys, include_params=False, substitutions={'temperature': st}, constants=Constants())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "init_conc = defaultdict(lambda: 0, HNO2=1, H2O=55)\n", + "params = dict(\n", + " Tbase=300,\n", + " Tamp=10,\n", + " Tangvel=2*math.pi/10,\n", + " Tphase=-math.pi/2,\n", + " dH1=85e3,\n", + " dS1=10,\n", + " dH2=70e3,\n", + " dS2=20\n", + ")\n", + "duration = 60" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(system):\n", + " result = system.integrate(duration, init_conc, params, integrator='cvode', nsteps=2000)\n", + " fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n", + " result.plot(names='NO HNO2 N2O4'.split(), ax=axes[0])\n", + " result.plot(names='NO2'.split(), ax=axes[1])\n", + " print({k: v for k, v in sorted(result.info.items()) if not k.startswith('internal')})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(odesys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.param_names" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "len(odesys.exprs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys = odesys.as_autonomous()\n", + "len(asys.exprs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[a - o for a, o in zip(asys.exprs[:-1],odesys.exprs)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.exprs[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.get_jac()[:-1,:-1] - odesys.get_jac()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sym\n", + "sym.init_printing()\n", + "args = _x, _y, _p = asys.pre_process(*asys.to_arrays(1, init_conc, params))\n", + "args" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.f_cb(*args)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.j_cb(*args)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "argsode = odesys.pre_process(*odesys.to_arrays(1, init_conc, params))\n", + "argsode" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "argsode[0] - args[0], argsode[1] - args[1][:-1], argsode[2] - args[2]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.f_cb(*argsode)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.j_cb(*argsode)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(asys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.ny, asys.ny" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.pre_process(1, [0,1,2,3,4], [5,6,7,8])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_kinetics_ode_time_dependent_autonomous_units.ipynb b/examples/_kinetics_ode_time_dependent_autonomous_units.ipynb new file mode 100644 index 00000000..3e8c8d82 --- /dev/null +++ b/examples/_kinetics_ode_time_dependent_autonomous_units.ipynb @@ -0,0 +1,275 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "from collections import defaultdict\n", + "import matplotlib.pyplot as plt\n", + "from chempy import ReactionSystem\n", + "from chempy.units import (\n", + " default_constants,\n", + " default_units as u,\n", + " SI_base_registry as ureg\n", + ")\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.rates import SinTemp\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys = ReactionSystem.from_string(\"\"\"\n", + "2 HNO2 -> H2O + NO + NO2; MassAction(EyringHS.fk('dH1', 'dS1'))\n", + "2 NO2 -> N2O4; MassAction(EyringHS.fk('dH2', 'dS2'))\n", + "\"\"\") # fictitious thermodynamic parameters" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "st = SinTemp(unique_keys='Tbase Tamp Tangvel Tphase'.split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys, extra = get_odesys(rsys, include_params=False, substitutions={'temperature': st},\n", + " unit_registry=ureg, constants=default_constants)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "init_conc = defaultdict(lambda: 0*u.M, HNO2=1*u.M, H2O=55*u.M)\n", + "params = dict(\n", + " Tbase=300*u.K,\n", + " Tamp=10*u.K,\n", + " Tangvel=2*math.pi/(10*u.s),\n", + " Tphase=-math.pi/2,\n", + " dH1=85e3*u.J/u.mol,\n", + " dS1=10*u.J/u.K/u.mol,\n", + " dH2=70e3*u.J/u.mol,\n", + " dS2=20*u.J/u.K/u.mol\n", + ")\n", + "duration = 60*u.s" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(system):\n", + " result = system.integrate(duration, init_conc, params, integrator='cvode', nsteps=2000)\n", + " fig, axes = plt.subplots(1, 2, figsize=(14, 4))\n", + " result.plot(names='NO HNO2 N2O4'.split(), ax=axes[0])\n", + " result.plot(names='NO2'.split(), ax=axes[1])\n", + " print({k: v for k, v in sorted(result.info.items()) if not k.startswith('internal')})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(odesys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.param_names" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "len(odesys.exprs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys = odesys.as_autonomous()\n", + "len(asys.exprs)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "[a - o for a, o in zip(asys.exprs[:-1],odesys.exprs)]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.exprs[-1]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.get_jac()[:-1,:-1] - odesys.get_jac()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sympy as sym\n", + "sym.init_printing()\n", + "args = _x, _y, _p = asys.pre_process(*asys.to_arrays(1*u.s, init_conc, params))\n", + "args" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.f_cb(*args)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.j_cb(*args)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "argsode = odesys.pre_process(*odesys.to_arrays(1*u.s, init_conc, params))\n", + "argsode" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "argsode[0] - args[0], argsode[1] - args[1][:-1], argsode[2] - args[2]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.f_cb(*argsode)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.j_cb(*argsode)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(asys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys.ny, asys.ny" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "asys.pre_process(1, [0,1,2,3,4], [5,6,7,8])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_kinetics_ode_time_dependent_units.ipynb b/examples/_kinetics_ode_time_dependent_units.ipynb new file mode 100644 index 00000000..6320559f --- /dev/null +++ b/examples/_kinetics_ode_time_dependent_units.ipynb @@ -0,0 +1,146 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import math\n", + "from collections import defaultdict\n", + "from chempy import ReactionSystem\n", + "from chempy.units import (\n", + " default_units as u,\n", + " SI_base_registry as ureg\n", + ")\n", + "from chempy.kinetics.ode import get_odesys\n", + "from chempy.kinetics.rates import SinTemp\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "rsys = ReactionSystem.from_string(\"\"\"\n", + "2 HNO2 -> H2O + NO + NO2; EyringParam(dH=85e3*J/mol, dS=10*J/K/mol)\n", + "2 NO2 -> N2O4; EyringParam(dH=70e3*J/mol, dS=20*J/K/mol)\n", + "\"\"\"\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "st = SinTemp(unique_keys='Tbase Tamp Tangvel Tphase'.split())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "odesys, extra = get_odesys(rsys, include_params=False, substitutions={\n", + " 'temperature': st}, unit_registry=ureg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "init_conc = defaultdict(lambda: 0*u.M, HNO2=1*u.M, H2O=55*u.M)\n", + "params = dict(\n", + " Tbase=300*u.K,\n", + " Tamp=10*u.K,\n", + " Tangvel=2*math.pi/(10*u.s),\n", + " Tphase=-math.pi/2\n", + ")\n", + "duration = 60*u.s" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def integrate_and_plot(system):\n", + " result = system.integrate(duration, init_conc, params, integrator='cvode', nsteps=2000)\n", + " result.plot(names='NO HNO2 N2O4 NO2'.split())\n", + " print({k: v for k, v in sorted(result.info.items()) if not k.startswith('internal')})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(odesys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "autsys = odesys.as_autonomous()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "autsys.exprs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "integrate_and_plot(autsys)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.2" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/_ode_system_sympy_symbols.ipynb b/examples/_ode_system_sympy_symbols.ipynb index 590c5582..689c143a 100644 --- a/examples/_ode_system_sympy_symbols.ipynb +++ b/examples/_ode_system_sympy_symbols.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", @@ -20,11 +16,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "K1, K3, R1f, R3f = sp.symbols('K_1 K_3 R_1f R_3f')\n", @@ -44,11 +36,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys = mono_rev_single()" @@ -57,11 +45,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys" @@ -70,11 +54,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "type(rsys.rxns[1].param)" @@ -83,24 +63,16 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ - "odesys = get_odesys(rsys, include_params=False)[0]" + "odesys = get_odesys(rsys, include_params=False, params=True)[0]" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "assert len(odesys.params) == 4\n", @@ -110,11 +82,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys.get_jac()" @@ -123,11 +91,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "c0 = {'E': 5, 'A': 2, 'P': 0, 'X': 0}\n", @@ -137,11 +101,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys.f_cb(42, [c0[k] for k in rsys.substances], [params[k] for k in odesys.params])" @@ -150,11 +110,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys.names, odesys.param_names # get_odesys or SymbolicSys could determine names from Symbol.name" @@ -163,11 +119,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "result = odesys.integrate(42, c0, params, integrator='cvode')" @@ -176,29 +128,13 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "result.plot(xscale='log')\n", "plt.gca().set_xlim([1e-6, plt.gca().get_xlim()[1]])\n", "_ = plt.legend(loc='best')" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/examples/_regression.ipynb b/examples/_regression.ipynb index 33632095..09b200d6 100644 --- a/examples/_regression.ipynb +++ b/examples/_regression.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.units import default_units as u\n", @@ -18,11 +14,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "help(least_squares)" @@ -31,11 +23,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x = [0, 1, 2, 3, 4+1e-9]\n", @@ -55,11 +43,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res = least_squares(x, y)\n", @@ -79,11 +63,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_irls = irls(x, y, irls.gaussian, itermax=20)\n", @@ -93,11 +73,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "x" @@ -106,11 +82,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "y*u.m" @@ -119,11 +91,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_units = least_squares_units(x*u.s, y*u.m)\n", @@ -133,11 +101,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_units" @@ -146,11 +110,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "err = [.1, .2, .2, .8, 3]*u.m\n", @@ -161,11 +121,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] } diff --git a/examples/_switch_between_reduced_ode_systems.ipynb b/examples/_switch_between_reduced_ode_systems.ipynb index e81937a7..da094c4e 100644 --- a/examples/_switch_between_reduced_ode_systems.ipynb +++ b/examples/_switch_between_reduced_ode_systems.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -25,11 +21,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys = ReactionSystem.from_string(\n", @@ -45,11 +37,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsys.composition_balance_vectors()" @@ -58,11 +46,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "dep_scaling = 1e8\n", @@ -73,11 +57,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "orisys.linear_invariants, orisys.names, orisys.latex_names" @@ -86,11 +66,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "indep_end = 1e18\n", @@ -148,11 +124,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_results([orisys, scaledsys], integrate_systems([orisys, scaledsys], first_step=1e-23),\n", @@ -162,11 +134,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "psys = [PartiallySolvedSystem.from_linear_invariants(scaledsys, [k], description=k) for k in 'ABC']\n", @@ -187,11 +155,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "orisys['A']" @@ -200,11 +164,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rssys_A = PartiallySolvedSystem.from_linear_invariants(\n", @@ -215,11 +175,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_A = rssys_A.integrate(indep_end, init_conc, return_on_root=True, **integrate_kw)\n", @@ -229,11 +185,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n", @@ -262,11 +214,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "linresults = integrate_systems(scaled_linear)" @@ -275,11 +223,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_results(scaled_linear, linresults, xlim=(1e-9, indep_end))" @@ -288,11 +232,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_result('scaled A/C ', rssys_A.integrate(\n", @@ -313,11 +253,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "indep0 = 1e-26\n", @@ -331,11 +267,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "stlogresults = integrate_systems(stlogsystems, nsteps=500)" @@ -344,11 +276,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "plot_results(stlogsystems, stlogresults, vline_post_proc=lambda x: np.abs(np.exp(x) - indep0), xlim=[indep0*10, indep_end])" @@ -357,12 +285,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "stlogsystems[0].exprs" @@ -371,11 +294,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "unscaled_linear = [orisys] + [PartiallySolvedSystem.from_linear_invariants(\n", @@ -388,11 +307,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "utlogsystems[0].exprs" @@ -401,12 +316,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "LogLogSys_apart = symmetric_factory(None)\n", @@ -418,11 +328,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "ualogsystems = [LogLogSys_apart.from_other(ls, description='ual ' + ls.description) for ls in unscaled_linear]\n", @@ -433,11 +339,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "rsalogsys_A = LogLogSys_apart.from_other(rssys_A, autonomous_interface=True)\n", @@ -448,11 +350,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "res_log_switch = rsalogsys_A.integrate(indep_end, init_conc, return_on_root=True, **integrate_kw)\n", @@ -469,11 +367,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots(1, 1, figsize=(16, 4))\n", diff --git a/examples/ammonical_cupric_solution.ipynb b/examples/ammonical_cupric_solution.ipynb index 142050fb..231f2796 100644 --- a/examples/ammonical_cupric_solution.ipynb +++ b/examples/ammonical_cupric_solution.ipynb @@ -18,11 +18,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from collections import defaultdict\n", @@ -52,11 +48,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "NH3_complexes = ['CuNH3+2', 'Cu(NH3)2+2', 'Cu(NH3)3+2', 'Cu(NH3)4+2', 'Cu(NH3)5+2']\n", @@ -81,11 +73,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "show(', '.join([s.latex_name for s in substances]))" @@ -104,11 +92,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "init_conc = defaultdict(float, {'H+': 1e-7, 'OH-': 1e-7, 'NH4+': 0,\n", @@ -128,11 +112,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "H2O_c = init_conc['H2O']\n", @@ -163,11 +143,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "show(', '.join([s.latex_name for s in substances]))\n", @@ -196,11 +172,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "new_eqs = CuOH2_s - CuOH_B3, CuOH2_s - CuOH_B4\n", @@ -220,11 +192,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "#skip_subs, skip_eq = (4, 4) # (0, 0), (1, 1), (3, 3), (4, 4), (11, 9)\n", @@ -247,11 +215,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import sympy as sp\n", @@ -277,11 +241,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "numsys_log = NumSysLog(eqsys, backend=sp)\n", @@ -302,12 +262,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "sp.Matrix(1, len(f), lambda _, q: f[q]).jacobian(x)" @@ -326,11 +281,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "len(f), eqsys.ns" @@ -351,12 +302,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "C, sol, sane = eqsys.root(simpl_c0, NumSys=NumSysLog)\n", @@ -377,12 +323,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", @@ -409,11 +350,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sol_prod = Cout_logC[:, eqsys.as_substance_index('Cu+2')]*Cout_logC[:, eqsys.as_substance_index('OH-')]**2\n", @@ -442,11 +379,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(1, 2, figsize=(20, 8))\n", @@ -466,11 +399,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from scipy.optimize import newton\n", @@ -503,11 +432,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "numsys_log_rref = NumSysLog(eqsys, rref_equil=True, rref_preserv=True, backend=sp)\n", @@ -528,11 +453,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sp.Matrix(1, len(rf), lambda _, q: rf[q]).jacobian(x)" @@ -551,11 +472,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "plt.figure(figsize=(20,8))\n", @@ -577,12 +494,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "full_eqsys = EqSystem(equilibria[:-2] + new_eqs, substances)\n", @@ -595,11 +507,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def solve_and_plot_full(NumSys, plot_kwargs={'latex_names': True}, **kwargs):\n", @@ -632,11 +540,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "xout_log, sols_log, sane_log = solve_and_plot_full(NumSysLog, solver='scipy', tol=1e-10, conditional_maxiter=20,\n", @@ -649,11 +553,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "xout_log, sols_log, sane_log = solve_and_plot_full(NumSysLog, solver='scipy', tol=1e-12,\n", @@ -667,11 +567,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def sum_species(x, species, substance_names, weights=None):\n", @@ -707,12 +603,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "xout_static, sols_static, sane_static = solve_and_plot_full((NumSysLog), solver='scipy', tol=1e-12, # , NumSysSquare\n", @@ -729,11 +620,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# Not yet, see https://github.com/sympy/sympy/issues/10375\n", @@ -744,12 +631,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "#out_lin = solve_and_plot_full((NumSysLin,), x0=xout_log[0, :], ) # <-- this currently fails" @@ -758,12 +640,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "# out_log_lin = solve_and_plot_full((NumSysLog, NumSysLin)) # <-- this currently fails" diff --git a/examples/aqueous_radiolysis.ipynb b/examples/aqueous_radiolysis.ipynb index 6a112cd4..db2d50ae 100644 --- a/examples/aqueous_radiolysis.ipynb +++ b/examples/aqueous_radiolysis.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from __future__ import absolute_import, division, print_function\n", @@ -30,11 +26,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# Never mind the next row, it contains all the reaction data of aqueous radiolysis at room temperature\n", @@ -44,11 +36,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "species = [Substance.from_formula(s) for s in _species]\n", @@ -63,11 +51,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "# radiolytic yields for gamma radiolysis of neat water\n", @@ -86,12 +70,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": false - }, + "metadata": {}, "outputs": [], "source": [ "# The productions reactions have hardcoded rates corresponding to 1 Gy/s\n", @@ -102,11 +81,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "odesys, extra = get_odesys(rsys)\n", @@ -116,11 +91,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "j = odesys.get_jac()\n", @@ -130,11 +101,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot(odesys, c0_dict, integrator, ax=None, zero_conc=0, log10_t0=-12, tout=None,\n", @@ -156,11 +123,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "c0_dict = defaultdict(float, {'H+': 1e-7, 'OH-': 1e-7, 'H2O': 55.5})\n", @@ -173,11 +136,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def integrate_and_plot_scaled(rsys, dep_scaling, *args, **kwargs):\n", @@ -188,11 +147,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(2, 2, figsize=(14, 14))\n", @@ -209,12 +164,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true, - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ "integrate_and_plot(odesys, c0_dict, 'cvode', first_step=1e-12)" @@ -223,11 +173,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "integrate_and_plot(odesys, c0_dict, 'scipy')" @@ -236,11 +182,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from pyodesys.symbolic import symmetricsys\n", @@ -253,11 +195,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "loglogsys.exprs[:3]" @@ -266,11 +204,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "integrate_and_plot(loglogsys, c0_dict, 'gsl', zero_conc=1e-26, first_step=1e-3, log10_t0=-13)" @@ -279,11 +213,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "ssys, _ = get_odesys(rsys, SymbolicSys=ScaledSys, dep_scaling=1e16)" @@ -292,11 +222,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.kinetics._native import get_native\n", @@ -306,11 +232,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "integrate_and_plot(nsys, c0_dict, 'cvode', nsteps=96000,\n", @@ -320,11 +242,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "integrate_and_plot(nsys, c0_dict, 'cvode', nsteps=96000, tout=(1e-16, 1e6),\n", @@ -335,11 +253,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "print(''.join(open(next(filter(lambda s: s.endswith('.cpp'), nsys._native._written_files))).readlines()[:42]))" @@ -348,42 +262,20 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "init_conc_H2O2 = c0_dict.copy()\n", "init_conc_H2O2['H2O2'] = 0.01\n", - "odesys2 = odesys.__class__.from_other(odesys, lower_bounds=0)\n", - "res = integrate_and_plot(odesys2, init_conc_H2O2, integrator='cvode',\n", - " first_step=1e-16, atol=1e-12, rtol=1e-13, autorestart=5)" + "odesys2 = ScaledSys.from_other(odesys, lower_bounds=0, dep_scaling=1e8, indep_scaling=1e-6)\n", + "res = integrate_and_plot(odesys2, init_conc_H2O2, integrator='cvode', nsteps=500,\n", + " first_step=1e-16, atol=1e-5, rtol=1e-10, autorestart=1)" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, - "outputs": [], - "source": [ - "%load_ext autoreload\n", - "%autoreload 2\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from chempy.kinetics.analysis import plot_reaction_contributions" @@ -392,11 +284,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sks = ['H2O2', 'OH', 'HO2']\n", @@ -410,22 +298,14 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": true, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [] } diff --git a/examples/demo_debye_huckel.ipynb b/examples/demo_debye_huckel.ipynb index f9dc07f9..4bdfd4a2 100644 --- a/examples/demo_debye_huckel.ipynb +++ b/examples/demo_debye_huckel.ipynb @@ -3,11 +3,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "from math import log as ln\n", @@ -25,11 +21,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "def DH_A_B(T):\n", @@ -44,11 +36,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "DH_A_B(298.15*u.kelvin)" @@ -57,11 +45,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "T = (np.linspace(0, 40)+273.15)*u.Kelvin\n", @@ -81,11 +65,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "import sympy as sp\n", @@ -95,11 +75,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "constants = chempy.symbolic.get_constant_symbols()\n", @@ -109,11 +85,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "var_symbs = A, B, eps_r, T, rho, b0, I, z, a, I0, C, gamma = sp.symbols('A B epsilon_r T rho b^o I z a I^o C gamma')\n", @@ -124,11 +96,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "A_expr = dh.A(eps_r, T, rho, b0, constants, backend='sympy')\n", @@ -138,11 +106,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "B_expr = dh.B(eps_r, T, rho, b0, constants, backend=sp)\n", @@ -152,11 +116,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sp.Eq(sp.log(gamma), dh.limiting_log_gamma(I, z, A, I0=I0, backend=sp))" @@ -165,11 +125,7 @@ { "cell_type": "code", "execution_count": null, - "metadata": { - "collapsed": false, - "deletable": true, - "editable": true - }, + "metadata": {}, "outputs": [], "source": [ "sp.Eq(sp.log(gamma), dh.extended_log_gamma(I, z, a, A, B, C=C, I0=I0, backend=sp))" diff --git a/scripts/ci.sh b/scripts/ci.sh index ad48b492..3aecda22 100755 --- a/scripts/ci.sh +++ b/scripts/ci.sh @@ -5,9 +5,6 @@ if [[ "$CI_BRANCH" =~ ^v[0-9]+.[0-9]?* ]]; then echo ${CI_BRANCH} | tail -c +2 > __conda_version__.txt fi -python2 -m pip install git+https://github.com/bjodah/sympy@master -python3 -m pip install git+https://github.com/bjodah/sympy@master - git archive -o /tmp/$PKG_NAME.zip HEAD # test pip installable zip (symlinks break) python3 -m pip install /tmp/$PKG_NAME.zip diff --git a/scripts/rasterize.js b/scripts/rasterize.js deleted file mode 100644 index d25c2986..00000000 --- a/scripts/rasterize.js +++ /dev/null @@ -1,74 +0,0 @@ -// From: -// https://raw.githubusercontent.com/ariya/phantomjs/a86ae1948ca5ab20efc5c89ee39d183c33f59232/examples/rasterize.js -// -// Redistribution and use in source and binary forms, with or without -// modification, are permitted provided that the following conditions are met: - -// * Redistributions of source code must retain the above copyright -// notice, this list of conditions and the following disclaimer. -// * Redistributions in binary form must reproduce the above copyright -// notice, this list of conditions and the following disclaimer in the -// documentation and/or other materials provided with the distribution. -// * Neither the name of the nor the -// names of its contributors may be used to endorse or promote products -// derived from this software without specific prior written permission. - -// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" -// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE -// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE -// ARE DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY -// DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES -// (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; -// LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND -// ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT -// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF -// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - -var page = require('webpage').create(), - system = require('system'), - address, output, size; - -if (system.args.length < 3 || system.args.length > 5) { - console.log('Usage: rasterize.js URL filename [paperwidth*paperheight|paperformat] [zoom]'); - console.log(' paper (pdf output) examples: "5in*7.5in", "10cm*20cm", "A4", "Letter"'); - console.log(' image (png/jpg output) examples: "1920px" entire page, window width 1920px'); - console.log(' "800px*600px" window, clipped to 800x600'); - phantom.exit(1); -} else { - address = system.args[1]; - output = system.args[2]; - page.viewportSize = { width: 600, height: 600 }; - if (system.args.length > 3 && system.args[2].substr(-4) === ".pdf") { - size = system.args[3].split('*'); - page.paperSize = size.length === 2 ? { width: size[0], height: size[1], margin: '0px' } - : { format: system.args[3], orientation: 'portrait', margin: '1cm' }; - } else if (system.args.length > 3 && system.args[3].substr(-2) === "px") { - size = system.args[3].split('*'); - if (size.length === 2) { - pageWidth = parseInt(size[0], 10); - pageHeight = parseInt(size[1], 10); - page.viewportSize = { width: pageWidth, height: pageHeight }; - page.clipRect = { top: 0, left: 0, width: pageWidth, height: pageHeight }; - } else { - console.log("size:", system.args[3]); - pageWidth = parseInt(system.args[3], 10); - pageHeight = parseInt(pageWidth * 3/4, 10); // it's as good an assumption as any - console.log ("pageHeight:",pageHeight); - page.viewportSize = { width: pageWidth, height: pageHeight }; - } - } - if (system.args.length > 4) { - page.zoomFactor = system.args[4]; - } - page.open(address, function (status) { - if (status !== 'success') { - console.log('Unable to load the address!'); - phantom.exit(1); - } else { - window.setTimeout(function () { - page.render(output); - phantom.exit(); - }, 200); - } - }); -} diff --git a/scripts/rasterize.js.bz2 b/scripts/rasterize.js.bz2 new file mode 100644 index 00000000..7e002225 Binary files /dev/null and b/scripts/rasterize.js.bz2 differ diff --git a/scripts/render_index.sh b/scripts/render_index.sh index 36615a9f..c5495b76 100755 --- a/scripts/render_index.sh +++ b/scripts/render_index.sh @@ -7,6 +7,8 @@ mkdir -p thumbs tmpdir=$(mktemp -d) trap "rm -r $tmpdir" INT TERM EXIT +SCRIPTS_PATH=$(unset CDPATH && cd "$(dirname "$0")" && echo $PWD) +bzcat $SCRIPTS_PATH/rasterize.js.bz2 >$SCRIPTS_PATH/rasterize.js cat <index.html @@ -23,7 +25,7 @@ for f in $@; do continue # don't include bokeh apps fi img=$(basename $f .html).png - QT_QPA_PLATFORM=offscreen phantomjs $(unset CDPATH && cd "$(dirname "$0")" && echo $PWD)/rasterize.js $f $tmpdir/$img 1200px*900px + phantomjs $SCRIPTS_PATH/rasterize.js $f $tmpdir/$img 1200px*900px convert $tmpdir/$img -resize 400x300 thumbs/$img cat <>index.html

diff --git a/setup.py b/setup.py index c94616d7..4608ff65 100644 --- a/setup.py +++ b/setup.py @@ -97,7 +97,7 @@ def _path_under_setup(*args): _author, _author_email = open(_path_under_setup('AUTHORS'), 'rt').readline().split('<') extras_req = { - 'integrators': ['scipy>=0.16.1', 'pyodeint>=0.7.0', 'pycvodes>=0.6.1', 'pygslodeiv2>=0.6.1'], + 'integrators': ['scipy>=0.16.1', 'pyodeint>=0.7.0', 'pycvodes>=0.8.3', 'pygslodeiv2>=0.6.1'], 'solvers': ['pykinsol'], 'native': ['pycompilation>=0.4.3', 'pycodeexport>=0.1.1', 'appdirs'], 'docs': ['Sphinx', 'sphinx_rtd_theme', 'numpydoc'], @@ -121,7 +121,7 @@ def _path_under_setup(*args): install_requires=[ 'numpy>1.7', 'scipy>=0.16.1', 'matplotlib>=1.3.1', 'sympy>=1.0', 'quantities>=0.11.1', 'pyneqsys>=0.4.4', - 'pyodesys>=0.7.0', 'pyparsing>=2.0.3', 'sym' + 'pyodesys>=0.10.3', 'pyparsing>=2.0.3', 'sym' # 'dot2tex>=2.9.0' ], extras_require=extras_req