#coding:utf-8
from sage.matroids.advanced import BasisMatroid
import copy
from itertools import chain,combinations

def powerset(s):
    """ (s)

        return a chain-iterator for the powerset of s """
    return chain.from_iterable(combinations(s, r) for r in xrange(len(s)+1))

def powersetSubsetIndex(S, X):
    """ (S, X)

        S __ subset of X
        X __ universe (list or implicitly sorted)

        returns the iteration index where the subset S occurs in the powerset
        iterator
    """
    if type(X) != list:
        X = sorted(X)
    idx0 = 0
    for i in range(len(S)):
        idx0 += binomial(len(X),i)
    idx00 = idx0
    idxs = sorted((X.index(s) for s in S))
    x0 = 0
    s0 = 0
    for i in idxs:
        for j in range(x0,i):
            # add number of choices where j is fixed at position s0;
            # consequently, the subsequent numbers must be bigger than j.
            idx0 += binomial(len(X) - j -1, len(S) - s0 - 1)
        s0 += 1
        x0 = i + 1
    return idx0


def allFlats(M):
    """ (M)

        return a chain-iterator of all flats of M, with ascending rank """
    return chain.from_iterable(M.flats(r) for r in xrange(0,M.rank()+1))

def leqWrtFlats(l,r,flats):
    """ (l,r, flats)

        l,r __ operands

        flats __ set/family of frozensets

        compares whether l is a subset of r, and if l is a proper subset,
                whether l in flats;
                returns True if l <=(flats) r.
    """
    l = frozenset(l)
    r = frozenset(r)
    if l == r:
        return True
    if not l.issubset(r):
        return False
    if not l in flats:
        return False
    return True

def downsetWrtFlats(r,flats):
    """ (r,flats)

       r __ operand

       flats __ set/family of frozensets

       returns the subsets l of r, for which l <=(flats) r holds.
    """
    r = frozenset(r)
    return frozenset([r] + [l for l in flats if l.issubset(r)])


def alphaPosetLeq(M):
    """ (M)

        M __ matroid

        returns a function that compares whether l <= r with respect to
                the alpha_M poset, and then returns either True of False
    """
    return lambda l,r,f=frozenset(allFlats(M)): leqWrtFlats(l,r,f)

def alphaPosetDownsets(M):
    """ (M)

        M __ matroid

        returns a function that gives the downsets of the alpha_M poset
    """
    return lambda r,f=frozenset(allFlats(M)): downsetWrtFlats(r,f)


def moebiusMatrixForAlphaPoset(M):
    """ (M)

        M __ matroid

        returns the Moebius-function of the alpha_M poset.
    """
    G = sorted(M.groundset())
    nG = len(G)
    mu = Matrix(2**nG,2**nG,sparse=True)
    idx1 = -1
    leq = alphaPosetLeq(M)
    downset = alphaPosetDownsets(M)
    for L in powerset(G):
        idx1 += 1
        idx2 = -1
        for R in powerset(G):
            idx2 += 1
            if L == R:
                mu[idx1,idx2] = 1
            elif leq(L,R):
                s = 0
                for P in downset(R):
                    if P == R:
                        continue
                    s += mu[idx1,powersetSubsetIndex(P,G)]
                mu[idx1,idx2] = -s
    return mu

def alphaVector(M):
    """ (M)

        M __ matroid

        returns the alpha vector of M
    """

    G = sorted(M.groundset())
    nG = len(G)
    alpha = Matrix(1,2**nG,sparse=True)
    idx1 = -1
    leq = alphaPosetLeq(M)
    downset = alphaPosetDownsets(M)
    for L in powerset(G):
        idx1 += 1
        aval = len(L) - M.rank(L)
        for R in downset(L):
            if R == L:
                continue
            aval -= alpha[0,powersetSubsetIndex(R,G)]
        if aval:
            alpha[0,idx1] = aval
    return alpha

