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
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]])