import itertools
import numpy as np

def generate_solutions_randomly(M):
    n, m = M.shape
    row_totals = M.sum(axis = 1)
    col_totals = M.sum(axis = 0)
    perms = []
    while 1:
        np.random.shuffle(M)
        if np.equal(M.sum(axis = 1), row_totals).all() and \
           np.equal(M.sum(axis = 0), col_totals).all():
            yield M


if __name__ == '__main__':

    # Demo
    # Generate random binary 9 by 14 matrix
    M = ( np.random.rand(9,14) > .5) * 1
    print M
    print

    # Find first 10 solutions using the generate_solutions iterator
    solutions = set()
    i = 0
    for sol in generate_solutions_randomly(M):
        solutions.add(str(sol))
        if len(solutions) != i:
            i = len(solutions)
            print sol
            print

        if len(solutions) == 10:
            break