code: def poissreg(x0, rho, x, y):
      
  
      # objective and gradient
      n = float(x.shape[0])
      f = lambda w: np.mean(np.exp(x.dot(w)) - y * x.dot(w))
      df = lambda w: (x.T.dot(np.exp(x.dot(w))) - x.T.dot(y)) / n
  
      # minimize via BFGS
      return bfgs(x0, rho, f, df)
docs: Proximal operator for Poisson regression
  
      Computes the proximal operator of the negative log-likelihood loss assumping a Poisson noise distribution.
  
      Parameters
      ----------
      x0 : array_like
          The starting or initial point used in the proximal update step
  
      rho : float
          Momentum parameter for the proximal step (larger value -> stays closer to x0)
  
      x : (n, k) array_like
          A design matrix consisting of n examples of k-dimensional features (or input).
  
      y : (n,) array_like
          A vector containing the responses (outupt) to the n features given in x.
  
      Returns
      -------
      theta : array_like
          The parameter vector found after running the proximal update step