<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def kaczmarz(ops, x, rhs, niter, omega=1, projection=None, random=False, callback=None, callback_loop='outer'):
 r"""Optimized implementation of Kaczmarz's method. Solves the inverse problem given by the set of equations:: A_n(x) = rhs_n This is also known as the Landweber-Kaczmarz's method, since the method coincides with the Landweber method for a single operator. Parameters ops : sequence of `Operator`'s Operators in the inverse problem. ``op[i].derivative(x).adjoint`` must be well-defined for ``x`` in the operator domain and for all ``i``. x : ``op.domain`` element Element 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. rhs : sequence of ``ops[i].range`` elements Right-hand side of the equation defining the inverse problem. niter : int Number of iterations. omega : positive float or sequence of positive floats, optional Relaxation parameter in the iteration. If a single float is given the same step is used for all operators, otherwise separate steps are used. projection : callable, optional Function that can be used to modify the iterates in each iteration, for example enforcing positivity. The function should take one argument and modify it in-place. random : bool, optional If `True`, the order of the operators is randomized in each iteration. callback : callable, optional Object executing code per iteration, e.g. plotting each iterate. callback_loop : {'inner', 'outer'} Whether the callback should be called in the inner or outer loop. Notes ----- This method calculates an approximate least-squares solution of the inverse problem of the first kind .. math:: \mathcal{A}_i (x) = y_i \quad 1 \leq i \leq n, for a given :math:`y_n \in \mathcal{Y}_n`, i.e. an approximate solution :math:`x^*` to .. math:: \min_{x\in \mathcal{X}} \sum_{i=1}^n \| \mathcal{A}_i(x) - y_i \|_{\mathcal{Y}_i}^2 for a (Frechet-) differentiable operator :math:`\mathcal{A}: \mathcal{X} \to \mathcal{Y}` between Hilbert spaces :math:`\mathcal{X}` and :math:`\mathcal{Y}`. The method starts from an initial guess :math:`x_0` and uses the iteration .. math:: x_{k+1} = x_k - \omega_{[k]} \ \partial \mathcal{A}_{[k]}(x_k)^* (\mathcal{A}_{[k]}(x_k) - y_{[k]}), where :math:`\partial \mathcal{A}_{[k]}(x_k)` is the Frechet derivative of :math:`\mathcal{A}_{[k]}` at :math:`x_k`, :math:`\omega_{[k]}` is a relaxation parameter and :math:`[k] := k \text{ mod } n`. For linear problems, a choice :math:`0 < \omega_i < 2/\lVert \mathcal{A}_{i}^2\rVert` guarantees convergence, where :math:`\|\mathcal{A}_{i}\|` stands for the operator norm of :math:`\mathcal{A}_{i}`. This implementation uses a minimum amount of memory copies by applying re-usable temporaries and in-place evaluation. The method is also described in a `Wikipedia article <https://en.wikipedia.org/wiki/Kaczmarz_method>`_. and in Natterer, F. Mathematical Methods in Image Reconstruction, section 5.3.2. See Also -------- landweber """

    domain = ops[0].domain
    if any(domain != opi.domain for opi in ops):
        raise ValueError('domains of `ops` are not all equal')

    if x not in domain:
        raise TypeError('`x` {!r} is not in the domain of `ops` {!r}'
                        ''.format(x, domain))

    if len(ops) != len(rhs):
        raise ValueError('`number of `ops` {} does not match number of '
                         '`rhs` {}'.format(len(ops), len(rhs)))

    omega = normalized_scalar_param_list(omega, len(ops), param_conv=float)

    # Reusable elements in the range, one per type of space
    ranges = [opi.range for opi in ops]
    unique_ranges = set(ranges)
    tmp_rans = {ran: ran.element() for ran in unique_ranges}

    # Single reusable element in the domain
    tmp_dom = domain.element()

    # Iteratively find solution
    for _ in range(niter):
        if random:
            rng = np.random.permutation(range(len(ops)))
        else:
            rng = range(len(ops))

        for i in rng:
            # Find residual
            tmp_ran = tmp_rans[ops[i].range]
            ops[i](x, out=tmp_ran)
            tmp_ran -= rhs[i]

            # Update x
            ops[i].derivative(x).adjoint(tmp_ran, out=tmp_dom)
            x.lincomb(1, x, -omega[i], tmp_dom)

            if projection is not None:
                projection(x)

            if callback is not None and callback_loop == 'inner':
                callback(x)
        if callback is not None and callback_loop == 'outer':
            callback(x)