def calculateAlphaOfExtension(e, C, M, alpha=None, mu=None):
    """ (e, C, M, alpha=None, mu=None)


            e __ name of the new element
            C __ modular cut of the new element (whole cut, not just minimal elts!)
            M __ matroid

            alpha __ alpha vector of M
            mu  __ Möbius function of the alpha_M-poset

        returns the alpha vector of M+e that corresponds to C
    """
    G0 = sorted(M.groundset())
    G1 = sorted(G0+[e])
    if alpha is None:
        alpha = alphaVector(M)
    if mu is None:
        mu = moebiusMatrixForAlphaPoset(M)

    alphaE = Matrix(1,2**len(G1),sparse=True)
    idxE = -1
    deltaAlphaE = {}

    for X in powerset(G1):
        X = frozenset(X)
        idxE += 1
        val = 0
        if e in X:
            X0 = X.difference([e])
            clX0 = M.closure(X0)
            if not clX0 in C:
                if clX0 == X0:
                    val = 0
                else:
                    val = alpha[0,powersetSubsetIndex(X0,G0)]
            else:
                d = 1
                for Z in C:
                    if Z.issubset(X0):
                        if Z == X0:
                            continue
                        #print setStr(Z), "<=", setStr(X0)
                        d -= deltaAlphaE[Z]
                deltaAlphaE[X0] = d
                val = alpha[0,powersetSubsetIndex(X0,G0)] + d
        else:
            s = 0
            for Y in powerset(X):
                Y = frozenset(Y)
                if Y == X:
                    continue
                nuY = len(Y) - M.rank(Y)
                for Z in C:
                    Z = frozenset(Z)
                    if not Y.issubset(Z):
                        continue
                    if not Z.issubset(X):
                        continue
                    if Z == X:
                        continue
                    s += nuY * mu[powersetSubsetIndex(Y,G0),powersetSubsetIndex(Z,G0)]
            val = alpha[0,powersetSubsetIndex(X,G0)] + s # Lemma 2.4.16
        if val:
            alphaE[0,idxE] = val
    return alphaE

def getAllModularCuts(M):
    """ (M)

        M __ matroid

        returns all modular cuts of M, except the empty modular cut.
        (SAGE does not allow an extension with a coloop through M.extension !!)
    """
    H = frozenset(M.hyperplanes())
    L = [frozenset(l) for l in M.linear_subclasses()]
    Ms = [set() for l in L]
    for F in allFlats(M):
        for m,l in zip(Ms,L):
            contained = True
            for h in H:
                if F.issubset(h):
                    if not h in l:
                        contained = False
                        break
            if contained:
                m.add(F)

    return tuple(( frozenset(x) for x in Ms ))

def hasMinorWith(self, N, E0, certificate=False):
    """
    (M, N, E0, bint certificate=False)
    Test if matroid M has the specified minor, that contains the set E0
    and optionally return frozensets ``X`` and ``Y`` so that ``N`` is isomorphic to ``self.minor(X, Y)``.

    INPUT:

    - ``N`` -- An instance of a ``Matroid`` object,
    - ``certificate`` -- Boolean (Defalt: ``False``) If ``True``, returns
      ``True, (X, Y, dic) where ``N`` is isomorphic to ``self.minor(X, Y)``,
      and ``dic`` is an isomorphism between ``N`` and ``self.minor(X, Y)``.

    OUTPUT:

    boolean or tuple.

    EXAMPLES::

        sage: M = matroids.named_matroids.Vamos()
        sage: M._has_minor(matroids.Whirl(3))
        False
        sage: M._has_minor(matroids.Uniform(2, 4))
        True
        sage: M._has_minor(matroids.Uniform(2, 4), certificate=True)
        (True, (frozenset({'a', 'c'}), frozenset({'b', 'e'}),
            {0: 'h', 1: 'd', 2: 'g', 3: 'f'}))

    .. TODO::

        This important method can (and should) be optimized considerably.
        See [Hli2006]_ p.1219 for hints to that end.
    """
    if self is N:
        if certificate:
           return True, (frozenset(), frozenset(), {x: x for x in self.groundset()})
        return True
    rd = self.full_rank() - N.full_rank()
    cd = self.full_corank() - N.full_corank()
    if rd < 0 or cd < 0:
        if certificate:
            return False, None
        return False
    YY = []
    for Y in self.dual().independent_r_sets(cd):
        if not Y.intersection(E0):
            YY.append(Y)
    for X in self.independent_r_sets(rd):
        if X.intersection(E0):
            continue
        for Y in YY:
            if X.isdisjoint(Y):
                if N._is_isomorphic(self._minor(contractions=X, deletions=Y)):
                    if certificate:
                        return True, (X, Y, N._isomorphism(self._minor(contractions=X, deletions=Y)))
                    return True
    if certificate:
        return False, None
    return False


