<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def osmlem(op, x, data, niter, callback=None, **kwargs):
 r"""Ordered Subsets Maximum Likelihood Expectation Maximation algorithm. This solver attempts to solve:: max_x L(x | data) where ``L(x, | data)`` is the likelihood of ``x`` given ``data``. The that (approximately):
: op[i](x) = data[i] Parameters op : sequence of `Operator` Forward operators in the inverse problem. x : ``op.domain`` element Vector to which the result is written. Its initial value is used as starting point of the iteration, and its values are updated in each iteration step. The initial value of ``x`` should be non-negative. data : sequence of ``op.range`` `element-like` Right-hand sides of the equation defining the inverse problem. niter : int Number of iterations. callback : callable, optional Function called with the current iterate after each iteration. Other Parameters sensitivities : float or ``op.domain`` `element-like`, optional The algorithm contains an ``A^T 1`` term, if this parameter is given, it is replaced by it. Default: ``op[i].adjoint(op[i].range.one())`` Notes ----- Given forward models :math:`A_i`, and data :math:`g_i`, the algorithm attempts to find an :math:`x` that maximizes: .. math:: \prod_{i=1}^M P(g_i | g_i \text{ is } Poisson(A_i(x)) \text{ distributed}). The algorithm is explicitly given by partial updates: .. math:: x_{n + m/M} = \frac{x_{n + (m - 1)/M}}{A_i^* 1} A_i^* (g_i / A_i(x_{n + (m - 1)/M})) The algorithm is not guaranteed to converge, but works for many practical problems. References Natterer, F. Mathematical Methods in Image Reconstruction, section 5.3.2. See Also -------- mlem : Ordinary MLEM algorithm without subsets. loglikelihood : Function for calculating the logarithm of the likelihood """

    n_ops = len(op)
    if len(data) != n_ops:
        raise ValueError('number of data ({}) does not match number of '
                         'operators ({})'.format(len(data), n_ops))
    if not all(x in opi.domain for opi in op):
        raise ValueError('`x` not an element in the domains of all operators')

    # Convert data to range elements
    data = [op[i].range.element(data[i]) for i in range(len(op))]

    # Parameter used to enforce positivity.
    # TODO: let users give this.
    eps = 1e-8

    if np.any(np.less(x, 0)):
        raise ValueError('`x` must be non-negative')

    # Extract the sensitivites parameter
    sensitivities = kwargs.pop('sensitivities', None)
    if sensitivities is None:
        sensitivities = [np.maximum(opi.adjoint(opi.range.one()), eps)
                         for opi in op]
    else:
        # Make sure the sensitivities is a list of the correct size.
        try:
            list(sensitivities)
        except TypeError:
            sensitivities = [sensitivities] * n_ops

    tmp_dom = op[0].domain.element()
    tmp_ran = [opi.range.element() for opi in op]

    for _ in range(niter):
        for i in range(n_ops):
            op[i](x, out=tmp_ran[i])
            tmp_ran[i].ufuncs.maximum(eps, out=tmp_ran[i])
            data[i].divide(tmp_ran[i], out=tmp_ran[i])

            op[i].adjoint(tmp_ran[i], out=tmp_dom)
            tmp_dom /= sensitivities[i]

            x *= tmp_dom

            if callback is not None:
                callback(x)