def birkhoff_von_neumann_decomposition(D):
m, n = D.shape<EOL>if m != n:<EOL><INDENT>raise ValueError('<STR_LIT>'.format(m, n))<EOL><DEDENT>indices = list(itertools.product(range(m), range(n)))<EOL>coefficients = []<EOL>permutations = []<EOL>S = D.astype('<STR_LIT:float>')<EOL>while not np.all(S == <NUM_LIT:0>):<EOL><INDENT>W = to_pattern_matrix(S)<EOL>X = to_bipartite_matrix(W)<EOL>G = from_numpy_matrix(X)<EOL>left_nodes = range(n)<EOL>M = maximum_matching(G, left_nodes)<EOL>M = {u: v % n for u, v in M.items() if u < n}<EOL>P = to_permutation_matrix(M)<EOL>q = min(S[i, j] for (i, j) in indices if P[i, j] == <NUM_LIT:1>)<EOL>coefficients.append(q)<EOL>permutations.append(P)<EOL>S -= q * P<EOL>S[np.abs(S) < TOLERANCE] = <NUM_LIT:0.0><EOL><DEDENT>return list(zip(coefficients, permutations))<EOL>
Returns the Birkhoff--von Neumann decomposition of the doubly
    stochastic matrix `D`.

    The input `D` must be a square NumPy array representing a doubly
    stochastic matrix (that is, a matrix whose entries are nonnegative
    reals and whose row sums and column sums are all 1). Each doubly
    stochastic matrix is a convex combination of at most ``n ** 2``
    permutation matrices, where ``n`` is the dimension of the input
    array.

    The returned value is a list of pairs whose length is at most ``n **
    2``. In each pair, the first element is a real number in the interval **(0,
    1]** and the second element is a NumPy array representing a permutation
    matrix. This represents the doubly stochastic matrix as a convex
    combination of the permutation matrices.

    The input matrix may also be a scalar multiple of a doubly
    stochastic matrix, in which case the row sums and column sums must
    each be *c*, for some positive real number *c*. This may be useful
    in avoiding precision issues: given a doubly stochastic matrix that
    will have many entries close to one, multiply it by a large positive
    integer. The returned permutation matrices will be the same
    regardless of whether the given matrix is a doubly stochastic matrix
    or a scalar multiple of a doubly stochastic matrix, but in the
    latter case, the coefficients will all be scaled by the appropriate
    scalar multiple, and their sum will be that scalar instead of one.

    For example::

        >>> import numpy as np
        >>> from birkhoff import birkhoff_von_neumann_decomposition as decomp
        >>> D = np.ones((2, 2))
        >>> zipped_pairs = decomp(D)
        >>> coefficients, permutations = zip(*zipped_pairs)
        >>> coefficients
        (1.0, 1.0)
        >>> permutations[0]
        array([[ 1.,  0.],
               [ 0.,  1.]])
        >>> permutations[1]
        array([[ 0.,  1.],
               [ 1.,  0.]])
        >>> zipped_pairs = decomp(D / 2)  # halve each value in the matrix
        >>> coefficients, permutations = zip(*zipped_pairs)
        >>> coefficients  # will be half as large as before
        (0.5, 0.5)
        >>> permutations[0]  # will be the same as before
        array([[ 1.,  0.],
               [ 0.,  1.]])
        >>> permutations[1]
        array([[ 0.,  1.],
               [ 1.,  0.]])

    The returned list of pairs is given in the order computed by the algorithm
    (so in particular they are not sorted in any way).