code: def proj_l1(v, gamma, axis=None, method=None):
      r
  
      if method is None:
          if axis is None:
              method = 'scalarroot'
          else:
              method = 'sortcumsum'
  
      if method == 'scalarroot':
          if axis is not None:
              raise ValueError('Method scalarroot only supports axis=None')
          return _proj_l1_scalar_root(v, gamma)
      elif method == 'sortcumsum':
          if isinstance(axis, tuple):
              vtr, rsi = ndto2d(v, axis)
              xtr = _proj_l1_sortsum(vtr, gamma, axis=1)
              return ndfrom2d(xtr, rsi)
          else:
              return _proj_l1_sortsum(v, gamma, axis)
      else:
          raise ValueError('Unknown solver method %s' % method)
docs: r"""Projection operator of the :math:`\ell_1` norm.
  
  
      Parameters
      ----------
      v : array_like
        Input array :math:`\mathbf{v}`
      gamma : float
        Parameter :math:`\gamma`
      axis : None or int or tuple of ints, 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`.
      method : None or str, optional (default None)
        Solver method to use. If `None`, the most appropriate choice is
        made based on the `axis` parameter. Valid methods are
  
           - 'scalarroot'
              The solution is computed via the method of Sec. 6.5.2 in
              :cite:`parikh-2014-proximal`.
           - 'sortcumsum'
              The solution is computed via the method of
              :cite:`duchi-2008-efficient`.
  
      Returns
      -------
      x : ndarray
        Output array