How to use the chempy.units.default_units function in chempy

from pyodesys.results import Result
    if isinstance(xyp, Result):
        xyp = xyp.odesys.to_arrays(xyp.xout, xyp.yout, xyp.params, reshape=False)
    if varied is None:
        varied = xyp[0]
    if xyp[1].shape[-2] != varied.size:
        raise ValueError("Size mismatch between varied and yout")
    if substance_keys is None:
        substance_keys = rsys.substances.keys()
    if axes is None:
        _fig, axes = plt.subplots(len(substance_keys))
    rates = rate_exprs_cb(*xyp)
    if unit_registry is not None:
        time_unit = get_derived_unit(unit_registry, 'time')
        conc_unit = get_derived_unit(unit_registry, 'concentration')
        rates = to_unitless(rates*conc_unit/time_unit, u.molar/u.second)

    eqk1, eqk2, eqs = _combine_rxns_to_eq(rsys) if combine_equilibria else ([], [], [])

    for sk, ax in zip(substance_keys, axes):
        data, tot = _dominant_reaction_effects(sk, rsys, rates, linthreshy, eqk1, eqk2, eqs)
        factor = 1/xyp[1][:, rsys.as_substance_index(sk)] if relative else 1
        if total:
            ax.plot(varied, factor*tot, c='k', label='Total', linewidth=2, ls=':')
        for y, rxn in sorted(data, key=lambda args: args[0][-1], reverse=True):
            ax.plot(varied, factor*y,
                    label=r'$\mathrm{%s}$' % rxn.latex(rsys.substances))

        if rsys.substances[sk].latex_name is None:
            ttl = rsys.substances[sk].name
            ttl_template = '%s'
        else:

try:
        import numpy
    except ImportError:
        def _numpy_not_installed_raise(*args, **kwargs):
            raise ImportError("numpy not installed, no such method")

        class numpy:
            array = staticmethod(_numpy_not_installed_raise)
            log = staticmethod(_numpy_not_installed_raise)
            exp = staticmethod(_numpy_not_installed_raise)

    _update(numpy, keys='array log exp'.split())  # could of course add more
    _update(rates)
    _update(chempy)
    for df in [default_units, default_constants]:
        if df is not None:
            globals_.update(df.as_dict())
    return globals_

if backend.__name__ != 'sympy':
            warnings.warn("Backend != SymPy, provide your own transform function.")

        def transform(arg):
            expr = backend.logcombine(arg, force=True)
            v, w = map(backend.Wild, 'v w'.split())
            expr = expr.replace(backend.log(w**v), v*backend.log(w))
            return expr

    args = [symbols[key] for key in conditions]
    seen = [False]*len(args)
    rates = {}
    for k, v in rsys.rates(symbols).items():
        expr = transform(v)
        if expr == 0:
            rate = 0 * u.molar/u.second
        else:
            rate = backend.lambdify(args, expr)(*conditions.values())
            to_unitless(rate, u.molar/u.second)
        rates[k] = rate
        seen = [b or a in expr.free_symbols for b, a in zip(seen, args)]
    not_seen = [a for s, a in zip(seen, args) if not s]
    for k in conditions:
        if k not in odesys.param_names and k not in odesys.names and k not in ignore:
            raise KeyError("Unknown param: %s" % k)
    return {'not_seen': not_seen, 'rates': rates}

Start by runing:
    $ bokeh serve interactive.py
Add --show argument or navigate to:
    http://localhost:5006/interactive
