def _proj_l1_sortsum(v, gamma, axis=None):
    r"""Projection operator of the :math:`\ell_1` norm. The solution is
    computed via the method of :cite:`duchi-2008-efficient`.


    Parameters
    ----------
    v : array_like
      Input array :math:`\mathbf{v}`
    gamma : float
      Parameter :math:`\gamma`
    axis : None or int, optional (default None)
      Axes of `v` over which to compute the :math:`\ell_1` norm. If
      `None`, an entire multi-dimensional array is treated as a
      vector. If axes are specified, then distinct norm values are
      computed over the indices of the remaining axes of input array
      `v`. **Note:** specifying a tuple of ints is not supported by
      this function.

    Returns
    -------
    x : ndarray
      Output array
    """

    if axis is None and norm_l1(v) <= gamma:
        return v
    if axis is not None and axis < 0:
        axis = v.ndim + axis
    av = np.abs(v)
    vs = np.sort(av, axis=axis)
    if axis is None:
        N = v.size
        c = 1.0 / np.arange(1, N + 1, dtype=v.dtype).reshape(v.shape)
        vs = vs[::-1].reshape(v.shape)
    else:
        N = v.shape[axis]
        ns = [v.shape[k] if k == axis else 1 for k in range(v.ndim)]
        c = 1.0 / np.arange(1, N + 1, dtype=v.dtype).reshape(ns)
        vs = vs[(slice(None),) * axis + (slice(None, None, -1),)]
    t = c * (np.cumsum(vs, axis=axis).reshape(v.shape) - gamma)
    K = np.sum(vs >= t, axis=axis, keepdims=True)
    t = (np.sum(vs * (vs >= t), axis=axis, keepdims=True) - gamma) / K
    t = np.asarray(np.maximum(0, t), dtype=v.dtype)
    return np.sign(v) * np.where(av > t, av - t, 0)