docs: def prox_gradf_lim12(x, step, j=None, Xs=None, boundary=None):
      """1D projection operator"""
      # TODO: split boundary in x1 and x2 and use appropriate operator
      if j == 0:
          x -= step*grad_fx(Xs[0][0], Xs[1][0])
      if j == 1:
          y = x
          y -= step*grad_fy(Xs[0][0], Xs[1][0])
      return prox_lim12(x, step, j=j, Xs=Xs, boundary=boundary)