r"""
Relational Structures: search for finite algebras and relations

AUTHOR:
    -- Peter Jipsen (2007-02-25): initial version

  TODO: convert output of the search routines to RelStructure objects.

  NOTE: currently this file is Python code with no SAGE dependencies
        (except the use of sloane_find to test the searches)

TUTORIAL:

    I. The Basics
    
        1. Relational Structure Format
        
            A. The SAGE RelStructure Class
            
            Finite (universal) algebras and relational structures in SAGE are
            implemented as dictionaries of operation tables and relations.
            The underlying set for a relational structure of size n is {0,1,...,n-1}.
            A key "labels" (if defined) contains a dictionary with a
            printable label for each element.
            A key "algtype" contains a dictionary with the arity of each operation.
            A key "reltype" contains a dictionary with the arity of each relation.
            Operations and relations are accessed by keys that represent their names.
            E.g. binary operations   "+", "*", "//", "\\", "^^", "vv"
                 unary operations    "~", "-", "inv"
                 constant operations "e", "top", "bot"
                 binary relations    "<=", "-<" (covering)
            
        2. Search routines

            The current file contains code to search for nonisomorphic groupoids
            and binary relations.
            A groupoid is given by an nxn (python) matrix with entires 0,...,n-1
            A binary relation is given by an nxn (python) 0-1-matrix
            A partial groupoid or binary relation may also have -1 = undefined as entry
            Permutations of {0,...,n-1} are represented by lists of length n.
            Searches are done in a simple minded way by backtracking.

            Examples: Find all nonisomorphic algebras of a given size.

            find_groupoids(3)
            find_semigroups(5)             # associative groupoids
            find_commutative_semigroups(5)
            find_bands(6)                  # idempotent semigroups
            find_monoids(6)                # semigroups with identity
            find_commutative_monoids(7) 
            find_semilattices(7)           # commutative bands
            find_lattices(8)               # semilattices with identity

            find_algebras_test(True)       # test all find_ methods and compare
                                           # lengths of lists with sloane_find

              Find all nonisomorphic binary relations of a given size.
              
            find_binary_relations(5)
            find_simple_digraphs(6)  #binary relations without loops
            find_simple_graphs(7)    #simple digraphs with arrows in both directions

            find_relations_test(True)      # test all find_ methods and compare
                                           # lengths of lists with sloane_find

"""

#*****************************************************************************
#           Copyright (C) 2007 Peter Jipsen <jipsen@chapman.edu>
#
# Distributed  under  the  terms  of  the  GNU  General  Public  License (GPL)
#                         http://www.gnu.org/licenses/
#*****************************************************************************

from sage.structure.sage_object import SageObject
import datetime

class RelStructure(SageObject):

    def __len__(self):
        return len(self[0])

    def __str__(self):
        if self.name != None:
            return self.name
        else: return repr(self)

    def _repr_(self):
        return "{"+",\n".join([str(k)+":"+\
                               (opstr(v) if isinstance(v,list) else str(v))\
                               for k,v in vars(self).iteritems()])+"}"

    def _latex_(self):
        return repr(self)

    ### General properties

    def order(self):
        """
        Returns the number of elements.
        """
        return len(self[0])

    def size(self):
        """
        Returns the number of elements.
        """
        return len(self[0])

def permutations(m,n):
    """
    return generator for all permutations of {m,...,n}
    """
    p = [m+i for i in range(n-m)]
    yield p
    n = len(p)
    j = 1
    while j>=0:
        q = range(n)
        j = n-2
        while j>=0 and p[j]>p[j+1]: j = j-1
        if j>=0:
            for k in range(j): q[k] = p[k]
            k = n-1
            while p[j]>p[k]: k = k-1
            q[j] = p[k]
            i = n-1
            while i>j:
                q[i] = p[j+n-i]
                i = i-1
            q[j+n-k]=p[j]
            p = q
            yield q

def inverse_permutation(p):
    q = range(len(p))
    for i in range(len(p)):
        q[p[i]]=i
    return q

def check_permutations(m):
    """
    Return False if every completion of some isomorphic copy q(m) of m
    will be lexicographically smaller than m
    Return True otherwise
    """
    global perms, invperms, pflags, pflags_list
    n = len(m)
    for k in range(len(perms)):
        if pflags[k]:
            q = perms[k]
            qi = invperms[k]
            equal = True
            for i in range(n):
                if equal: #equal means q[m[qi[i]][qi[j]]] == m[i][j] != -1
                    for j in range(n):
                        if equal:
                            mqi = m[qi[i]][qi[j]]
                            mij = m[i][j]
                            equal = (mqi>-1 and mij==q[mqi])
                            if mqi>-1 and q[mqi]<mij: return False
                            if mqi>-1 and mij>-1 and q[mqi]>mij: #get rid of q
                                pflags[k] = False
                                pflags_list += [k]
    return True

def check_associativity(m):
    """
    Return False if the partial operation table  m  violates associativity
    Return True otherwise
    """
    global entries_list
    n = len(m)
    for x in range(n):
        for y in range(n):
            xy = m[x][y]
            if xy>-1:
                for z in range(n):
                    xyz = m[xy][z];
                    if xyz>-1:
                        if m[y][z]>-1:
                            if m[x][m[y][z]]>-1:
                                if xyz!=m[x][m[y][z]]: return False
                            else:
                                m[x][m[y][z]] = xyz
                                entries_list += [(x,m[y][z])]
                    else:
                        if m[y][z]>-1 and m[x][m[y][z]]>-1:
                            m[xy][z] = m[x][m[y][z]]
                            entries_list += [(xy,z)]
    return True

def check_selfdual(m):
    """
    Return False if the partial meet table  m  is not selfdual wrt  sd  map
    Return True otherwise
    sdl is a global list s.t. sdl[i][sdl[i][x]]=x
    selfdual means m[x][y]=x => m[s[x]][s[y]]=s[y] for some s in sdl
    """
    global entries_list,sd
    n = len(m)
    for x in range(n):
        for y in range(n):
            xy = m[x][y]
            if xy==x:
                if x>y: return False  #force linear extension
                sxsy = m[sd[x]][sd[y]]
                if sxsy>-1:
                    if sxsy!=sd[y]: return False
                else:
                    m[sd[x]][sd[y]] = sd[y]
                    entries_list += [(sd[x],sd[y])]
    return True