def canonicalMatroidLabel(M):
    """ (M)

        M __ matroid

        returns the canonical label of the matroid M;
               which is a tuple
    """
    IS = IncidenceStructure(M.bases())
    phi = IS.canonical_label()

    return (frozenset((frozenset((phi[x] for x in B)) for B in M.bases())), len(M.groundset()))

class AlphaMatroid(object):
    def __init__(self, M):
        """ (M)
            M __ matroid

            creates a decorated matroid object with added functionality for
               gammoid and strict gammoid backtracking
        """
        self.M = M
        self.E = tuple(sorted(M.groundset()))
        self._alpha = None
        self._mu = None
        self._modcuts = None
        for name in set(dir(M))-set(dir(self)):
            setattr(self, name, getattr(M,name))

    def __cmp__(self,R):
        if type(R) is AlphaMatroid:
            return cmp(self.M,R.M)
        else:
            return cmp(self.M,R)

    def modular_cuts(self):
        """ ()

            returns a tuple of all modular cuts of this matroid
        """
        if self._modcuts is None:
            self._modcuts = getAllModularCuts(self.M)
        return self._modcuts

    def alpha(self):
        """ ()

            returns the alpha vector of M, which is sorted in the power-set
                    enumeration order of powerset(self.E)
        """
        if self._alpha is None:
            self._alpha = alphaVector(self.M)
        return self._alpha

    def mu(self):
        """ ()

            returns the Moebius-function of the alpha_M poset of this matroid.
        """
        if self._mu is None:
            self._mu = moebiusMatrixForAlphaPoset(self)
        return self._mu

    def is_strict_gammoid(self):
        """ ()

            returns True, if M is a strict gammoid
        """
        return min(self.alpha()[0]) >= 0

def decorateAlphaMatroid(M):
    """ (M)
        M __ matroid

        returns an AlphaMatroid(M) or M, depending whether M has been decorated
        or not.
    """
    if type(M) == AlphaMatroid:
        return M
    return AlphaMatroid(M)

class Rank3Decider(object):
    def __init__(self):
        self.NonGammoids = []
        self.statistics = [0,0,0,1]
        self.Gammoids = [AlphaMatroid(BasisMatroid(bases=[[0,1,2]],groundset=[0,1,2]))]
        self.nMax = 3

    def __call__(self,M, E0=None):
        """ (M, E0=None)

         M __ matroid

         E0 __ only test minors of M that contain all elements from E0

        tests whether M has an excluded rank-3 minor with at most self.nMax
        elemests. If such one is found, it returns an isomorphic minor, otherwise
        it returns False
        """
        if E0:
            for G in self.NonGammoids:
                if hasMinorWith(M,G.M,E0):
                    return G
        else:
            for G in self.NonGammoids:
                if M.has_minor(G.M):
                    return G
        return False

    def prepare(self, n):
        """ (n)

            n __ nbr of elements

            prepares (up to isomorphism) all excluded minors for Gammoids
            with rank 3 and up to (n) elements.
        """
        if n <= self.nMax:
            return
        for i in range(self.nMax+1,n+1):
            NextGammoids = []
            Ei = frozenset([i])
            found = set()
            for G in self.Gammoids:
                for C in G.modular_cuts():
                    M = G.extension(element=i,subsets=C)
                    L = canonicalMatroidLabel(M)
                    if not L in found:
                        found.add(L)
                        if self(M,Ei) == False: #has no known excluded minor
                            M = AlphaMatroid(M)
                            M._alpha = calculateAlphaOfExtension(i,C,G,G.alpha(),G.mu())
                            if M.is_strict_gammoid():
                                NextGammoids.append(M)
                            else:
                                self.NonGammoids.append(M)
            self.Gammoids = NextGammoids
            self.statistics.append(len(NextGammoids))
            self.nMax = i

    def __repr__(self):
        return "Collection of "+str(len(self.NonGammoids))+" excluded rank-3 minors for gammoids up to "+str(self.nMax)+" elements\n" +\
        "  and "+str(len(self.Gammoids))+" non-isomorphic rank-3 gammoids with "+str(self.nMax)+" elements."

R3D=Rank3Decider()

def a(n):
    """ returns the n-th element of the sequence A300985 """
    R3D.prepare(n)
    return R3D.statistics[n]