<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def logsumexp(a, axis=None, b=None, use_numexpr=True):
 """Compute the log of the sum of exponentials of input elements. Parameters a : array_like Input array. axis : None or int, optional, default=None Axis or axes over which the sum is taken. By default `axis` is None, and all elements are summed. b : array-like, optional Scaling factor for exp(`a`) must be of the same shape as `a` or broadcastable to `a`. use_numexpr : bool, optional, default=True If True, use the numexpr library to speed up the calculation, which can give a 2-4X speedup when working with large arrays. Returns ------- res : ndarray The result, ``log(sum(exp(a)))`` calculated in a numerically more stable way. If `b` is given then ``log(sum(b*exp(a)))`` is returned. See Also -------- numpy.logaddexp, numpy.logaddexp2, scipy.misc.logsumexp (soon to be replaced with scipy.special.logsumexp) Notes ----- This is based on scipy.misc.logsumexp but with optional numexpr support for improved performance. """


    a = np.asarray(a)

    a_max = np.amax(a, axis=axis, keepdims=True)

    if a_max.ndim > 0:
        a_max[~np.isfinite(a_max)] = 0
    elif not np.isfinite(a_max):
        a_max = 0

    if b is not None:
        b = np.asarray(b)
        if use_numexpr and HAVE_NUMEXPR:
            out = np.log(numexpr.evaluate("b * exp(a - a_max)").sum(axis))
        else:
            out = np.log(np.sum(b * np.exp(a - a_max), axis=axis))
    else:
        if use_numexpr and HAVE_NUMEXPR:
            out = np.log(numexpr.evaluate("exp(a - a_max)").sum(axis))
        else:
            out = np.log(np.sum(np.exp(a - a_max), axis=axis))

    a_max = np.squeeze(a_max, axis=axis)
    out += a_max

    return out