def complete_operation(m,i,j,associative=False,commutative=False,selfdual=False):
    """
    find next i,j where alg[i][j] = -1 = undefined; for each val=0..n-1
    set alg[i][j] = val, check permutations
    if ok, complete_operation(m,i,j+1) then restore and return
    """
    global algebra_list,entries_list,pflags,pflags_list
    ok = True
    n = len(m)
    while ok and i<n:
        while j<n and m[i][j]>-1: j = j+1
        if j>=n:
            j = 0 
            i = i+1
        else: ok = False
    if ok:
        if not selfdual or check_selfdual(m):
            algebra_list.append([m[i][:] for i in range(n)])
    else:
        for val in range(n):
            ok = True
            el = len(entries_list)
            m[i][j] = val
            entries_list += [(i,j)]
            if commutative:
                m[j][i] = val
                entries_list += [(j,i)]
            if selfdual: ok = check_selfdual(m)
            if ok and associative: ok = check_associativity(m)
            pl = len(pflags_list)
            if ok: ok = check_permutations(m)
            if ok: complete_operation(m,i,j+1,associative,commutative,selfdual)
            while len(entries_list)>el:
                e = entries_list.pop()
                m[e[0]][e[1]] = -1
            while len(pflags_list)>pl:
                p = pflags_list.pop()
                pflags[p] = True

def find_groupoids(n,associative=False,commutative=False,idempotent=False,identity=False):
    """
    Return list of nonisomorphic groupoids of size  n  as nxn matrices
    """
    global algebra_list,perms,invperms,entries_list,pflags,pflags_list
    algebra_list=[]
    entries_list=[]
    m = [[-1 for x in range(n)] for x in range(n)]
    if idempotent:
        for x in range(n):
            m[x][x] = x
    if identity:
        perms = [[0]+p for p in permutations(1,n)]
        for x in range(n):
            m[0][x] = x
            m[x][0] = x
    elif associative and commutative and idempotent: 
        perms = [[0]+p for p in permutations(1,n)] #requires lex order
        for x in range(n):
            m[0][x] = 0
            m[x][0] = 0
    else:
        perms = [p for p in permutations(0,n)]
    invperms = [inverse_permutation(p) for p in perms]
    pflags = [True for i in range(len(perms))]
    pflags_list = []
    complete_operation(m,0,0,associative,commutative)
    return algebra_list

def find_commutative_groupoids(n):
    """
    Return list of nonisomorphic commutative groupoids of size  n  as nxn matrices
    """
    return find_groupoids(n,commutative=True)

def find_semigroups(n):
    """
    Return list of nonisomorphic semigroups of size  n  as nxn matrices
    """
    return find_groupoids(n,associative=True)

def find_commutative_semigroups(n):
    """
    Return list of nonisomorphic commutative semigroups of size  n  as nxn matrices
    """
    return find_groupoids(n,associative=True,commutative=True)

def find_bands(n):
    """
    Return list of nonisomorphic bands of size  n  as nxn matrices
    """
    return find_groupoids(n,associative=True,idempotent=True)

def find_semilattices(n):
    """
    Return list of nonisomorphic semilattices of size  n  as nxn matrices
    0 is absorbtive so best to consider these algebras as meet-semilattices
    """
    return find_groupoids(n,associative=True,commutative=True,idempotent=True)

def find_lattices(n):
    """
    Return list of nonisomorphic lattices of size  n  as nxn matrices
    """
    if n==0: return [[]]
    alg_list = find_semilattices(n-1)
    alg_list = [[m[i]+[i] for i in range(n-1)]+[range(n)] for m in alg_list]
    return [Poset(meet2uc(m)) for m in alg_list]

def sdmap(p,n):
    k = len(p)
    return [n-1-p[k-1-i] for i in range(k)]

def find_selfdual_lattices(n):
    """
    Return list of nonisomorphic selfdual lattices of size  n  as nxn matrices
    """
    global algebra_list,perms,invperms,entries_list,pflags,pflags_list,sd
    algebra_list=[]
    entries_list=[]
    associative = True
    commutative = True
    m = [[-1 for x in range(n)] for x in range(n)]
    sd = range(n)
    pairs = int(n)/2
    for x in range(pairs):
        sd[x] = n-1-x
        sd[n-1-x] = x
    perms = [[0]+p for p in permutations(1,pairs)]
    perms =[p+q+sdmap(p,n) for p in perms for q in permutations(pairs,n-pairs)]
    for x in range(n):
        m[x][x] = x
        m[n-1][x] = x
        m[x][n-1] = x
        m[0][x] = 0
        m[x][0] = 0
    invperms = [inverse_permutation(p) for p in perms]
    pflags = [True for i in range(len(perms))]
    pflags_list = []
    complete_operation(m,0,0,associative,commutative,True)
    return [Poset(meet2uc(m)) for m in algebra_list]

def find_monoids(n):
    """
    Return list of nonisomorphic monoids of size  n  as nxn matrices
    """
    return find_groupoids(n,associative=True,identity=True)

def find_commutative_monoids(n):
    """
    Return list of nonisomorphic commutative monoids of size  n  as nxn matrices
    """
    return find_groupoids(n,associative=True,commutative=True,identity=True)

def find_algebras_test(check_sloane=False):
    if check_sloane: n = 7
    else: n = 6
    tl = [
        [len(find_groupoids(i)) for i in range(1,4)],
        [len(find_semigroups(i)) for i in range(1,5)],  #associative groupoids
        [len(find_commutative_semigroups(i)) for i in range(1,5)],
        [len(find_bands(i)) for i in range(1,6)],       #idempotent semigroups
        [len(find_monoids(i)) for i in range(1,n)],     #semigroups with identity
        [len(find_commutative_monoids(i)) for i in range(1,7)],
        [len(find_semilattices(i)) for i in range(1,7)],#commutative bands
        [len(find_lattices(i)) for i in range(1,7)]]    #semilattices with identity
    if not check_sloane: return tl
    sl = [[t]+sloane_find(t,1)[0] for t in tl]
    return [[s[0],s[3][:8],s[2]] for s in sl]

##############################################################################
# Search for linearly ordered structures

def is_ordered(m):
    for x in range(len(m)-1):
        for y in range(len(m)):
            if m[x+1][y]>-1 and m[x][y]>m[x+1][y]: return False
            if m[y][x+1]>-1 and m[y][x]>m[y][x+1]: return False
    return True

def is_associative(m):
    for x in range(len(m)):
        for y in range(len(m)):
            for z in range(len(m)):
                if m[m[x][y]][z]!=m[x][m[y][z]]: return False
    return True

def is_commutative(m):
    for x in range(len(m)):
        for y in range(len(m)):
            if m[x][y]!=m[y][x]: return False
    return True

def is_idempotent(m):
    for x in range(len(m)):
        if m[x][x]!=x: return False
    return True

def is_partial(m):
    for x in range(len(m)):
        for y in range(len(m)):
            if m[x][y]<0 or m[x][y]>=len(m): return True
    return False

def is_ordered_semigroup(m):
    return is_ordered(m) and is_associative(m)

