#!/usr/bin/env python
##
## Name:     perm.py
## Purpose:  Various computations related to permutation.
## Author:   M. J. Fromberger <http://www.dartmouth.edu/~sting/>
## Info:     $Id: Permutation.py 441 2007-10-19 14:59:52Z sting $
##

from operator import add, itemgetter
import random

class Permutation (object):
    """This class allows the manipulation of permutations as
    computational objects.  The basic operation on a permutation
    is multiplication, which applies the permutation on the RHS
    to the permutation or sequence on the LHS and returns the
    result.
    """
    def __init__(self, seq):
        """Construct a new permutation from a list of indices.  The
        indices must cover the range (0,..len(seq)).
        """
        self._perm = tuple(seq)
    
    def __len__(self):
        """Return the number of elements of the permutation."""
        return len(self._perm)
    
    def __getitem__(self, itm):
        """Read elements or slices of the permutation."""
        return self._perm[itm]

    def __iter__(self):
        """Iterate over the indices of the permutation."""
        return iter(self._perm)
    
    def __repr__(self):
        return 'Permutation(%s)' % (self._perm,)
    
    def __mul__(self, other):
        """Multiply (compose) this permutation with another permutation.
        The dimensions of the two arguments must match.
        """
        if len(other) <> len(self):
            raise ValueError("Arguments have incompatible dimensions")
        
        try:
            return Permutation(self[i] for i in other)
        except TypeError:
            raise ValueError("Incorrect type for permutation index")
    
    def __rmul__(self, other):
        """Right-multiply another sequence type by this permutation,
        assuming compatible dimensions.
        """
        if len(other) <> len(self):
            raise ValueError("Arguments have incompatible dimensions")
        
        if isinstance(other, str):
            return ''.join(other[i] for i in self)
        else:
            return type(other)(other[i] for i in self)

    def __add__(self, other):
        """Add this permutation to either another permutation or to an
        index.  The effect is to add the indices of the permutation,
        modulo the number of total permutations.  If the right-hand
        argument is a permutation, it must have compatible dimensions.
        """
        my_idx = self.index()
        try:
            other_idx = other.index()
            if len(self) <> len(other):
                raise ValueError("Arguments have incompatible dimensions")
        except AttributeError:
            try:
                other_idx = int(other)
            except ValueError:
                raise TypeError("Right-hand argument must be permutation"
                                " or integer")
        
        new_idx = my_idx + other_idx
        try:
            return Permutation.from_index(len(self), new_idx)
        except IndexError:
            return Permutation.from_index(len(self),
                                          new_idx % factorial(len(self)))
    
    def __pow__(self, z):
        """Raise a permutation to a power by repeated composition."""
        if not isinstance(z, (int, long)) or z < 0:
            raise TypeError("power must be a non-negative integer")

        # This is basically the square-and-multiply algorithm that is
        # used for ordinary exponentiation, only using composition of
        # permutations instead
        
        out = Permutation.identity(len(self))
        cur = self
        while z <> 0:
            if z % 2 == 1:
                out *= cur

            cur *= cur
            z /= 2

        return out
    
    def inverse(self):
        """Return the inverse of the permutation."""
        inv = [0] * len(self)
        for i in xrange(len(self)):
            inv[self[i]] = i
        
        return Permutation(inv)
    
    def reversed(self):
        """Return the reverse of the permutation."""
        return Permutation(reversed(self))
    
    def invert(self):
        """Invert the permutation in place."""
        inv = [0] * len(self)
        for i in xrange(len(self)):
            inv[self[i]] = i
        
        self._perm = tuple(inv)
    
    def reverse(self):
        """Reverse the permutation in place."""
        self._perm = tuple(reversed(self._perm))
        return self
    
    def cycles(self):
        """Return a list of permutation indices in cycle normal form:
        Each cycle begins with its least index, and all cycles are
        ordered decreasing by their primary indices.
        """
        def tindex(tup, elt):  # Since tuples have no .index().
            for pos in xrange(len(tup)):
                if tup[pos] == elt:
                    return pos
            else:
                return -1
        
        P = self._perm
        not_done = list(P)
        cycles = []
        
        while not_done:
            pos = tindex(P, not_done[0])
            cycle = [ pos ]
            not_done.remove(pos)
            
            cursor = P[pos]
            while cursor <> pos:
                cycle.append(cursor)
                not_done.remove(cursor)
                cursor = P[cursor]
        
            cycles.append(cycle)
        
        return reduce(add, sorted(cycles, key = itemgetter(0),
                                  reverse = True))

    def index(self):
        """Return the index of this permutation in the canonical
        factoradic mapping of permutations to indices in (0,..n!).
        """
        tmp = range(len(self))
        idx = [None] * len(self)

        for pos, elt in enumerate(self):
            idx[pos] = tmp.index(elt)
            tmp.pop(idx[pos])

        return fact2int(idx)
    
    @staticmethod
    def identity(n):
        """Generate an identity permutation of (0,..n)."""
        return Permutation(xrange(n))

    @staticmethod
    def random(n):
        """Generate a random permutation of (0,..n)."""
        tmp = range(n)
        random.shuffle(tmp)
        return Permutation(tmp)
    
    @staticmethod
    def from_index(n, z):
        """Construct a permutation of (0,..n) from an integer index in
        the range (0,..n!).  Raises IndexError if z >= n!
        """
        tmp = range(n)
        out = [None] * n
        idx = int2fact(z)
        idx = [0] * (n - len(idx)) + idx
        
        for pos, elt in enumerate(idx):
            out[pos] = tmp[elt]
            tmp.pop(elt)
        
        return Permutation(out)
    
    @staticmethod
    def from_cycles(cyc):
        """Construct a permutation from a list in cycle normal form."""
        breaks = [ p for p in xrange(len(cyc))
                   if p > 0 and not [ x for x in cyc[:p]
                                      if x < cyc[p] ] ]
        
        res = [None] * len(cyc)
        last = 0
        
        for b in breaks + [len(cyc)]:
            cur = cyc[last:b]
            last = b
            
            save = prev = cur[0]
            for i in cur[1:]:
                res[prev] = i
                prev = i
            
            res[prev] = save
        
        return Permutation(res)

