code: def dual_obj_grad(alpha, beta, a, b, C, regul):
      
      obj = np.dot(alpha, a) + np.dot(beta, b)
      grad_alpha = a.copy()
      grad_beta = b.copy()
  
      # X[:, j] = alpha + beta[j] - C[:, j]
      X = alpha[:, np.newaxis] + beta - C
  
      # val.shape = len(b)
      # G.shape = len(a) x len(b)
      val, G = regul.delta_Omega(X)
  
      obj -= np.sum(val)
      grad_alpha -= G.sum(axis=1)
      grad_beta -= G.sum(axis=0)
  
      return obj, grad_alpha, grad_beta
docs: Compute objective value and gradients of dual objective.
  
      Parameters
      ----------
      alpha: array, shape = len(a)
      beta: array, shape = len(b)
          Current iterate of dual potentials.
      a: array, shape = len(a)
      b: array, shape = len(b)
          Input histograms (should be non-negative and sum to 1).
      C: array, shape = len(a) x len(b)
          Ground cost matrix.
      regul: Regularization object
          Should implement a delta_Omega(X) method.
  
      Returns
      -------
      obj: float
          Objective value (higher is better).
      grad_alpha: array, shape = len(a)
          Gradient w.r.t. alpha.
      grad_beta: array, shape = len(b)
          Gradient w.r.t. beta.