def complete_ordered_operation(m,i,j,associative=False,commutative=False):
    """
    find next i,j where alg[i][j] = -1 = undefined; for each val=0..n-1
    set alg[i][j] = val, check axioms and permutations
    if ok, complete_ordered_operation(m,i,j+1) then restore and return
    """
    global entries_list
    ok = True
    n = len(m)
    while ok and i<n:
        while j<n and m[i][j]>-1: j = j+1
        if j>=n:
            j = 0 
            i = i+1
        else: ok = False
    if ok:
        algebra_list.append([m[i][:] for i in range(n)])
    else:
        ok = True
        k = 0
        if i>0: k = m[i-1][j]
        if j>0: k = max(k,m[i][j-1])
        for val in range(k,n):
            el = len(entries_list)
            m[i][j] = val
            entries_list += [(i,j)]
            if commutative:
                m[j][i] = val
                entries_list += [(j,i)]
            if associative: ok = check_associativity(m)
            if ok: ok = is_ordered(m)
            if ok: complete_ordered_operation(m,i,j+1,associative,commutative)
            while len(entries_list)>el:
                e = entries_list.pop()
                m[e[0]][e[1]] = -1

def find_ordered_groupoids(n,associative=False,commutative=False,idempotent=False,identity=False,id=None,zero=False):
    """
    Return list of linearly ordered groupoids of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    global algebra_list,entries_list
    algebra_list=[]
    entries_list=[]
    m = [[-1 for x in range(n)] for x in range(n)]
    if idempotent:
        for x in range(n):
            m[x][x] = x
    if zero:
        for x in range(n):
            m[x][0] = 0
            m[0][x] = 0
    if identity:
        if id==None:
            t = n
            if zero: s = 1
            else: s = 0
        else: s = id; t = id+1
        for e in range(s,t):
            for x in range(n):
                m[e][x] = x
                m[x][e] = x
            complete_ordered_operation(m,0,0,associative,commutative)
            for x in range(n):
                m[e][x] = -1
                m[x][e] = -1
    else:
        complete_ordered_operation(m,0,0,associative,commutative)
    return algebra_list

def find_o_semigroups(n):
    """
    Return list of linearly ordered semigroups of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    return find_ordered_groupoids(n,associative=True)

def find_commutative_o_semigroups(n):
    """
    Return list of linearly ordered commutative semigroups of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    return find_ordered_groupoids(n,associative=True,commutative=True)

def find_ordered_bands(n):
    """
    Return list of linearly ordered monoids of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    return find_ordered_groupoids(n,associative=True,idempotent=True)

def find_o_semilattices(n):
    """
    Return list of linearly ordered semilattices of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.

    Wall time: 103.34
    [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862] this looks like the Catalan numbers, no idea why
    """
    return find_ordered_groupoids(n,associative=True,commutative=True,idempotent=True)

def find_o_monoids(n,e=None):
    """
    Return list of linearly ordered monoids of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    return find_ordered_groupoids(n,associative=True,identity=True,id=e)

def find_commutative_o_monoids(n,e=None):
    """
    Return list of linearly ordered commutative monoids of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.

    Is this possible:
    'Counts the Bell numbers associated with the partition number [n].',
    [1, 1, 2, 6, 22, 92, 426, 2146, 11624]]]
    """
    return find_ordered_groupoids(n,associative=True,commutative=True,identity=True,id=e)

def find_o_monoids_with_zero(n,e=None):
    """
    Return list of linearly ordered monoids with zero of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    return find_ordered_groupoids(n,associative=True,identity=True,id=e,zero=True)

def find_commutative_o_monoids_with_zero(n,e=None):
    """
    Return list of linearly ordered commutative monoids with zero of size  n  as nxn matrices.
    The ordering is the natural one 0<1<2<...<n-1.
    The operation is compatible with the order.
    """
    return find_ordered_groupoids(n,associative=True,commutative=True,identity=True,id=e,zero=True)

def find_integral_o_monoids_with_zero(n):
    return find_ordered_monoids_with_zero(n,n-1)

def find_commutative_integral_o_monoids_with_zero(n):
    return find_ordered_commutative_monoids_with_zero(n,n-1)

def find_ordered_algebras_test():
    tl = [
        [len(find_ordered_groupoids(i)) for i in range(1,4)],
        [len(find_ordered_semigroups(i)) for i in range(1,5)],  #associative groupoids
        [len(find_ordered_commutative_semigroups(i)) for i in range(1,5)],
        [len(find_ordered_bands(i)) for i in range(1,5)],       #idempotent semigroups
        [len(find_ordered_monoids(i)) for i in range(1,5)],     #semigroups with identity
        [len(find_ordered_commutative_monoids(i)) for i in range(1,6)],
        [len(find_ordered_monoids_with_zero(i)) for i in range(1,5)],
        [len(find_ordered_commutative_monoids_with_zero(i)) for i in range(1,6)],
        [len(find_ordered_integral_monoids_with_zero(i)) for i in range(1,5)], #identity is top element
        [len(find_ordered_commutative_integral_monoids_with_zero(i)) for i in range(1,6)],
        [len(find_ordered_semilattices(i)) for i in range(1,6)]]#commutative bands
    return tl

############################################################################
# Search for partially ordered structures

#import poset
#poset.py jipsen@chapman.edu 2006-11-25

# a poset is an instance of the Poset class

# The base of a poset is always the set {0,1,...,size-1}.
# uc[x] is a sorted list of upper covers of x.
# Elements in uc[x] are (numerically) bigger than x.
# This means all posets are topologically sorted.
# 0 is always a minimal element.
# size-1 is always a maximal element.

# optional attributes:
# pos = list of [x,y] coordinates for element position in Hasse diagram
# labels = list of strings for element names

# The following command is useful to see the data in a poset:
# vars(p) shows the attributes of instance p

def permuted_binary_op(m,q):
    qi = inverse_permutation(q)
    return [[q[m[qi[x]][qi[y]]] for y in range(len(m))] for x in range(len(m))]

def leq2uc(le): # assumes le[x][y] => x <= y (topologically sorted)
    n = len(le)
    uc = []
    for a in range(n):
        S = []   # accumulate upper covers of a
        for x in range(a+1,n):
            if le[a][x]==1:
                y = len(S)-1
                while y>=0 and le[S[y]][x]==0:
                    y = y-1
                if y<0: S.append(x)
        uc.append(S)
    return uc

def meet2uc(m):
    n = len(m)
    uc = []
    for a in range(n):
        S = []
        for x in range(a+1,n):
            if m[a][x]==a:
                y = len(S)-1
                while y>=0 and m[S[y]][x]!=S[y]:
                    y = y-1
                if y<0: S.append(x)
        uc.append(S)
    return uc

