<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def proximal_convex_conj_kl(space, lam=1, g=None):
 r"""Proximal operator factory of the convex conjugate of the KL divergence. Function returning the proximal operator of the convex conjugate of the functional F where F is the entropy-type Kullback-Leibler (KL) divergence:: F(x) = sum_i (x_i - g_i + g_i ln(g_i) - g_i ln(pos(x_i))) + ind_P(x) with ``x`` and ``g`` elements in the linear space ``X``, and ``g`` non-negative. Here, ``pos`` denotes the nonnegative part, and ``ind_P`` is the indicator function for nonnegativity. Parameters space : `TensorSpace` Space X which is the domain of the functional F lam : positive float, optional Scaling factor. g : ``space`` element, optional Data term, positive. If None it is take as the one-element. Returns ------- prox_factory : function Factory for the proximal operator to be initialized. See Also -------- proximal_convex_conj_kl_cross_entropy : proximal for releated functional Notes ----- The functional is given by the expression .. math:: F(x) = \sum_i (x_i - g_i + g_i \ln(g_i) - g_i \ln(pos(x_i))) + I_{x \geq 0}(x) The indicator function :math:`I_{x \geq 0}(x)` is used to restrict the domain of :math:`F` such that :math:`F` is defined over whole space :math:`X`. The non-negativity thresholding :math:`pos` is used to define :math:`F` in the real numbers. Note that the functional is not well-defined without a prior g. Hence, if g is omitted this will be interpreted as if g is equal to the one-element. The convex conjugate :math:`F^*` of :math:`F` is .. math:: F^*(p) = \sum_i (-g_i \ln(\text{pos}({1_X}_i - p_i))) + I_{1_X - p \geq 0}(p) where :math:`p` is the variable dual to :math:`x`, and :math:`1_X` is an element of the space :math:`X` with all components set to 1. The proximal operator of the convex conjugate of F is .. math:: \mathrm{prox}_{\sigma (\lambda F)^*}(x) = \frac{\lambda 1_X + x - \sqrt{(x - \lambda 1_X)^2 + 4 \lambda \sigma g}}{2} where :math:`\sigma` is the step size-like parameter, and :math:`\lambda` is the weighting in front of the function :math:`F`. KL based objectives are common in MLEM optimization problems and are often used when data noise governed by a multivariate Poisson probability distribution is significant. The intermediate image estimates can have negative values even though the converged solution will be non-negative. Non-negative intermediate image estimates can be enforced by adding an indicator function ind_P the primal objective. This functional :math:`F`, described above, is related to the Kullback-Leibler cross entropy functional. The KL cross entropy is the one described in `this Wikipedia article <https://en.wikipedia.org/wiki/Kullback%E2%80%93Leibler_divergence>`_, and the functional :math:`F` is obtained by switching place of the prior and the varialbe in the KL cross entropy functional. See the See Also section. """

    lam = float(lam)

    if g is not None and g not in space:
        raise TypeError('{} is not an element of {}'.format(g, space))

    class ProximalConvexConjKL(Operator):

        """Proximal operator of the convex conjugate of the KL divergence."""

        def __init__(self, sigma):
            """Initialize a new instance.

            Parameters
            ----------
            sigma : positive float
            """
            super(ProximalConvexConjKL, self).__init__(
                domain=space, range=space, linear=False)
            self.sigma = float(sigma)

        def _call(self, x, out):
            """Return ``self(x, out=out)``."""
            # (x + lam - sqrt((x - lam)^2 + 4*lam*sig*g)) / 2

            # out = (x - lam)^2
            if x is out:
                # Handle aliased `x` and `out` (need original `x` later on)
                x = x.copy()
            else:
                out.assign(x)
            out -= lam
            out.ufuncs.square(out=out)

            # out = ... + 4*lam*sigma*g
            # If g is None, it is taken as the one element
            if g is None:
                out += 4.0 * lam * self.sigma
            else:
                out.lincomb(1, out, 4.0 * lam * self.sigma, g)

            # out = x - sqrt(...) + lam
            out.ufuncs.sqrt(out=out)
            out.lincomb(1, x, -1, out)
            out += lam

            # out = 1/2 * ...
            out /= 2

    return ProximalConvexConjKL