automatic "lotjes trekken"

[Tim Peters]

[Carel Fellinger wades thru a mass of Tim's chessboard prose, and
 cleverly deduces that he still doesn't have a usable answer]

[Tim]
> ...
> The next twist in the tale is too involved to explain here:  it
> turns out that it's possible to deduce the result from the numbers of
> ways to place 0, 1, 2, ... and N non-capturing rooks on the "cheese-
> inverted" chessboard obtained by putting cheese on all the empty squares
> and removing the cheese from the cheesy squares.  In most applications,
> there are many more empty squares than cheesy squares, so inverting
> their roles can slash the amount of computation.

[Carel]
> and how then can this be deduced, what simple function are we to apply?

Who said "simple" <wink>?  If H(i) is the # of ways to place i non-capturing
rooks on the inverted chessboard, the number of ways to place N
non-capturing rooks on the original chessboard is the sum of the N+1 terms
of this form, where i ranges from 0 thru N inclusive:

    (-1)**i * H(i) * (N-i)! =

    H(0)*N! - H(1)*(N-1)! + H(2)*(N-2)! - ... + (-1)**N*H(N)

H(0) is 1 by definition (why?  because no other value works <wink>).

By way of non-trivial example, let's go back to the number of derangements
of N things.  Then the inverted chessboard is a degenerate one consisting of
cheese everywhere *except* on the main diagonal.  Since all the empty
squares are in distinct rows and columns, the number of ways to place i
non-capturing rooks on the N free diagonal squares is simply C(N, i), where
C(a, b) is the number of combinations of a things taken b at a time =
a!/b!/(a-b)!.  Plugging H(i)=C(N,i) into the above gives the sum of terms of
the form

    (-1)**i * C(N, i) * (N-i)! =
    (-1)**i * N! / i!

which is N! times the sum of N+1 terms of the form

    (-1)**i / i!

Note that the latter sum consists of the first N+1 terms of the Taylor
series expansion for exp(-1), so is an excellent approximation to 1/e.  And
that's "why" the number of derangements is approximately N!/e (the exact
number is gotten by cutting off the approximation to 1/e at the N+1'st
term).

I haven't explained why any of that *works*, and that's still too involved
to go into here.  Any good text on combinatorics will explain it in detail,
though!  In brief, it's a straightforward application of the "principle of
inclusion/exclusion" (which may be called that, or hiding in a chapter on
"sieve formulas") to the "rook polynomial" of the cheesy squares (which is a
simple flavor of generating function).  You already know where the "rook"
comes from, so the rest is a small matter of filling in details <wink>.

Anyway, forget all that unless it interests you.  The news->mail gateway has
disgorged other old msgs in this thread, from which I deduce you don't have
an application in mind for more than 10 people anyway.  I'll attach
something unrelated to the theory above that simply works, and indeed works
faster the more restrictions there are.

Here's an example showing why you should *not* take the other approaches
that have been suggested:

>>> b = Board()
>>> for name in "Carel", "Guido", "Martijn", "Tim":
        b.register_name(name)

>>> b.diagonalize()  # force derangement
>>> b.place_cheese("Carel", "Guido")   # Guido has enough wooden shoes
>>> b.place_cheese("Martijn", "Guido") # ditto
>>> for solution in b.solve():
        print solution

[('Carel', 'Martijn'), ('Guido', 'Carel'),
 ('Martijn', 'Tim'),   ('Tim', 'Guido')]
[('Carel', 'Martijn'), ('Guido', 'Tim'),
 ('Martijn', 'Carel'), ('Tim', 'Guido')]
[('Carel', 'Tim'),     ('Guido', 'Martijn'),
 ('Martijn', 'Carel'), ('Tim', 'Guido')]
>>>

That means there are exactly three permutations that avoid the forbidden
positions.  Note that in two of them, you have to tease Martijn, and in only
one do I enjoy the pleasure of your gift-selecting powers.

That leads to a (possibly) surprising observation:  if you draw lots "at
random" one at a time, you do *not* (in general) have an equal chance of
selecting each of the possible solutions!  In the above, if you draw your
victim first, you have an equal chance of picking me or Martijn.  But you
*should* pick Martijn *twice* as often as me if you're to have equal chances
of ending up with each of the 3 solutions!  If you choose Guido's victim
first, there's no such problem, and ditto if you pick my victim first
(albeit for a different reason) -- but you have no way to know that without
analyzing the set of solutions first.  So what's required isn't backtracking
but forwardtracking <wink>.