class Poset(RelStructure):
    def __init__(self,uppercovers):
        self.uc = dict([(x,uppercovers[x]) for x in range(len(uppercovers))])

    def size(self):
        return len(self.uc)

    def get_leq(self):
        if hasattr(self,'leq'): return self.leq
        n = self.size()
        leq = [[0 for y in range(n)] for x in range(n)]
        for i in range(n):
            leq[i][i] = 1
            for j in self.uc[i]:
                leq[i][j] = 1
        for i in range(n):
            for j in range(i+1,n):
                if leq[i][j]:
                    for k in range(j+1,n):
                        if leq[j][k]: leq[i][k] = 1
        self.leq = leq
        return leq

    def get_downsets(self):
        if hasattr(self,'down'): return self.down
        n = self.size()
        le = self.get_leq()
        down = [[x for x in range(n) if le[x][y]==1] for y in range(n)]
        self.down = down
        return down

    def get_lowercovers(self):
        if hasattr(self,'lc'): return self.lc
        lc = [[] for x in range(self.size())]
        for i in range(self.size()):
            for j in self.uc[i]:
                lc[j].append(i)
        self.lc = lc
        return lc

    def get_join(self): # Freese-Jezek-Nation p217
        if hasattr(self,'join'): return self.join
        n = self.size()
        join = [[0 for x in range(n)] for x in range(n)]
        le = self.get_leq()
        if not all([le[x][n-1]==1 for x in range(n)]):
            return "poset has no top element"
        p = range(n-1,-1,-1)
        uc = [sorted([p[y] for y in self.uc[x]]) for x in p]
        S = []
        for x in range(n): # x=x_k
            join[x][x] = x
            for y in S:
                T = []
                for z in uc[x]:
                    T.append(join[y][z]) # T = {x_i \vee z : z>-x_k}
                q = T[0]
                for z in T[1:]:
                    if z>q: q = z #error in Listing 11.9
                for z in T:
                    if not le[p[q]][p[z]]:
                        return "poset is not a semilattice: x="+x+" y="+y
                join[x][y] = q
                join[y][x] = q
            S.append(x)
        self.join = permuted_binary_op(join,p)
        return self.join

    def get_meet(self): # Freese-Jezek-Nation p217
        if hasattr(self,'meet'): return self.meet
        n = self.size()
        meet = [[0 for x in range(n)] for x in range(n)]
        le = self.get_leq()
        if not all([le[0][x]==1 for x in range(n)]):
            return "poset has no bottom element"
        lc = self.get_lowercovers()
        S = []
        for x in range(n): # x=x_k
            meet[x][x] = x
            for y in S:
                T = []
                for z in lc[x]:
                    T.append(meet[y][z]) # T = {x_i \wedge z : z>-x_k}
                q = T[0]
                for z in T[1:]:
                    if z>q: q = z
                for z in T:
                    if not le[z][q]:
                        return "poset is not a semilattice: x="+x+" y="+y
                meet[x][y] = q
                meet[y][x] = q
            S.append(x)
        self.meet = meet
        return meet

    def is_distributive(self):
        jn = self.get_join()
        mt = self.get_meet()
        n = self.size()
        for x in range(n):
            for y in range(n):
                for z in range(n):
                    if mt[x][jn[y][z]]!=jn[mt[x][y]][mt[x][z]]: return False
        return True

    def is_selfdual(self):
	lc = self.get_lowercovers()
        n = self.size()
        perms = [[n-1]+p+[0] for p in permutations(1,n-1)]
        for p in perms:
            plc = {}
            for i in range(n):
                plc[i] = [p[y] for y in lc[p[i]]]
                plc[i].sort()
            if plc == self.uc: return True
	return False

    def get_uniqueLC(self): # find subset of elements with unique lower cover
        lc = self.get_lowercovers()
        return [x for x in range(self.size()) if len(lc[x])==1]

    def get_uniqueUC(self):
        return [x for x in range(self.size()) if len(self.uc[x])==1]

    def get_basicsets(self):
        n = self.size()
        le = self.get_leq()
        ji = self.get_uniqueLC()
        mi = self.get_uniqueUC()
        return [[x for x in ji if le[x][y]==1] for y in mi]

    def subposet(self,S): # S a subset of the elements of self
        leq = self.get_leq()
        leq = [[leq[x][y] for y in S] for x in S]
        Q = Poset(leq2uc(leq))
        if hasattr(self,'labels'): Q.label = [self.labels[x] for x in S]
        return Q

    def get_labels(self):
        if hasattr(self,'labels'): return self.labels
        return range(self.size())

    def heights(self): # assumes P is topologically sorted
        l=[0 for x in range(self.size())]
        for i in range(self.size()):
            for j in self.uc[i]:
                if l[j]<l[i]+1: l[j]=l[i]+1
        return l

    def depths(self): # assumes P is topologically sorted
        l=[0 for x in range(self.size())]
        lc=self.get_lowercovers();
        for i in range(self.size()-1,-1,-1):
            for j in lc[i]:
                if l[j]<l[i]+1: l[j]=l[i]+1
        return l

    def levels(self): # P topological
        d=self.depths()
        m=max(d)
        l=[[] for x in range(m+1)]
        for x in range(self.size()):
            l[m-d[x]].append(x)
        return l

    def get_pos(self): # get [x,y] position for each element
        if hasattr(self,'pos'): return self.pos
        lev = self.levels();
        lengths = [len(x) for x in lev]
        m = max(lengths)
        irreg = any([x%2==0 for x in lengths]) and \
                any([x%2==1 for x in lengths])
        if irreg and m%2==0: m = m+1
        pos = [[] for i in range(self.size())]
        for l in range(len(lev)):
            n = len(lev[l])
            for i in range(n):
                if irreg and n%2==0:
                    if 2*i<=n:
                        pos[lev[l][i]] = [(m-n-1)/2+i,l]
                    else:
                        pos[lev[l][i]] = [(m-n-1)/2+i+1,l]
                else:
                    pos[lev[l][i]] = [(m-n)/2+i,l]
        self.pos = pos
        return pos

    def html(self):
        emb = self.get_pos()
        nam = self.get_labels()
        if type(nam[0]) is list and type(nam[0][0]) is str:
            nam = [join(nam[x],';') for x in range(self.size())]
        st="""<html><body>\n<applet codebase=\"http://www.chapman.edu/~jipsen/gap\" code=\"Poset.class\" width=\"800\" height=\""""
        maxY = max([x[1] for x in emb])
        st += str(30*maxY+2*20+38)
        st += """\">\n<param name=fix value=\"true\">\n<param name=repel value=\"10\">\n<param name=attract value=\"0\">\n<param name=maxX value=\""""
        st += str(max(25,max([x[0] for x in emb]))) #width
        st += """\">\n<param name=maxY value=\""""
        st += str(maxY) #height
        st += """\">\n<param name=poset value='\n"""
        uc = self.uc
        n = len(uc)
        for i in range(n): # for each elt at this level
            st += """["""+str(i)+""", """+str(emb[i][0])+""","""+\
                  str(emb[i][1])+""", \""""+str(nam[i])+"""\", """+\
                  str(uc[i]).replace(' ','')+"""],\n"""
        st = st[:-2]
        st += """'>\n</applet>\n</body></html>\n"""
        print st

def Chain(n):
    c = [[x+1] for x in range(n)]
    c[n-1] = []
    return Poset(c)

def Antichain(n):
    c = [[] for x in range(n)]
    return Poset(c)

def Pentagon():
    p = Poset([[1,2],[4],[3],[4],[]])
    p.pos = [[2,0],[0,2],[3,1],[3,3],[2,4]]
    return p

