def spdhg_generic(x, f, g, A, tau, sigma, niter, **kwargs):
    r"""Computes a saddle point with a stochastic PDHG.

    This means, a solution (x*, y*), y* = (y*_1, ..., y*_n) such that

    (x*, y*) in arg min_x max_y sum_i=1^n <y_i, A_i> - f*[i](y_i) + g(x)

    where g : X -> IR_infty and f[i] : Y[i] -> IR_infty are convex, l.s.c. and
    proper functionals. For this algorithm, they all may be non-smooth and no
    strong convexity is assumed.

    Parameters
    ----------
    x : primal variable
        This variable is both input and output of the method.
    f : functions
        Functionals Y[i] -> IR_infty that all have a convex conjugate with a
        proximal operator, i.e.
        f[i].convex_conj.proximal(sigma[i]) : Y[i] -> Y[i].
    g : function
        Functional X -> IR_infty that has a proximal operator, i.e.
        g.proximal(tau) : X -> X.
    A : functions
        Operators A[i] : X -> Y[i] that possess adjoints: A[i].adjoint
    tau : scalar / vector / matrix
        Step size for primal variable. Note that the proximal operator of g
        has to be well-defined for this input.
    sigma : scalar
        Scalar / vector / matrix used as step size for dual variable. Note that
        the proximal operator related to f (see above) has to be well-defined
        for this input.
    niter : int
        Number of iterations

    Other Parameters
    ----------------
    y : dual variable, optional
        Dual variable is part of a product space. By default equals 0.
    z : variable, optional
        Adjoint of dual variable, z = A^* y. By default equals 0 if y = 0.
    mu_g : scalar
        Strong convexity constant of g.
    theta : scalar
        Global extrapolation factor.
    extra: list
        List of local extrapolation paramters for every index i. By default
        extra_i = 1.
    fun_select : function
        Function that selects blocks at every iteration IN -> {1,...,n}. By
        default this is serial uniform sampling, fun_select(k) selects an index
        i \in {1,...,n} with probability 1/n.
    callback : callable, optional
        Function called with the current iterate after each iteration.

    References
    ----------
    [CERS2017] A. Chambolle, M. J. Ehrhardt, P. Richtarik and C.-B. Schoenlieb,
    *Stochastic Primal-Dual Hybrid Gradient Algorithm with Arbitrary Sampling
    and Imaging Applications*. ArXiv: http://arxiv.org/abs/1706.04957 (2017).

    [E+2017] M. J. Ehrhardt, P. J. Markiewicz, P. Richtarik, J. Schott,
    A. Chambolle and C.-B. Schoenlieb, *Faster PET reconstruction with a
    stochastic primal-dual hybrid gradient method*. Wavelets and Sparsity XVII,
    58 (2017) http://doi.org/10.1117/12.2272946.
    """
    # Callback object
    callback = kwargs.pop('callback', None)
    if callback is not None and not callable(callback):
        raise TypeError('`callback` {} is not callable'
                        ''.format(callback))

    # Dual variable
    y = kwargs.pop('y', None)
    if y is None:
        y = A.range.zero()

    # Adjoint of dual variable
    z = kwargs.pop('z', None)
    if z is None:
        if y.norm() == 0:
            z = A.domain.zero()
        else:
            z = A.adjoint(y)

    # Strong convexity of g
    mu_g = kwargs.pop('mu_g', None)
    if mu_g is None:
        update_proximal_primal = False
    else:
        update_proximal_primal = True

    # Global extrapolation factor theta
    theta = kwargs.pop('theta', 1)

    # Second extrapolation factor
    extra = kwargs.pop('extra', None)
    if extra is None:
        extra = [1] * len(sigma)

    # Selection function
    fun_select = kwargs.pop('fun_select', None)
    if fun_select is None:
        def fun_select(x):
            return [int(np.random.choice(len(A), 1, p=1 / len(A)))]

    # Initialize variables
    z_relax = z.copy()
    dz = A.domain.element()
    y_old = A.range.element()

    # Save proximal operators
    proximal_dual_sigma = [fi.convex_conj.proximal(si)
                           for fi, si in zip(f, sigma)]
    proximal_primal_tau = g.proximal(tau)

    # run the iterations
    for k in range(niter):

        # select block
        selected = fun_select(k)

        # update primal variable
        # tmp = x - tau * z_relax; z_relax used as tmp variable
        z_relax.lincomb(1, x, -tau, z_relax)
        # x = prox(tmp)
        proximal_primal_tau(z_relax, out=x)

        # update extrapolation parameter theta
        if update_proximal_primal:
            theta = float(1 / np.sqrt(1 + 2 * mu_g * tau))

        # update dual variable and z, z_relax
        z_relax.assign(z)
        for i in selected:

            # save old yi
            y_old[i].assign(y[i])

            # tmp = Ai(x)
            A[i](x, out=y[i])

            # tmp = y_old + sigma_i * Ai(x)
            y[i].lincomb(1, y_old[i], sigma[i], y[i])

            # y[i]= prox(tmp)
            proximal_dual_sigma[i](y[i], out=y[i])

            # update adjoint of dual variable
            y_old[i].lincomb(-1, y_old[i], 1, y[i])
            A[i].adjoint(y_old[i], out=dz)
            z += dz

            # compute extrapolation
            z_relax.lincomb(1, z_relax, 1 + theta * extra[i], dz)

        # update the step sizes tau and sigma for acceleration
        if update_proximal_primal:
            for i in range(len(sigma)):
                sigma[i] /= theta
            tau *= theta

            proximal_dual_sigma = [fi.convex_conj.proximal(si)
                                   for fi, si in zip(f, sigma)]
            proximal_primal_tau = g.proximal(tau)

        if callback is not None:
            callback([x, y])