"""
from collections import defaultdict

from chempy.util.bkh import integration_with_sliders
from chempy.units import SI_base_registry, default_units as u

from bokeh_interactive_arrhenius import get_rsys


if __name__.startswith('bk_'):
    from bokeh.io import curdoc
    Af, Ab, Ea, Er = 1e16/u.molar/u.s, 1.5e15/u.s, 72e3*u.J/u.mol, -12e3*u.J/u.mol
    curdoc().add_root(integration_with_sliders(
        get_rsys(Af, Ab, Ea, Er), tend=3*u.s,
        c0=defaultdict(lambda: 0*u.molar, {'Fe+3': 3e-3*u.molar, 'SCN-': 1.5e-3*u.molar}),
        parameters={'temperature': 298.15*u.K},
        slider_kwargs={'temperature': dict(start=273.15*u.K, end=313.15*u.K, step=.05*u.K)},
        unit_registry=SI_base_registry,
        output_conc_unit=u.molar,
        output_time_unit=u.second
    ))
else:
    import warnings
    warnings.warn("Run using 'bokeh serve %s'" % __file__)

    def __call__(self, T, units=default_units, backend=None):
        """ Evaluates Henry's constant for provided temperature """
        return super(HenryWithUnits, self).__call__(T, units, backend)

def from_rateconst_at_T(cls, *args, **kwargs):
        if 'constants' not in kwargs:
            kwargs['constants'] = default_constants
        if 'units' not in kwargs:
            kwargs['units'] = default_units
        if 'backend' not in kwargs:
            kwargs['backend'] = patched_numpy
        return super(ArrheniusParamWithUnits, cls).from_rateconst_at_T(*args, **kwargs)

r_str = ('{}, 0 \u21D2 {}' if ratcoeff > 0 else '{}, {} \u21D2 0').format(
            parmap[pk], substmap[prod])
    else:
        raise NotImplementedError("Whats that?")
    return r_str, pk, abs(ratcoeff)


class DiffEqBioJl:
    _template_body = """\
{name} = @{crn_macro} begin
    {reactions}
end {parameters}
{post}
"""

    defaults = dict(unit_conc=u.molar, unit_time=u.second)

    def __init__(self, *, rxs, pars, substance_key_map, parmap, **kwargs):
        self.rxs = rxs
        self.pars = pars
        self.substance_key_map = substance_key_map
        self.parmap = parmap
        self.unit_conc = kwargs.get('unit_conc', self.defaults['unit_conc'])
        self.unit_time = kwargs.get('unit_time', self.defaults['unit_time'])
        self.latex_names = kwargs.get('latex_names', None)

    @classmethod
    def from_rsystem(cls, rsys, par_vals, *, variables=None, substance_key_map=lambda i, sk: 'y%d' % i, **kwargs):
        if not isinstance(substance_key_map, dict):
            substance_key_map = {sk: substance_key_map(si, sk) for si, sk in enumerate(rsys.substances)}
        parmap = dict([(r.param.unique_keys[0], 'p%d' % i) for i, r in enumerate(rsys.rxns)])
        rxs, pars = [], OrderedDict()

def check_consistent_units(self, throw=False):
        if is_quantity(self.param):  # This will assume mass action
            exponent = sum(self.prod.values()) - sum(self.reac.values())
            unit_param = unit_of(self.param, simplified=True)
            unit_expected = unit_of(default_units.molar**exponent, simplified=True)
            if unit_param == unit_expected:
                return True
            else:
                if throw:
                    raise ValueError("Inconsistent units in equilibrium: %s vs %s" %
                                     (unit_param, unit_expected))
                else:
                    return False
        else:
            return True  # the user might not be using ``chempy.units``

    def __call__(self, state, constants=default_constants, units=default_units,
                 backend=None):
        """ See :func:`chempy.eyring.eyring_equation`. """
        if backend is None:
            backend = Backend()
        return super(EyringParamWithUnits, self).__call__(
            state, constants, units, backend)

def molar_mass(self, units=None):
        """ Returns the molar mass (with units) of the substance

        Examples
        --------
        >>> nh4p = Substance.from_formula('NH4+')  # simpler
        >>> from chempy.units import default_units as u
        >>> nh4p.molar_mass(u)
        array(18.0384511...) * g/mol

        """
        if units is None:
            units = default_units
        return self.mass*units.g/units.mol