def Diamond(n):
    c = [[n-1] for x in range(n)]
    c[0] = [x for x in range(1,n-1)]
    c[n-1] = []
    return Poset(c)

def opstr(m):  # convert 2-dim list to a compact string for display
        nr = len(m)
        if nr == 0: return "[]"
        nc = len(m[0])
        s = [[str(y).replace(' ','') for y in x] for x in m]
        width = [max([len(s[x][y]) for x in range(nr)]) for y in range(nc)]
        s = [[" "*(width[y]-len(s[x][y]))+s[x][y] for y in range(nc)]\
             for x in range(nr)]
        rows = ["["+",".join(x)+"]" for x in s]
        s = "[\n"+",\n".join(rows)+"]"
        return s

def is_automorphism(p,m): #p[m[i][j]] == m[p[i]][p[j]]
    n = len(m)
    for i in range(n):
        for j in range(n):
            if p[m[i][j]] != m[p[i]][p[j]]: return False
    return True

def flat(m):
    return reduce(lambda x,y:x+y,m)

def printto(fname,obj):
    fh = open(fname,'w')
    fh.write(str(obj))
    fh.close()

def savelist(fname,lst,doc=""):
    fh = open(fname,'w')
    if doc!="": fh.write('# '+doc)
    fh.write('[\n%s]'%(',\n'.join([repr(a) for a in lst])))
    fh.close()

def findandsave(clsname,clsabbr,clsfind,n):
    algs = [clsfind(i) for i in range(1,n+1)]
    count = [len(x) for x in algs]
    algs = reduce(lambda x,y:x+y,algs)
    for i in range(len(algs)): algs[i].num = i
    fh = open(clsabbr+str(n)+'.txt','w')
    fh.write('# %s up to size %s\n\n%s_count = %s # total = %s\n%s_index = %s\
    \n\n%s_list = [\n%s]\n\n# Computed by Peter Jipsen with %s on %s\
    \n# This file can be used in Python with execfile("filename")'\
             %(clsname,n,clsabbr,count,len(algs),clsabbr,\
               [sum(count[:i]) for i in range(n)],clsabbr,\
               ',\n'.join([repr(a) for a in algs]),version()[0:version()\
               .index(',')],datetime.datetime.today().date()))
    fh.close()
    return algs

class POGroupoid(Poset):
    def __init__(self,uc,mult):
        self.uc = dict([(x,uc[x]) for x in range(len(uc))])
        self.mult = mult

    def get_backslash(self): #calculate the left residual (assumed to exist)
        if hasattr(self,'bs'): return self.bs  # xz <= y iff z <= x\y
        mult = self.mult
        join = self.get_join()
        n = self.size()
        bs = [[0 for x in range(n)] for y in range(n)]
        for x in range(n):
            for y in range(n):
                z = n-1
#                print (join,mult)
                while z>=0 and join[mult[x][z]][y]!=y: z = z-1
                if z<0: raise "x\y not defined for x="+str(x)+" y="+str(y)
                bs[x][y] = z
        self.bs = bs
        return bs

    def get_slash(self): #calculate the right residual (assumed to exist)
        if hasattr(self,'sl'): return self.sl  # zy <= x iff z <= x/y
        mult = self.mult
        join = self.get_join()
        n = self.size()
        sl = [[0 for x in range(n)] for y in range(n)]
        for x in range(n):
            for y in range(n):
                z = n-1
                while z>=0 and join[mult[z][y]][x]!=x: z = z-1
                if z<0: raise "x/y not defined for x="+x+" y="+y
                sl[x][y] = z
        self.sl = sl
        return sl

    def is_lattice_ordered(self):
        jn = self.get_join()
        m = self.mult
        n = self.size()
        for x in range(n):
            for y in range(n):
                for z in range(n):
                    if m[x][jn[y][z]]!=jn[m[x][y]][m[x][z]]: return False
                    if m[jn[x][y]][z]!=jn[m[x][z]][m[y][z]]: return False
        return True

    def is_GBL(self):
        bs = self.get_backslash()
        sl = self.get_slash()
        n = self.size()
        for x in range(n):
            for y in range(n):
                if self.join[x][y]==y and \
                   not (x==self.mult[y][bs[y][x]] and \
                        x==self.mult[sl[x][y]][y]): return False
        return True

    def is_GMV(self):
        bs = self.get_backslash()
        sl = self.get_slash()
        n = self.size()
        for x in range(n):
            for y in range(n):
                if self.join[x][y]==y and \
                   not (y==sl[x][bs[y][x]] and \
                        y==bs[sl[x][y]][x]): return False
        return True

    def is_MV(self):
        return self.is_commutative() and self.is_GMV()

    def is_involutive(self,d=None):
        bs = self.get_backslash()
        sl = self.get_slash()
        if d == None: d = self.d
        for x in range(self.size()):
            if sl[d][bs[x][d]]!=x: return False
            if bs[sl[d][x]][d]!=x: return False
        return True

    def is_cyclic(self,d=None):
        bs = self.get_backslash()
        sl = self.get_slash()
        if d == None: d = self.d
        for x in range(self.size()):
            if bs[x][d]!=sl[d][x]: return False
        return True

    def get_automorphisms(self): # assumes 0=bottom, n-1=top
        n = self.size()
        join = self.get_join()
        perms = [[0]+p+[n-1] for p in permutations(1,n-1)]
        return [p for p in perms if is_automorphism(p,join) and\
                is_automorphism(p,self.mult)]

    def all_involutive(self):
        # return all involutive FL-algs with the RL reduct "self"
        n = self.size()
        infl_list = []
        aut = self.get_automorphisms()
        for d in range(n):
            if self.is_involutive(d) and all([p[d]>=d for p in aut]):
                infl_list.append(InFL(self.uc,self.mult,d))
        return infl_list

    def is_reduct_of_involutive(self):
        return self.all_involutive()!=[]

class RL(POGroupoid):
    
    def __init__(self,uc,mult,e=None,num=None):
        self.mult = mult
        self.uc = dict([(x,uc[x]) for x in range(len(uc))])
        self.e = mult.index(range(len(mult)))
        if num!=None: self.num = num

    def _repr_(self):
        return "\nRL(e = %s, "%self.e+\
               ("num = %s,"%self.num if hasattr(self,'num') else "")+\
               "\nuc = %s,\nmult = %s)"%(self.uc,opstr(self.mult))

    @staticmethod
    def find(n):
        return \
           findandsave('Residuated lattices','RL',find_residuated_lattices,n)

    def get_rl_frame(self):
        le = self.get_leq()
        W = self.get_uniqueLC()
        E = [z for z in W if le[z][self.e]]
        bclo = self.get_basicsets()
        circ = [[[z for z in W if le[z][self.mult[x][y]]] for y in W]\
                for x in W]
        return RLFrame(W,bclo,circ,E,self.num)

