Merge pull request #35 from bjodah/more-features

More features

CHANGES.rst

@@ -1,3 +1,10 @@
+v0.4.1
+======
+- Fixes for enhanced robustness:
+ - ``.kinetics.ode.get_odesys``
+ - ``.chemistry.as_per_substance_array``
+- Minor changes.
+
 v0.4.0
 ======
 - Multiple fixes throughout

chempy/chemistry.py

@@ -1301,7 +1301,7 @@ def as_per_substance_array(self, cont, dtype='float64', unit=None, raise_on_unk=
  raise KeyError("Unkown substance key: %s" % k)
  cont = [cont[k] for k in substance_keys]
 
- cont = np.asarray(cont, dtype=dtype)
+ cont = np.atleast_1d(np.asarray(cont, dtype=dtype).squeeze())
  if cont.shape[-1] != self.ns:
  raise ValueError("Incorrect size")
  return cont*(unit if unit is not None else 1)

chempy/kinetics/_rates.py

@@ -16,7 +16,7 @@ class TPolyMassAction(TPoly, MassAction):
  argument_names = ('temperature_offset', Ellipsis)
  parameter_keys = TPoly.parameter_keys
 
- def rate_coeff(self, variables, backend):
+ def rate_coeff(self, variables, backend=math):
  return self.eval_poly(variables, backend)
 
 
@@ -27,7 +27,7 @@ def rate_coeff(self, variables, backend):
 class PiecewiseTPolyMassAction(TPiecewisePoly, MassAction):
  parameter_keys = TPoly.parameter_keys
 
- def rate_coeff(self, variables, backend):
+ def rate_coeff(self, variables, backend=math):
  return self.eval_poly(variables, backend)
 
 
@@ -39,19 +39,27 @@ def g_value(self, variables, backend):
  return self.eval_poly(variables, backend)
 
 
+class RTPolyRadiolytic(RTPoly, Radiolytic):
+ nargs = None
+ parameter_keys = Radiolytic.parameter_keys + RTPoly.parameter_keys
+
+ def g_value(self, variables, backend):
+ return self.eval_poly(variables, backend)
+
+
 class RTPolyMassAction(RTPoly, MassAction):
  """ Arguments: temperature_offset, c0, c1, ... """
  parameter_keys = RTPoly.parameter_keys
  nargs = None
 
- def rate_coeff(self, variables, backend):
+ def rate_coeff(self, variables, backend=math):
  return self.eval_poly(variables, backend)
 
 
 class _Log10XPolyMassAction(MassAction):
  skip_poly = 1 # kunit
 
- def rate_coeff(self, variables, backend):
+ def rate_coeff(self, variables, backend=math):
  k_unit = self.arg(variables, 0)
  return 10**self.eval_poly(variables, backend)*k_unit
 

chempy/kinetics/ode.py

@@ -162,7 +162,6 @@ def get_odesys(rsys, include_params=True, substitutions=None,
  if ratex.unique_keys is not None:
  unique_keys.extend(ratex.unique_keys)
  p_defaults.extend(ratex.args)
-
  if unit_registry is None:
  def pre_processor(x, y, p):
  return (
@@ -184,9 +183,9 @@ def post_processor(x, y, p):
  p_units = [get_derived_unit(unit_registry, k) for k in param_keys]
  new_r_exprs = []
  for ratex in r_exprs:
- _pu, _new_rates = ratex.dedimensionalisation(unit_registry)
+ _pu, _new_rate = ratex._recursive_as_RateExpr().dedimensionalisation(unit_registry)
  p_units.extend(_pu)
- new_r_exprs.append(_new_rates)
+ new_r_exprs.append(_new_rate)
  r_exprs = new_r_exprs
 
  time_unit = get_derived_unit(unit_registry, 'time')

chempy/kinetics/rates.py

@@ -31,6 +31,22 @@ def rxn(self, value):
  if isinstance(arg, RateExpr):
  arg.rxn = value
 
+ def _recursive_as_RateExpr(self):
+ new_args = []
+ for arg in self.args:
+ if isinstance(arg, Expr):
+ new_args.append(arg)
+ else:
+ if hasattr(arg, '_as_RateExpr'):
+ new_args.append(arg._as_RateExpr(self.rxn))
+ else:
+ new_args.append(arg)
+ if self.kw is None:
+ kw = {}
+ else:
+ kw = {k: getattr(self, k) for k in self.kw}
+ return self.__class__(new_args, self.unique_keys, **kw)
+
  @classmethod
  @deprecated(use_instead=Expr.from_callback)
  def subclass_from_callback(cls, cb, cls_attrs=None):
@@ -91,7 +107,7 @@ class _Radiolytic(RadiolyticBase):
  parameter_keys = (doserate_name, 'density')
  print_name = 'Radiolytic' if doserate_name == 'doserate' else ('Radiolytic{'+doserate_name+'}')
 
- def g_value(self, variables, backend): # for subclasses
+ def g_value(self, variables, backend=math): # for subclasses
  return self.arg(variables, 0, backend=backend)
 
  def __call__(self, variables, backend=math):
@@ -168,6 +184,22 @@ def as_mass_action(self, variables, backend=math):
 
 
 class ArrheniusMassAction(MassAction):
+ """ Rate expression for a Arrhenius-type of rate
+
+ Examples
+ --------
+ >>> from math import exp
+ >>> from chempy import Reaction
+ >>> from chempy.units import allclose, default_units as u
+ >>> A = 1e11 / u.second
+ >>> Ea_over_R = 42e3/8.3145 * u.K**-1
+ >>> ratex = ArrheniusMassAction([A, Ea_over_R])
+ >>> rxn = Reaction({'R'}, {'P'}, ratex)
+ >>> dRdt = rxn.rate({'R': 3*u.M, 'temperature': 298.15*u.K})['R']
+ >>> allclose(dRdt, -3*1e11*exp(-42e3/8.3145/298.15)*u.M/u.s)
+ True
+
+ """
  argument_names = ('A', 'Ea_over_R')
  parameter_keys = ('temperature',)
 

chempy/tests/test_chemistry.py

@@ -187,8 +187,7 @@ def test_ReactionSystem__as_per_substance_array_dict():
  m = default_units.metre
  M = default_units.molar
  rs = ReactionSystem([], [Substance('H2O')])
- c = rs.as_per_substance_array({'H2O': 1*M},
- unit=M)
+ c = rs.as_per_substance_array({'H2O': 1*M}, unit=M)
  assert c.dimensionality == M.dimensionality
  assert abs(c[0]/(1000*mol/m**3) - 1) < 1e-16
 

chempy/tests/test_units.py

@@ -11,9 +11,9 @@
 
 from ..util.testing import requires
 from ..units import (
- allclose, get_derived_unit, is_unitless, linspace,
+ allclose, get_derived_unit, is_unitless, linspace, logspace_from_lin,
  SI_base_registry, unitless_in_registry, format_string, get_physical_quantity,
- to_unitless, magnitude, default_unit_in_registry, Backend,
+ to_unitless, magnitude, default_unit_in_registry, Backend, latex_of_unit,
  unit_of, unit_registry_to_human_readable, units_library,
  unit_registry_from_human_readable, _sum, UncertainQuantity,
  default_units as u
@@ -104,6 +104,7 @@ def test_to_unitless():
  assert abs(to_unitless(3/(u.second*u.molar),
  u.metre**3/u.mole/u.second) - 3e-3) < 1e-12
  assert abs(to_unitless(2*u.dm3, u.cm3) - 2000) < 1e-12
+ assert abs(to_unitless(2*u.m3, u.dm3) - 2000) < 1e-12
  assert (float(to_unitless(UncertainQuantity(2, u.dm3, .3), u.cm3)) - 2000) < 1e-12
 
  g1 = UncertainQuantity(4.46, u.per100eV, 0)
@@ -137,6 +138,13 @@ def test_linspace():
  assert abs(to_unitless(ls[0], u.hour) - 2/3600.) < 1e-15
 
 
+@requires(units_library)
+def test_logspace_from_lin():
+ ls = logspace_from_lin(2*u.second, 3*u.second)
+ assert abs(to_unitless(ls[0], u.hour) - 2/3600.) < 1e-15
+ assert abs(to_unitless(ls[-1], u.hour) - 3/3600.) < 1e-15
+
+
 @requires(units_library)
 def test_get_derived_unit():
  registry = SI_base_registry.copy()
@@ -323,3 +331,8 @@ def test_joule_html():
  joule_htm = 'kg&sdot;m<sup>2</sup>/s<sup>2</sup>'
  joule = u.J.dimensionality.simplified
  assert joule.html == joule_htm
+
+
+@requires(units_library)
+def test_latex_of_unit():
+ assert latex_of_unit(u.gram/u.metre**2) == r'\mathrm{\frac{g}{m^{2}}}'

chempy/units.py

@@ -5,19 +5,21 @@
 - ``chempy.units.default_constants``
 - ``chempy.units.SI_base_registry``
 
+together with some functions.
+
 Currently `quantities <https://pypi.python.org/pypi/quantities>`_ is used as
 the underlying package to handle units. If it is possible you should try to
-only use the `chempy.units` module (in case ``ChemPy`` changes this backend).
+only use the `chempy.units` module (in case ``ChemPy`` changes this backend),
+and avoid relying on any attributes of the Quantity instances (and rather use
+functions in `chempy.units`).
 
 """
 from __future__ import (absolute_import, division, print_function)
 
-from operator import mul
 from functools import reduce
+from operator import mul
 
 from ._util import NameSpace
-# Currently we use quantities for units. This may change, therefore use this
-# file for all units. A requirement is first-class numpy support.
 
 units_library = 'quantities' # info used for selective testing.
 
@@ -40,6 +42,7 @@
  default_units = NameSpace(pq)
  default_units.decimetre = pq.UnitQuantity(
  'decimetre', default_units.m / 10.0, u_symbol='dm')
+ default_units.m3 = default_units.metre**3
  default_units.dm3 = default_units.decimetre**3
  default_units.cm3 = default_units.centimetre**3
  if not hasattr(default_units, 'molar'):
@@ -147,6 +150,18 @@ def unit_registry_to_human_readable(unit_registry):
  return new_registry
 
 
+def latex_of_unit(quant):
+ """ Returns LaTeX reperesentation of the unit of a quantity
+
+ Examples
+ --------
+ >>> print(latex_of_unit(1/default_units.kelvin))
+ \\mathrm{\\frac{1}{K}}
+
+ """
+ return quant.dimensionality.latex.strip('$')
+
+
 def unit_registry_from_human_readable(unit_registry):
  """ Deserialization of unit_registry. """
  if unit_registry is None:
@@ -308,7 +323,15 @@ def allclose(a, b, rtol=1e-8, atol=None):
 
 
 def linspace(start, stop, num=50):
- """ Analogous to ``numpy.linspace``. """
+ """ Analogous to ``numpy.linspace``.
+
+ Examples
+ --------
+ >>> abs(linspace(2, 8, num=3)[1] - 5) < 1e-15
+ True
+
+ """
+
  # work around for quantities v0.10.1 and NumPy
  import numpy as np
  unit = unit_of(start)
@@ -317,6 +340,22 @@ def linspace(start, stop, num=50):
  return np.linspace(start_, stop_, num)*unit
 
 
+def logspace_from_lin(start, stop, num=50):
+ """ Logarithmically spaced data points
+
+ Examples
+ --------
+ >>> abs(logspace_from_lin(2, 8, num=3)[1] - 4) < 1e-15
+ True
+
+ """
+ import numpy as np
+ unit = unit_of(start)
+ start_ = np.log2(to_unitless(start, unit))
+ stop_ = np.log2(to_unitless(stop, unit))
+ return np.exp2(np.linspace(start_, stop_, num))*unit
+
+
 def _sum(iterable):
  try:
  result = next(iterable)