def invert_perm_inplace(P):
    """Given a permutation P of (0,..n), specified as a list, compute
    the inverse permutation P' in place, replacing P.  The result P is
    also returned.
    """
    mark = [False] * len(P)
    
    for i in xrange(len(P)):
        if not mark[i]:
            j = P[i]
            save = i
            
            while not mark[j]:
                mark[j] = True
                next = P[j]
                P[j] = save
                save = j
                j = next
    
    return P

def generate_moves(cycles):
    """Given a permutation of (0,..n) as a list of disjoint cycles,
    return a minimal list of the "moves" necessary to apply that
    permutation to a sequence.  The moves are human-readable strings
    of the form
    
    move <src>, <dst>

    ... where <src> and <dst> are either locations in the sequence
    being permuted, or the word "spare" to indicate a location that is
    separate."""
    moves = []
    
    for cycle in cycles:
        if len(cycle) > 1:
            moves.append('move %d, spare' % cycle[0])
            
            prev = cycle[0]
            for elt in reversed(cycle[1:]):
                moves.append('move %d, %d' % (elt, prev))
                prev = elt
            
            moves.append('move spare, %d' % cycle[1])
    
    return moves

def permutations(seq):
    """Returns a generator that yields all the permutations of seq in
    lexicographic order.  The order of elements in seq determines the
    basic ordering of permutations.

    The generator does not copy its input, so modifying the source
    while the generator is idle will change the output.
    """

    # In-place permutation algorithm from E. W. Dikjstra, "A
    # Discipline of Programming" (1976).
    
    px = range(len(seq) + 1)  # permutation indices
    px[-1] = -1               # sentinel
    
    while True:
        yield list(seq[p] for p in px[:-1])
        
        for i in reversed(xrange(len(seq))):
            if px[i] < px[i + 1]:
                break
        else:
            return # no more permutations

        for j in reversed(xrange(len(seq))):
            if px[j] > px[i]:
                break

        px[i], px[j] = px[j], px[i]

        for j in reversed(xrange(len(seq))):
            i += 1
            if j <= i:
                break
            px[i], px[j] = px[j], px[i]

def count_breaks(n, k):
    """Returns the number of distinct ways to partition a sequence of
    length n into k + 1 groups of contiguous elements.  For instance,
    if n = 4 and k = 2, there are three ways: a:b:cd, a:bc:d, and
    ab:c:d."""
    
    def cb(n, k, c):
        assert n > 0
        if (n, k) in c:
            return c[n, k]
        
        if k == 0 or k == n - 1:
            return 1
        elif k == 1:
            return n - 1
        else:
            c[n, k] = sum([ sum([ cb(p, i, c) *
                                  cb(n - p, k - i, c)
                                  for i in xrange(min(k, p)) ])
                            for p in xrange(1, n) ])
            return c[n, k]
    
    return cb(n, k, {})

def int2fact(z):
    """Compute the factoradic representation of abs(z) for an integer
    z, returning the digits in an array of decreasing order of
    significance.
    """
    if z == 0:
        return [0]
    elif z < 0:
        z = -z

    out = []
    rdx = 1
    while z <> 0:
        out.append(z % rdx)
        z /= rdx
        rdx += 1

    out.reverse()
    return out

def fact2int(fr):
    """Convert a factoradic representation of an integer back into the
    integer it represents.  The input should be a list of digit values
    in decreasing order of significance.
    """
    out = 0
    pval = 1
    for pos, digit in enumerate(reversed(fr)):
        out += pval * digit
        pval *= pos + 1
    
    return out

def factorial(n):
    """Compute the integer factorial."""
    s = 1
    while n > 0:
        s *= n
        n -= 1

    return s

__all__ = ("Permutation", "int2fact", "fact2int", "permutations")

# Here there be dragons