class RLFrame(RelStructure):
    def __init__(self,W,bclo,circ,E,num=None):
        self.W = W
        self.E = E
        self.bclo = bclo # basic closed sets
        self.circ = dict([(W[x],dict([(W[y],circ[x][y])\
                          for y in range(len(W))])) for x in range(len(W))])
        if num!=None: self.num = num

    def _repr_(self):
        return "\nRLFrame(W = %s, E = %s, "%(self.W,self.E)+\
               ("num = %s,"%self.num if hasattr(self,'num') else "")+\
               "\nbclo = %s,\ncirc = %s)"%(self.bclo,\
               opstr([[self.circ[x][y] for y in self.W] for x in self.W]))

class InFL(POGroupoid):
    
    def __init__(self,uc,mult,d,e=None,num=None):
        self.mult = mult
        self.uc = dict([(x,uc[x]) for x in range(len(uc))])
        self.e = mult.index(range(len(mult)))
        self.d = d
        if num!=None: self.num = num

    def _repr_(self):
        bs = self.get_backslash()
        sl = self.get_slash()
        tilde = dict([(x,bs[x][self.d]) for x in range(self.size())])
        minus = dict([(x,sl[self.d][x]) for x in range(self.size())])
        return "\nInFL(e = %s, d = %s, "%(self.e,self.d)+\
               ("num = %s,"%self.num if hasattr(self,'num') else "")+\
               "\nuc = %s,\n# tilde = %s%s\nmult = %s)"\
               %(self.uc,tilde,"\n# minus = "+str(minus) \
                 if minus!=tilde else " = minus",opstr(self.mult))

    @staticmethod
    def find(n):
        return \
           findandsave('Involutive FL-algebras','InFL',find_involutive_FL,n)

    def get_dualmult(self):
        n = self.size()
        bs = self.get_backslash()
        sl = self.get_slash()
        return [[bs[self.mult[sl[self.d][y]][sl[self.d][x]]][self.d]\
                 for y in range(n)] for x in range(n)]

    def get_infl_frame(self):
        le = self.get_leq()
        W = self.get_uniqueLC()
        E = [z for z in W if le[z][self.e]]
        D = [z for z in W if le[z][self.d]]
        bclo = self.get_basicsets()
        circ = [[[z for z in W if le[z][self.mult[x][y]]] for y in W]\
                for x in W]
        bs = self.get_backslash()
        sl = self.get_slash()
        return InFLFrame(W,bclo,circ,E,D,(self.num if hasattr(self,'num') \
                         else None),\
            dict([(x,[z for z in W if le[z][bs[x][self.d]]]) for x in W]),\
            dict([(x,[z for z in W if le[z][sl[self.d][x]]]) for x in W]))

class InFLFrame(RLFrame):
    def __init__(self,W,bclo,circ,E,D,num=None,tilde=None,minus=None):
        RLFrame.__init__(self,W,bclo,circ,E,num)
        self.D = D
        self.tilde = tilde
        self.minus = minus

    def _repr_(self):
        return "\nInFLFrame(W = %s, E = %s, D = %s, "%(self.W,self.E,self.D)+\
               ("num = %s,"%self.num if hasattr(self,'num') else "")+\
               "\nbclo = %s,\n# tilde = %s%s\ncirc = %s)"%(self.bclo,\
               self.tilde,'n# minus = '+str(self.minus) \
               if self.tilde != self.minus else ' = minus',\
               opstr([[self.circ[x][y] for y in self.W] for x in self.W]))

    def get_external_negations(self):
        """
        find all pairs t,m such that (circ,t,m) is a bi-involutive monoid
        """
        W = self.W
        n = len(W)
        perms = permutations(0,n)
        biinv = []
        for p in perms:
            q = inverse_permutation(p)
            pw = dict([(W[i],W[p[i]]) for i in range(n)])
            if p <= q:
                if all([[pw[z] for z in self.circ[x][y]]==\
                        self.circ[pw[y]][pw[x]] for x in W for y in W]):
                    biinv.append([p,q])
        return biinv

def is_partially_ordered(m,uc,le):
    n = len(m)
    for x in range(n):
        for y in uc[x]:
            for z in range(n):
                if m[x][z]>-1 and m[y][z]>-1 and le[m[x][z]][m[y][z]]==0:
                    return False
                if m[z][x]>-1 and m[z][y]>-1 and le[m[z][x]][m[z][y]]==0:
                    return False
    return True

def is_op_automorphism(m,q):
    n = len(m)
    for i in range(n):
        for j in range(n):
            if m[q[i]][q[j]]!=q[m[i][j]]: return false
    return true

def is_rel_automorphism(m,q):
    n = len(m)
    for i in range(n):
        for j in range(n):
            if m[q[i]][q[j]]!=m[i][j]: return false
    return true

def is_po_automorphism(uc,q):
    n = len(uc)
    for i in range(n):
        for j in uc[i]:
            if not q[j] in uc[q[i]]: return false
    return true

def complete_po_operation(m,uc,le,i,j,associative=False,commutative=False):
    """
    find next i,j where alg[i][j] = -1 = undefined; for each val=0..n-1
    set alg[i][j] = val, check axioms and permutations
    if ok, complete_ordered_operation(m,i,j+1) then restore and return
    """
    global algebra_list,entries_list,pflags,pflags_list
    ok = True
    n = len(m)
    while ok and i<n:
        while j<n and m[i][j]>-1: j = j+1
        if j>=n:
            j = 0 
            i = i+1
        else: ok = False
    if ok: algebra_list.append([m[i][:] for i in range(n)])
    else:
        for val in range(n):
            ok = True
            el = len(entries_list)
            m[i][j] = val
            entries_list += [(i,j)]
            if commutative:
                m[j][i] = val
                entries_list += [(j,i)]
            if associative: ok = check_associativity(m)
            if ok: ok = is_partially_ordered(m,uc,le)
            pl = len(pflags_list)
            if ok: ok = check_permutations(m)
            if ok: complete_po_operation(m,uc,le,i,j+1,associative,commutative)
            while len(entries_list)>el:
                e = entries_list.pop()
                m[e[0]][e[1]] = -1
            while len(pflags_list)>pl:
                p = pflags_list.pop()
                pflags[p] = True

def find_po_groupoids(po,associative=False,commutative=False,idempotent=False,identity=False,id=None,zero=False):
    """
    Return list of partially ordered groupoids of size  n  as nxn matrices.
    The ordering is given by the list uc of upper covers: y in uc[x] iff x-<y.
    The operation is compatible with the order.
    """
    global algebra_list,perms,invperms,entries_list,pflags,pflags_list
    algebra_list=[]
    entries_list=[]
    n = po.size()
    le = po.get_leq()
    m = [[-1 for x in range(n)] for x in range(n)]
    perms = [p for p in permutations(0,n) if is_po_automorphism(po.uc,p)]
    if idempotent:
        for x in range(n):
            m[x][x] = x
    if zero:
        perms = [[0]+p for p in permutations(1,n)]
        for x in range(n):
            m[x][0] = 0
            m[0][x] = 0
    if identity:
        if id==None: #should only try all canconical identity positions
            t = n
            if zero: s = 1
            else: s = 0
        else: s = id; t = id+1
        for e in range(s,t):
            for x in range(n):
                m[e][x] = x
                m[x][e] = x
            perms = [p for p in permutations(0,n) if p[e]==e] #keep e fixed
            perms = [p for p in perms if is_po_automorphism(po.uc,p)]
            invperms = [inverse_permutation(p) for p in perms]
            pflags = [True for i in range(len(perms))]
            pflags_list = []
            complete_po_operation(m,po.uc,le,0,0,associative,commutative)
            for x in range(n):
                m[e][x] = -1
                m[x][e] = -1
    else:
        invperms = [inverse_permutation(p) for p in perms]
        pflags = [True for i in range(len(perms))]
        pflags_list = []
        complete_po_operation(m,po.uc,le,0,0,associative,commutative)
    return algebra_list

