docs: def solve_mle_rev(C, tol=1e-10, maxiter=100, show_progress=False, full_output=False,
                    return_statdist=True, **kwargs):
      """Number of states"""
      M = C.shape[0]
  
      """Initial guess for primal-point"""
      z0 = np.zeros(2*M)
      z0[0:M] = 1.0
  
      """Inequality constraints"""
      # G = np.zeros((M, 2*M))
      # G[np.arange(M), np.arange(M)] = -1.0
      G = -1.0*scipy.sparse.eye(M, n=2*M, k=0)
      h = np.zeros(M)
  
      """Equality constraints"""
      A = np.zeros((1, 2*M))
      A[0, M] = 1.0
      b = np.array([0.0])
  
      """Scaling"""
      c0 = C.max()
      C = C/c0
  
      """Symmetric part"""
      Cs = C + C.T
  
      """Column sum"""
      c = C.sum(axis=0)
  
      if scipy.sparse.issparse(C):
          Cs = Cs.tocsr()
          c = c.A1
          A = scipy.sparse.csr_matrix(A)
          F = objective_sparse.F
          DF = objective_sparse.DFsym
          convert_solution = objective_sparse.convert_solution
      else:
          F = objective_dense.F
          DF = objective_dense.DF
          convert_solution = objective_dense.convert_solution
  
      """PDIP iteration"""
      res = primal_dual_solve(F, z0, DF, A, b, G, h,
                              args=(Cs, c),
                              maxiter=maxiter, tol=tol,
                              show_progress=show_progress,
                              full_output=full_output)
      if full_output:
          z, info = res
      else:
          z = res
      pi, P = convert_solution(z, Cs)
  
      result = [P]
      if return_statdist:
          result.append(pi)
      if full_output:
          result.append(info)
  
      return tuple(result) if len(result) > 1 else result[0]