Here are three approaches that work:

1) A variant of the one you're already using:  generate a *permutation* at
random, in one gulp.  Then check whether it meets the restrictions.  If it
doesn't, try again.  Then each solution clearly has an equal chance of
getting picked.

Deep flaw #1:  As explained here a few months ago, Python's random-number
generator isn't capable of generating more than about 1E12 distinct
permutations, so even for a list of size 20 it's incapable of generating
"almost all" of the possible permutations.

Deep flaw #2:  As with your current approach, it will run "forever" if there
is no solution.  A variant of the birthday paradox kicks in here too:  even
if all permutations could be generated at random, "at random" means you get
duplicates.  The expected number of permutations you have to generate at
random before you see each of the N! possibilites at least once is about N!
* ln(N!) (ln == log to the base e).  Cut off the search earlier than that,
and chances are better than half that you never saw at least one permutation
(which may have been a solution).

2) The approach below:  generate all possible solutions, then simply apply
random.choice to the list of results.  For your size of problem, it will run
quickly, and in any case runs faster the more restrictions there are.

Deep flaw:  For larger problems, it's utterly impractical (too many
solutions).

3) Go back to the theory:  at each step, analyze the board according to the
method sketched yesterday, and put a rook on that step's distinguised square
with probability given by the number of solutions that have a rook on that
square divided by the total number of solutions (i.e., with or without a
rook on that square).  In the example above, this means that if your first
choice is wrt the Carel->Martijn square, you would first compute that there
are 2 solutions with a rook there and 1 without, and so you should put the
rook there with probability 2/3.  That's the way to do this properly without
backtracking -- and if you think it's worth the effort, your winter is far
too long <wink -- but it is a difficult programming problem!>.

van-lotje-getikt-ly y'rs  - tim

class Board:
    def __init__(self):
        # list of names
        self.names = []

        # map name to index in self.names
        self.name2i = {}

        # set of forbidden positions (pairs of ints)
        self.forbidden = {}

    def register_name(self, name):
        if not self.name2i.has_key(name):
            self.name2i[name] = len(self.names)
            self.names.append(name)

    # Remember that name x can't map to name y.
    def place_cheese(self, x, y):
        for name in x, y:
            if not self.name2i.has_key(name):
                raise ValueError("I don't know " + name + "!")
        self.forbidden[self.name2i[x], self.name2i[y]] = 1

    # Put cheese along the main diagonal (that is, don't let
    # anyone map to themself).
    def diagonalize(self):
        for i in range(len(self.names)):
            self.forbidden[i, i] = 1

    # Return list of all solutions.  Each solution is a list
    # of name pairs.
    def solve(self):
        self.board = board = []
        n = len(self.names)
        rangen = range(n)

        # Set board[i] to a list of its non-forbidden columns.
        cheesy = self.forbidden.has_key
        for row in rangen:
            goodcols = []
            for col in rangen:
                if not cheesy((row, col)):
                    goodcols.append(col)
            board.append(goodcols)
        results = []

        # Build set of all column indices
        columns = {}
        for i in rangen:
            columns[i] = 1

        # Find all solutions.
        _solve(board, 0, columns, [None]*n, results)

        # Translate results back to name pairs.
        names = self.names
        for i in xrange(len(results)):
            asolution = results[i]
            for row in rangen:
                col = asolution[row]
                asolution[row] = (names[row], names[col])
        return results

def _solve(board, row, cols, current, results):
    # board is the current board.
    # row is the index of the row to fill on this call,
    # and increases by one for each recursive call.
    # cols is the set of column indices still available.
    # current is a list whose i'th entry is the index of the
    # column to which row i maps (current solution so far).
    # results is the list of all complete solutions so far.
    if not cols:
        results.append(current[:])
        return
    # Try all ways of placing in this row (the set of columns not
    # forbidden to this row intersected with the set of columns not
    # already filled by shallower calls).
    still_available = cols.has_key
    for col in board[row]:
        if still_available(col):
            current[row] = col
            del cols[col]
            _solve(board, row+1, cols, current, results)
            cols[col] = 1