code: def solve_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500,
                 verbose=False):
      
  
      def _func(params):
          # Unpack alpha and beta.
          alpha = params[:len(a)]
          beta = params[len(a):]
  
          obj, grad_alpha, grad_beta = dual_obj_grad(alpha, beta, a, b, C, regul)
  
          # Pack grad_alpha and grad_beta.
          grad = np.concatenate((grad_alpha, grad_beta))
  
          # We need to maximize the dual.
          return -obj, -grad
  
      # Unfortunately, `minimize` only supports functions whose argument is a
      # vector. So, we need to concatenate alpha and beta.
      alpha_init = np.zeros(len(a))
      beta_init = np.zeros(len(b))
      params_init = np.concatenate((alpha_init, beta_init))
  
      res = minimize(_func, params_init, method=method, jac=True,
                     tol=tol, options=dict(maxiter=max_iter, disp=verbose))
  
      alpha = res.x[:len(a)]
      beta = res.x[len(a):]
  
      return alpha, beta, res
docs: Solve the "smoothed" dual objective.
  
      Parameters
      ----------
      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.
      method: str
          Solver to be used (passed to `scipy.optimize.minimize`).
      tol: float
          Tolerance parameter.
      max_iter: int
          Maximum number of iterations.
  
      Returns
      -------
      alpha: array, shape = len(a)
      beta: array, shape = len(b)
          Dual potentials.