def find_po_semigroups(po):
    return find_po_groupoids(po,associative=True)

def find_commutative_po_semigroups(po):
    return find_po_groupoids(po,associative=True,commutative=True)

def find_po_bands(po):
    return find_po_groupoids(po,associative=True,idempotent=True)

def find_po_semilattices(po):
    return find_po_groupoids(po,associative=True,commutative=True,idempotent=True)

def find_po_monoids(po,e=None):
    return find_po_groupoids(po,associative=True,identity=True,id=e)

def find_commutative_po_monoids(po,e=None):
    return find_po_groupoids(po,associative=True,commutative=True,identity=True,id=e)

def find_po_monoids_with_zero(po,e=None):
    return find_po_groupoids(po,associative=True,identity=True,id=e,zero=True)

def find_commutative_po_monoids_with_zero(po,e=None):
    return find_po_groupoids(po,associative=True,commutative=True,identity=True,id=e,zero=True)

def find_integral_po_monoids_with_zero(po):
    return find_po_monoids_with_zero(po,po.size()-1)

def find_commutative_integral_po_monoids_with_zero(po):
    return find_po_commutative_monoids_with_zero(po,po.size()-1)

###############################################################################
# Search for lattice ordered structures (join-preserving)

def check_lattice_ordered(m,join):
    """
    Return False if partial optable  m  does not distribute over join
    Return True otherwise
    """
    global entries_list
    n = len(m)
    for x in range(n):
        for y in range(n):
            xy = m[x][y]
            if xy>-1:
                for z in range(n):
                    xz = m[x][z]
                    if xz>-1:
                        yz = join[y][z]
                        xyz = m[x][yz]
                        if xyz>-1:
                            if xyz!=join[xy][xz]: return False
                        else:
                            m[x][yz] = join[xy][xz]
                            entries_list += [(x,yz)]
                    zy = m[z][y]
                    if zy>-1:
                        xz = join[x][z]
                        xyz = m[xz][y]
                        if xyz>-1:
                            if xyz!=join[xy][zy]: return False
                        else:
                            m[xz][y] = join[xy][zy]
                            entries_list += [(xz,y)]
    return True

def check_divisible(m,down):
    """
    Return False if partial optable  m  is not divisible
    i.e. x<=y implies x=m[y][u]=m[v][y] for some u,v
    Return True otherwise
    Compute list S of elements in y-row (col).
    Compute k = |down[y] - S|.
    if k>#undefined in y-row, return False
    if #undefined=1 and k=1 let undefined = this elt.
    """
    global entries_list
    n = len(m)
    for y in range(n):
        S = set([m[y][u] for u in range(n) if m[y][u]>-1])
        E = set(down[y]).difference(S)
        U = [u for u in range(n) if m[y][u]==-1]
        if len(E)>len(U): return False
        if len(U)==1 and len(E)==1:
            m[y][U[0]] = E.pop()
            entries_list += [(y,U[0])]
    for y in range(n):
        S = set([m[u][y] for u in range(n) if m[u][y]>-1])
        E = set(down[y]).difference(S)
        U = [u for u in range(n) if m[u][y]==-1]
        if len(E)>len(U): return False
        if len(U)==1 and len(E)==1:
            m[U[0]][y] = E.pop()
            entries_list += [(U[0],y)]
    return True

def complete_l_operation(m,join,down,uc,i,j,associative=False,commutative=False,divisible=False):
    """
    find next i,j where alg[i][j] = -1 = undefined; for each val=0..n-1
    set alg[i][j] = val, check axioms and permutations
    if ok, complete_ordered_operation(m,i,j+1) then restore and return
    """
    global algebra_list,entries_list,pflags,pflags_list
    ok = True
    n = len(m)
    while ok and i<n:
        while j<n and m[i][j]>-1: j = j+1
        if j>=n:
            j = 0 
            i = i+1
        else: ok = False
    if ok:
        if check_lattice_ordered(m,join) and \
               (not divisible or check_divisible(m,down)):
            algebra_list.append(POGroupoid(uc,[x[:] for x in m]))
    else:
        for val in range(n):
            if all([join[x][i]!=i or m[x][j]==-1 or join[m[x][j]][val]==val \
                    for x in range(i)]) and \
               all([join[x][j]!=j or m[i][x]==-1 or join[m[i][x]][val]==val \
                    for x in range(j)]):
                ok = True
                el = len(entries_list)
                m[i][j] = val
                entries_list += [(i,j)]
                if commutative:
                    m[j][i] = val
                    entries_list += [(j,i)]
                if ok: ok = check_lattice_ordered(m,join)
                if ok and divisible: ok = check_divisible(m,down)
                if ok and associative: ok = check_associativity(m)
                pl = len(pflags_list)
                if ok: ok = check_permutations(m)
                if ok: complete_l_operation(m,join,down,uc,i,j+1,associative,commutative,divisible)
                while len(entries_list)>el:
                    e = entries_list.pop()
                    m[e[0]][e[1]] = -1
                while len(pflags_list)>pl:
                    p = pflags_list.pop()
                    pflags[p] = True

def find_l_groupoids(l,associative=False,commutative=False,divisible=False,idempotent=False,identity=False,id=None,zero=False):
    """
    Return list of partially ordered groupoids of size  n  as nxn matrices.
    The ordering is given by the list uc of upper covers: y in uc[x] iff x-<y.
    The operation is compatible with the order.
    """
    global algebra_list,perms,invperms,entries_list,pflags,pflags_list
    algebra_list=[]
    entries_list=[]
    n = l.size()
    join = l.get_join()
    down = l.get_downsets()
    m = [[-1 for x in range(n)] for x in range(n)]
    perms = [[0]+p+[n-1] for p in permutations(1,n-1) \
             if is_po_automorphism(l.uc,[0]+p+[n-1])]
    if idempotent:
        for x in range(n):
            m[x][x] = x
    if zero:
        perms = [[0]+p for p in permutations(1,n)]
        for x in range(n):
            m[x][0] = 0
            m[0][x] = 0
    if identity:
        if id==None: #should only try all canconical identity positions
            t = n
            if zero and n>1: s = 1
            else: s = 0
        else: s = id; t = id+1
        for e in range(s,t):
            for x in range(n):
                m[e][x] = x
                m[x][e] = x
            perms = [[0]+p+[n-1] for p in permutations(1,n-1) \
                     if is_po_automorphism(l.uc,[0]+p+[n-1])]
            if all([p[e]>=e for p in perms]): # eliminate isomorphic (P,e)
                perms = [p for p in perms if p[e]==e] #keep e fixed
                invperms = [inverse_permutation(p) for p in perms]
                pflags = [True for i in range(len(perms))]
                pflags_list = []
                complete_l_operation(m,join,down,l.uc,0,0,associative,commutative,divisible)
            for x in range(n):
                m[e][x] = -1
                m[x][e] = -1
    else:
        invperms = [inverse_permutation(p) for p in perms]
        pflags = [True for i in range(len(perms))]
        pflags_list = []
        complete_l_operation(m,join,down,l.uc,0,0,associative,commutative,divisible)
    return algebra_list

