Merge pull request #52 from bjodah/bump-ci-image

Use a tagged dockerimage for CI

.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:

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

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)]

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

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)

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)

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,27 +205,83 @@ 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):
  """ Ramped temperature, pass as substitution to e.g. ``get_odesys`` """
  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)