How to use the chempy._util.get_backend function in chempy

def stoichs_constants(self, eq_params=None, rref=False, Matrix=None,
                          backend=None, non_precip_rids=()):
        if eq_params is None:
            eq_params = self.eq_constants()
        if rref:
            from pyneqsys.symbolic import linear_rref
            be = get_backend(backend)
            rA, rb = linear_rref(self.stoichs(non_precip_rids),
                                 list(map(be.log, eq_params)),
                                 Matrix)
            return rA.tolist(), list(map(be.exp, rb))
        else:
            return (self.stoichs(non_precip_rids),
                    eq_params)

def davies_activity_product(I, stoich, z, a, T, eps_r, rho, C=-0.3, backend=None):
    be = get_backend(backend)
    Aval = A(eps_r, T, rho)
    tot = 0
    for idx, nr in enumerate(stoich):
        tot += nr*davies_log_gamma(I, z[idx], Aval, C)
    return be.exp(tot)

def davies_log_gamma(I, z, A, C=-0.3, I0=1, backend=None):
    """ Davies formula """
    be = get_backend(backend)
    one = be.pi**0
    I_I0 = I/I0
    sqrt_I_I0 = (I_I0)**(one/2)
    return -A * z**2 * (sqrt_I_I0/(1 + sqrt_I_I0) + C*I_I0)

else:
            see source code for what attributes are used.
            Tip: pass quantities.constants

    Notes
    -----
    Remember to divide by ln(10) if you want to use the constant
    with log10 based expression.

    References
    ----------
    Atkins, De Paula, Physical Chemistry, 8th edition

    """
    b0 = _get_b0(b0, units)
    be = get_backend(backend)
    one = be.pi**0
    if constants is None:
        combined = 132871.85866393594
        if units is not None:
            m = units.meter
            K = units.Kelvin
            mol = units.mol
            combined *= (m*K)**(3*one/2) / mol**(one/2)
        return combined*(rho * b0 * T**-3 * eps_r**-3)**0.5
    F = constants.Faraday_constant
    NA = constants.Avogadro_constant
    eps0 = constants.vacuum_permittivity
    kB = constants.Boltzmann_constant
    pi = constants.pi
    A = F**3/(4*pi*NA)*(rho*b0/(2*(eps0*eps_r*kB*NA*T)**3))**(one/2)
    return A

Default is 'numpy', can also be e.g. ``sympy``.

    Returns
    -------
    length-2 tuple
        concentrations of reactant and product

    """
    # Mathematica source:
    # FortranForm[
    # DSolve[{Derivative[1][A][t] == a0 f - f A[t] - 2 k A[t]^2,
    #   Derivative[1][B][t] == b0 f + n*k A[t]^2 - f B[t], A[0] == x,
    #   B[0] == y}, {A[t], B[t]}, {t}]]
    # Post processed using sympy's cse function.
    # (see _derive_analytic_cstr_bireac.ipynb)
    be = get_backend(backend)
    atanh = getattr(be, 'atanh', be.arctanh)
    three = 3*be.cos(0)

    x0 = 1/k
    x1 = be.sqrt(fv)
    x2 = 8*k
    x3 = fr*x2
    x4 = be.sqrt(fv + x3)
    x5 = x1*x4
    x6 = x1*x4/2
    x7 = atanh((-fv**(three/2)*x4 - 4*k*r*x5)/(fv**2 + fv*x3))
    x8 = fv*t
    x9 = fp*x2
    x10 = 4*k*n
    x11 = fr*x10
    x12 = be.exp(x8)

kf : number or Symbol
        Forward (bimolecular) rate constant.
    kb : number or Symbol
        Backward (unimolecular) rate constant.
    prod : number or Symbol
        Initial concentration of the complex.
    major : number or Symbol
        Initial concentration of the more abundant reactant.
    minor : number or Symbol
        Initial concentration of the less abundant reactant.
    backend : module or str
        Default is 'numpy', can also be e.g. ``sympy``.

    """
    # see _integrated.ipynb for derivation
    be = get_backend(backend)
    X, Y, Z = prod, major, minor
    x0 = Y*kf
    x1 = Z*kf
    x2 = 2*X*kf
    x3 = -kb - x0 - x1
    x4 = -x2 + x3
    x5 = be.sqrt(-4*kf*(X**2*kf + X*x0 + X*x1 + Z*x0) + x4**2)
    x6 = kb + x0 + x1 + x5
    x7 = (x3 + x5)*be.exp(-t*x5)
    x8 = x3 - x5
    return (x4*x8 + x5*x8 + x7*(x2 + x6))/(2*kf*(x6 + x7))

Entropy of activation.
    T: float with unit
        temperature
    constants: object (optional, default: None)
        if None:
            T assumed to be in Kelvin, Ea in J/(K mol)
        else:
            attributes accessed: molar_gas_constant
            Tip: pass quantities.constants
    units: object (optional, default: None)
        attributes accessed: Joule, Kelvin and mol
    backend: module (optional)
        module with "exp", default: numpy, math

    """
    be = get_backend(backend)
    R = _get_R(constants, units)
    kB_over_h = _get_kB_over_h(constants, units)

    try:
        RT = (R*T).rescale(dH.dimensionality)
    except AttributeError:
        RT = R*T

    try:
        kB_over_h = kB_over_h.simplified
    except AttributeError:
        pass

    return kB_over_h*T*be.exp(dS/R)*be.exp(-dH/RT)

def limiting_activity_product(I, stoich, z, T, eps_r, rho, backend=None):
    """ Product of activity coefficients based on DH limiting law. """
    be = get_backend(backend)
    Aval = A(eps_r, T, rho)
    tot = 0
    for idx, nr in enumerate(stoich):
        tot += nr*limiting_log_gamma(I, z[idx], Aval)
    return be.exp(tot)

modules which contains "exp", default: numpy, math

    Returns
    -------
    Relative permittivity of water (dielectric constant)

    References
    ----------
    Bradley, D.J.; Pitzer, K.S. `Thermodynamics of electrolytes. 12. Dielectric
        properties of water and Debye--Hueckel parameters to 350/sup
        0/C and 1 kbar`, J. Phys. Chem.; Journal Volume 83 (12)
        (1979), pp. 1599-1603,
        http://pubs.acs.org/doi/abs/10.1021/j100475a009
        DOI: 10.1021/j100475a009
    """
    be = get_backend(backend)
    if units is None:
        K = 1
        bar = 1
    else:
        K = units.kelvin
        bar = units.bar
    if T is None:
        T = 298.15*K
    if P is None:
        P = 1*bar
    if U is None:
        U = (3.4279e2,
             -5.0866e-3 / K,
             9.4690e-7 / K**2,
             -2.0525,
             3.1159e3*K,

def extended_activity_product(I, stoich, z, a, T, eps_r, rho, C=0, backend=None):
    be = get_backend(backend)
    Aval = A(eps_r, T, rho)
    Bval = B(eps_r, T, rho)
    tot = 0
    for idx, nr in enumerate(stoich):
        tot += nr*extended_log_gamma(I, z[idx], a[idx], Aval, Bval, C)
    return be.exp(tot)