code: def solve_semi_dual(a, b, C, regul, method="L-BFGS-B", tol=1e-3, max_iter=500,
                      verbose=False):
      
  
      def _func(alpha):
          obj, grad = semi_dual_obj_grad(alpha, a, b, C, regul)
          # We need to maximize the semi-dual.
          return -obj, -grad
  
      alpha_init = np.zeros(len(a))
  
      res = minimize(_func, alpha_init, method=method, jac=True,
                     tol=tol, options=dict(maxiter=max_iter, disp=verbose))
  
      return res.x, res
docs: Solve the "smoothed" semi-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 max_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)
          Semi-dual potentials.