def find_l_semigroups(l):
    return find_l_groupoids(l,associative=True)

def find_l_commutative_semigroups(l):
    return find_l_groupoids(l,associative=True,commutative=True)

def find_l_bands(l):
    return find_l_groupoids(l,associative=True,idempotent=True)

def find_l_semilattices(l):
    return find_l_groupoids(l,associative=True,commutative=True,idempotent=True)

def find_l_monoids(l,e=None):
    return find_l_groupoids(l,associative=True,identity=True,id=e)

def find_l_commutative_monoids(l,e=None):
    return find_l_groupoids(l,associative=True,commutative=True,identity=True,id=e)

def find_l_monoids_with_zero(l,e=None):
    return find_l_groupoids(l,associative=True,identity=True,id=e,zero=True)

def find_commutative_l_monoids_with_zero(l,e=None):
    return find_l_groupoids(l,associative=True,commutative=True,identity=True,id=e,zero=True)

def find_integral_l_monoids_with_zero(l):
    return find_l_monoids_with_zero(l,l.size()-1)

def find_commutative_integral_l_monoids_with_zero(l):
    return find_l_commutative_monoids_with_zero(l,l.size()-1)

def find_residuated_lattices(n):
    algs = []
    for l in find_lattices(n):
        algs.extend(find_l_monoids_with_zero(l))
    return [RL(l.uc,l.mult) for l in algs]

def find_selfdual_lattices_slow(n):
    return [x for x in find_lattices(n) if x.is_selfdual()]

def find_selfdual_residuated_lattices(n):
    alg_list = []
    for l in find_selfdual_lattices_slow(n):
        alg_list.extend(find_l_monoids_with_zero(l))
    return alg_list

def find_involutive_FL(n):
    alg_list = []
    for l in find_selfdual_lattices_slow(n):
        alg_list.extend(find_l_monoids_with_zero(l))
    fl = []
    for l in alg_list:
        fl.extend(l.all_involutive())
    return fl

def find_GBL(l): # of course they are all commutative [Jipsen-Montagna 2006]
    return find_l_groupoids(l,associative=True,divisible=True,identity=True,id=l.size()-1,zero=True)

def find_basic_logic_algebras(l):
    return find_l_groupoids(l,associative=True,commutative=True,divisible=True,identity=True,id=l.size()-1,zero=True)

###############################################################################
# Search for binary relations

def check_relation_permutations(m):
    """
    Return False if every completion of some isomorphic copy q(m) of m
    will be lexicographically smaller than m
    Return True otherwise
    """
    global perms, invperms, pflags, pflags_list
    n = len(m)
    for k in range(len(perms)):
        if pflags[k]:
            q = perms[k]
            qi = invperms[k]
            equal = True
            for i in range(n):
                if equal: #equal means m[qi[i]][qi[j]] == m[i][j] != -1
                    for j in range(n):
                        if equal:
                            mqi = m[qi[i]][qi[j]]
                            mij = m[i][j]
                            equal = (mqi>-1 and mij==mqi)
                            if mqi>-1 and mqi<mij: return False
                            if mqi>-1 and mij>-1 and mqi>mij: #get rid of q
                                pflags[k] = False
                                pflags_list += [k]
    return True

def complete_relation(m,i,j,symmetric=False):
    """
    find next i,j where alg[i][j] = -1 = undefined; for each val=0..n-1
    set alg[i][j] = val, check permutations
    if ok, complete_operation(m,i,j+1) then restore and return
    """
    global pflags_list
    ok = True
    n = len(m)
    while ok and i<n:
        while j<n and m[i][j]>-1: j = j+1
        if j>=n:
            j = 0 
            i = i+1
        else: ok = False
    if ok: relation_list.append([m[i][:] for i in range(n)])
    else:
        for val in range(2):
            m[i][j] = val
            if symmetric: m[j][i] = val
            pl = len(pflags_list)
            ok = check_relation_permutations(m)
            if ok: complete_relation(m,i,j+1,symmetric)
            m[i][j] = -1
            if symmetric: m[j][i] = -1
            while len(pflags_list)>pl:
                p = pflags_list.pop()
                pflags[p] = True

def find_binary_relations(n,irreflexive=False,symmetric=False):
    """
    Return list of nonisomorphic binary relations of size  n  as 0-1-matrices
    """
    global relation_list,perms,invperms,entries_list,pflags,pflags_list
    relation_list=[]
    perms = [p for p in permutations(0,n)]
    invperms = [inverse_permutation(p) for p in perms]
    pflags = [True for i in range(len(perms))]
    pflags_list = []
    m = [[-1 for x in range(n)] for x in range(n)]
    if irreflexive:
        for x in range(n):
            m[x][x]=0
    complete_relation(m,0,0,symmetric)
    return relation_list

def find_simple_digraphs(n):
    """
    Return list of nonisomorphic simple digraphs of size  n  as adjacency matrices
    Simple digraphs are binary relations without loops
    """
    return find_binary_relations(n,irreflexive=True)

def find_symmetric_relations(n):
    """
    Return list of nonisomorphic symmetric binary relations of size  n  as adjacency matrices
    Symmetric binary relations correspond to graphs with loops allowed
    """
    return find_binary_relations(n,symmetric=True)

def find_simple_graphs(n):
    """
    Return list of nonisomorphic graphs of size  n  as adjacency matrices
    Here graphs are symmetric binary relations without loops
    """
    return find_binary_relations(n,irreflexive=True,symmetric=True)

def find_relations_test(check_sloane=False):
    tl = [
        [len(find_binary_relations(i)) for i in range(1,5)],
        [len(find_simple_digraphs(i)) for i in range(1,6)],      #binary relations without loops
        [len(find_symmetric_relations(i)) for i in range(1,6)],
        [len(find_simple_graphs(i)) for i in range(1,7)]]        #simple digraphs with arrows in both directions
    if not check_sloane: return tl
    sl = [[t]+sloane_find(t,1)[0] for t in tl]
    return [[s[0],s[3][:8],s[2]] for s in sl]