docs: def transition_matrix(C, reversible=False, mu=None, method='auto', **kwargs):
      r"""Estimate the transition matrix from the given countmatrix.
  
      Parameters
      ----------
      C : numpy ndarray or scipy.sparse matrix
          Count matrix
      reversible : bool (optional)
          If True restrict the ensemble of transition matrices
          to those having a detailed balance symmetry otherwise
          the likelihood optimization is carried out over the whole
          space of stochastic matrices.
      mu : array_like
          The stationary distribution of the MLE transition matrix.
      method : str
          Select which implementation to use for the estimation.
          One of 'auto', 'dense' and 'sparse', optional, default='auto'.
          'dense' always selects the dense implementation, 'sparse' always selects
          the sparse one.
          'auto' selects the most efficient implementation according to
          the sparsity structure of the matrix: if the occupation of the C
          matrix is less then one third, select sparse. Else select dense.
          The type of the T matrix returned always matches the type of the
          C matrix, irrespective of the method that was used to compute it.
      **kwargs: Optional algorithm-specific parameters. See below for special cases
      Xinit : (M, M) ndarray
          Optional parameter with reversible = True.
          initial value for the matrix of absolute transition probabilities. Unless set otherwise,
          will use X = diag(pi) t, where T is a nonreversible transition matrix estimated from C,
          i.e. T_ij = c_ij / sum_k c_ik, and pi is its stationary distribution.
      maxiter : 1000000 : int
          Optional parameter with reversible = True.
          maximum number of iterations before the method exits
      maxerr : 1e-8 : float
          Optional parameter with reversible = True.
          convergence tolerance for transition matrix estimation.
          This specifies the maximum change of the Euclidean norm of relative
          stationary probabilities (:math:`x_i = \sum_k x_{ik}`). The relative stationary probability changes
          :math:`e_i = (x_i^{(1)} - x_i^{(2)})/(x_i^{(1)} + x_i^{(2)})` are used in order to track changes in small
          probabilities. The Euclidean norm of the change vector, :math:`|e_i|_2`, is compared to maxerr.
      rev_pisym : bool, default=False
          Fast computation of reversible transition matrix by normalizing
          :math:`x_{ij} = \pi_i p_{ij} + \pi_j p_{ji}`. :math:`p_{ij}` is the direct
          (nonreversible) estimate and :math:`pi_i` is its stationary distribution.
          This estimator is asympotically unbiased but not maximum likelihood.
      return_statdist : bool, default=False
          Optional parameter with reversible = True.
          If set to true, the stationary distribution is also returned
      return_conv : bool, default=False
          Optional parameter with reversible = True.
          If set to true, the likelihood history and the pi_change history is returned.
      warn_not_converged : bool, default=True
          Prints a warning if not converged.
      sparse_newton : bool, default=False
          If True, use the experimental primal-dual interior-point solver for sparse input/computation method.
  
      Returns
      -------
      P : (M, M) ndarray or scipy.sparse matrix
         The MLE transition matrix. P has the same data type (dense or sparse)
         as the input matrix C.
      The reversible estimator returns by default only P, but may also return
      (P,pi) or (P,lhist,pi_changes) or (P,pi,lhist,pi_changes) depending on the return settings
      P : ndarray (n,n)
          transition matrix. This is the only return for return_statdist = False, return_conv = False
      (pi) : ndarray (n)
          stationary distribution. Only returned if return_statdist = True
      (lhist) : ndarray (k)
          likelihood history. Has the length of the number of iterations needed.
          Only returned if return_conv = True
      (pi_changes) : ndarray (k)
          history of likelihood history. Has the length of the number of iterations needed.
          Only returned if return_conv = True
  
      Notes
      -----
      The transition matrix is a maximum likelihood estimate (MLE) of
      the probability distribution of transition matrices with
      parameters given by the count matrix.
  
      References
      ----------
      .. [1] Prinz, J H, H Wu, M Sarich, B Keller, M Senne, M Held, J D
          Chodera, C Schuette and F Noe. 2011. Markov models of
          molecular kinetics: Generation and validation. J Chem Phys
          134: 174105
      .. [2] Bowman, G R, K A Beauchamp, G Boxer and V S Pande. 2009.
          Progress and challenges in the automated construction of Markov state models for full protein systems.
          J. Chem. Phys. 131: 124101
      .. [3] Trendelkamp-Schroer, B, H Wu, F Paul and F. Noe. 2015
          Estimation and uncertainty of reversible Markov models.
          J. Chem. Phys. 143: 174101
  
      Examples
      --------
  
      >>> import numpy as np
      >>> from msmtools.estimation import transition_matrix
  
      >>> C = np.array([[10, 1, 1], [2, 0, 3], [0, 1, 4]])
  
      Non-reversible estimate
  
      >>> T_nrev = transition_matrix(C)
      >>> T_nrev
      array([[ 0.83333333,  0.08333333,  0.08333333],
             [ 0.4       ,  0.        ,  0.6       ],
             [ 0.        ,  0.2       ,  0.8       ]])
  
      Reversible estimate
  
      >>> T_rev = transition_matrix(C, reversible=True)
      >>> T_rev
      array([[ 0.83333333,  0.10385551,  0.06281115],
             [ 0.35074677,  0.        ,  0.64925323],
             [ 0.04925323,  0.15074677,  0.8       ]])
  
      Reversible estimate with given stationary vector
  
      >>> mu = np.array([0.7, 0.01, 0.29])
      >>> T_mu = transition_matrix(C, reversible=True, mu=mu)
      >>> T_mu
      array([[ 0.94771371,  0.00612645,  0.04615984],
             [ 0.42885157,  0.        ,  0.57114843],
             [ 0.11142031,  0.01969477,  0.86888491]])
  
      """
      if issparse(C):
          sparse_input_type = True
      elif isdense(C):
          sparse_input_type = False
      else:
          raise NotImplementedError('C has an unknown type.')
  
      if method == 'dense':
          sparse_computation = False
      elif method == 'sparse':
          sparse_computation = True
      elif method == 'auto':
          # heuristically determine whether is't more efficient to do a dense of sparse computation
          if sparse_input_type:
              dof = C.getnnz()
          else:
              dof = np.count_nonzero(C)
          dimension = C.shape[0]
          if dimension*dimension < 3*dof:
              sparse_computation = False
          else:
              sparse_computation = True
      else:
          raise ValueError(('method="%s" is no valid choice. It should be one of'
                            '"dense", "sparse" or "auto".') % method)
  
      # convert input type
      if sparse_computation and not sparse_input_type:
          C = coo_matrix(C)
      if not sparse_computation and sparse_input_type:
          C = C.toarray()
  
      return_statdist = 'return_statdist' in kwargs
  
      if not return_statdist:
          kwargs['return_statdist'] = False
  
      sparse_newton = kwargs.pop('sparse_newton', False)
  
      if reversible:
          rev_pisym = kwargs.pop('rev_pisym', False)
  
          if mu is None:
              if sparse_computation:
                  if rev_pisym:
                      result = sparse.transition_matrix.transition_matrix_reversible_pisym(C, **kwargs)
                  elif sparse_newton:
                      from msmtools.estimation.sparse.newton.mle_rev import solve_mle_rev
                      result = solve_mle_rev(C, **kwargs)
                  else:
                      result = sparse.mle_trev.mle_trev(C, **kwargs)
              else:
                  if rev_pisym:
                      result = dense.transition_matrix.transition_matrix_reversible_pisym(C, **kwargs)
                  else:
                      result = dense.mle_trev.mle_trev(C, **kwargs)
          else:
              kwargs.pop('return_statdist') # pi given, keyword unknown by estimators.
              if sparse_computation:
                  # Sparse, reversible, fixed pi (currently using dense with sparse conversion)
                  result = sparse.mle_trev_given_pi.mle_trev_given_pi(C, mu, **kwargs)
              else:
                  result = dense.mle_trev_given_pi.mle_trev_given_pi(C, mu, **kwargs)
      else:  # nonreversible estimation
          if mu is None:
              if sparse_computation:
                  # Sparse,  nonreversible
                  result = sparse.transition_matrix.transition_matrix_non_reversible(C)
              else:
                  # Dense,  nonreversible
                  result = dense.transition_matrix.transition_matrix_non_reversible(C)
              # Both methods currently do not have an iterate of pi, so we compute it here for consistency.
              if return_statdist:
                  from msmtools.analysis import stationary_distribution
                  mu = stationary_distribution(result)
          else:
              raise NotImplementedError('nonreversible mle with fixed stationary distribution not implemented.')
  
      if return_statdist and isinstance(result, tuple):
          T, mu = result
      else:
          T = result
  
      # convert return type
      if sparse_computation and not sparse_input_type:
          T = T.toarray()
      elif not sparse_computation and sparse_input_type:
          T = csr_matrix(T)
  
      if return_statdist:
          return T, mu
      return T