CoCalc Public Files2021-01-17-124214.sagews
Author: Stephen Kauffman
Views : 93
Description: Algebra to solve logic puzzles Convert boolean polynomials to probability
Compute Environment: Ubuntu 20.04 (Default)
load('ProbRecurse.py')
import numpy
import scipy.optimize

load('ProbRecurse.py')
import numpy
import scipy.optimize
#from operator import mul

KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])
GIVEN = lambda A, B: ProbM(AND(A,B))/ProbM(B)
convert = lambda f: reduce(lambda a,b:a+b,[reduce(lambda a,b:a*b,[(m.degree(g)+1+g) for g in BPR.gens()]) for m in f.monomials()]) #self inverse convert poly in gens to poly in atoms (missing var x implies, interpret as ~x)
# convert monomial terms A^i*B^j*C^k --> (A+1+i)*(B+1+j)*(C*1*k) in Boolean Ring i,j,k=0,1

objects = 6 #rel_classes*types; print('objects = ' + str(objects))
BPR = BooleanPolynomialRing(objects^2,'X',order='degneglex')
BPR.inject_variables()
print(str(objects^2) + ' Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X' + str(objects^2-1) + ' = R[' + str(objects-1) + ',' + str(objects-1) + ']')
X = BPR.gens()
NewGens = str(X).lower()
NewGens = NewGens[1:len(NewGens) - 1]
PPR = PolynomialRing(QQ,objects^2,NewGens,order='degneglex')
PPR.inject_variables()
R = matrix(BPR, objects, objects, BPR.gens())
prems = [];reflex=[];symmtrc=[];trans2=[];direct=[]
for i in range(0, objects):    #reflexive
reflex.append(R[i,i])

for i in range(0,objects):
for j in range(0, objects):
#symmtrc.append(IF(R[i,j],R[j,i]))
direct.append(IF(R[i,j],NOT(R[j,i])))

for i in range(0,objects):    #symmetric
for j in range(i+1, objects):
symmtrc.append(IFF(R[i,j],R[j,i]))
symmtrc=list(set(symmtrc))

for j in range(0, objects):    #transitive
for i in range(0, objects):
for k in range(0, objects):
poly = IF(AND(R[min(i,j),max(i,j)],R[min(j,k),max(j,k)]),R[min(i,k),max(i,k)])
trans2.append(poly)
trans2 = list(set(trans2))

M=matrix([[125,25,5,1],[64,16,4,1],[27,9,3,1],[8,4,2,1]])#;M
var('t')
B=M.inverse()*vector([35,18,8,3])#;B
mypoly2=vector([t^3,t^2,t,1])*B;mypoly2
# mypoly2.subs(t=6);mypoly2.subs(t=7) #seems to give the number of polynomials in groebner basis for equiv relation with gen order degneglexx
prems=reflex+symmtrc+trans2
if BPR(0) in prems:
prems = list(set(prems)) # remove dup
premstrue = [p+1 for p in prems] #make true
REI = ideal(premstrue);# print('REI.gens() = ' + str(len(REI.gens())))
GB = REI.groebner_basis();GB
#compute groebner basis directly assuming equivalence relation
#correspond to trans2 objects=5: X1*(X2+X(2+1*5));X1*(X3+X(3+1*5));X1*(X4+X(4+1*5));X2*(X3+X(3+2*5));X2*(X4+X(4+2*5));X3*(X4+X(4+3*5));X7*(X8+X(8+1*5));X7*(X9+X(9+1*5));X8*(X9+X(9+2*5));X13*(X14+X(14+1*5)) then
#X2*(X1+X(2+1*5));X3*(X1+X(3+1*5));X4*(X1+X(4+1*5));X3*(X2+X(3+2*5));X4*(X2+X(4+2*5));X4*(X3+X(4+3*5));X8*(X7+X(8+1*5));X9*(X7+X(9+1*5));X9*(X8+X(9+2*5));X14*(X13+X(14+1*5)) that 20 polynomials
GB2 = []
PGB = []
for i in range(objects):
for j in range(i+1,objects):
for k in range(j+1,objects):
print(i);print(j);print(k)
print(R[i,j]*(R[i,k]+R[j,k]))
GB2.append(R[i,j]*(R[i,k]+R[j,k]))
print(PPR(ProbM(R[i,k])))
PGB.append(PPR(ProbM(R[i,k]))-PPR(ProbM(R[j,k])))

Defining X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20, X21, X22, X23, X24, X25, X26, X27, X28, X29, X30, X31, X32, X33, X34, X35 36 Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X35 = R[5,5] Defining x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31, x32, x33, x34, x35 t 1/3*t^3 - 1/2*t^2 + 7/6*t Polynomial Sequence with 61 Polynomials in 36 Variables 0 1 2 X1*X8 + X1*X2 Defining x2 x2 Defining x2 Defining x8 0 1 3 X1*X9 + X1*X3 Defining x3 x3 Defining x3 Defining x9 0 1 4 X1*X10 + X1*X4 Defining x4 x4 Defining x4 Defining x10 0 1 5 X1*X11 + X1*X5 Defining x5 x5 Defining x5 Defining x11 0 2 3 X2*X15 + X2*X3 Defining x3 x3 Defining x3 Defining x15 0 2 4 X2*X16 + X2*X4 Defining x4 x4 Defining x4 Defining x16 0 2 5 X2*X17 + X2*X5 Defining x5 x5 Defining x5 Defining x17 0 3 4 X3*X22 + X3*X4 Defining x4 x4 Defining x4 Defining x22 0 3 5 X3*X23 + X3*X5 Defining x5 x5 Defining x5 Defining x23 0 4 5 X4*X29 + X4*X5 Defining x5 x5 Defining x5 Defining x29 1 2 3 X8*X15 + X8*X9 Defining x9 x9 Defining x9 Defining x15 1 2 4 X8*X16 + X8*X10 Defining x10 x10 Defining x10 Defining x16 1 2 5 X8*X17 + X8*X11 Defining x11 x11 Defining x11 Defining x17 1 3 4 X9*X22 + X9*X10 Defining x10 x10 Defining x10 Defining x22 1 3 5 X9*X23 + X9*X11 Defining x11 x11 Defining x11 Defining x23 1 4 5 X10*X29 + X10*X11 Defining x11 x11 Defining x11 Defining x29 2 3 4 X15*X22 + X15*X16 Defining x16 x16 Defining x16 Defining x22 2 3 5 X15*X23 + X15*X17 Defining x17 x17 Defining x17 Defining x23 2 4 5 X16*X29 + X16*X17 Defining x17 x17 Defining x17 Defining x29 3 4 5 X22*X29 + X22*X23 Defining x23 x23 Defining x23 Defining x29

ProbM(X1*X8 + X1*X2)

Defining x1, x2, x8 -2*x1*x2*x8 + x1*x8 + x2*x8
load('ProbRecurse.py')
import numpy
import scipy.optimize
#from operator import mul

KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])
GIVEN = lambda A, B: ProbM(AND(A,B))/ProbM(B)
convert = lambda f: reduce(lambda a,b:a+b,[reduce(lambda a,b:a*b,[(m.degree(g)+1+g) for g in BPR.gens()]) for m in f.monomials()]) #self inverse convert poly in gens to poly in atoms (missing var x implies, interpret as ~x)
# convert monomial terms A^i*B^j*C^k --> (A+1+i)*(B+1+j)*(C*1*k) in Boolean Ring i,j,k=0,1

objects = 6 #rel_classes*types; print('objects = ' + str(objects))
BPR = BooleanPolynomialRing(objects^2,'X',order='degneglex')
BPR.inject_variables()
print(str(objects^2) + ' Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X' + str(objects^2-1) + ' = R[' + str(objects-1) + ',' + str(objects-1) + ']')
X = BPR.gens()
NewGens = str(X).lower()
NewGens = NewGens[1:len(NewGens) - 1]
PPR = PolynomialRing(QQ,objects^2,NewGens,order='degneglex')
PPR.inject_variables()
R = matrix(BPR, objects, objects, BPR.gens())
prems = [];reflex=[];symmtrc=[];trans2=[];direct=[]
for i in range(0, objects):    #reflexive
reflex.append(R[i,i])

for i in range(0,objects):
for j in range(0, objects):
#symmtrc.append(IF(R[i,j],R[j,i]))
direct.append(IF(R[i,j],NOT(R[j,i])))

for i in range(0,objects):    #symmetric
for j in range(i+1, objects):
symmtrc.append(IFF(R[i,j],R[j,i]))
symmtrc=list(set(symmtrc))

for j in range(0, objects):    #transitive
for i in range(0, objects):
for k in range(0, objects):
poly = IF(AND(R[min(i,j),max(i,j)],R[min(j,k),max(j,k)]),R[min(i,k),max(i,k)])
trans2.append(poly)
trans2 = list(set(trans2))

M=matrix([[125,25,5,1],[64,16,4,1],[27,9,3,1],[8,4,2,1]])#;M
var('t')
B=M.inverse()*vector([35,18,8,3])#;B
mypoly2=vector([t^3,t^2,t,1])*B;mypoly2
# mypoly2.subs(t=6);mypoly2.subs(t=7) #seems to give the number of polynomials in groebner basis for equiv relation with gen order degneglexx
prems=reflex+symmtrc+trans2
if BPR(0) in prems:
prems = list(set(prems)) # remove dup
premstrue = [p+1 for p in prems] #make true
REI = ideal(premstrue);# print('REI.gens() = ' + str(len(REI.gens())))
GB = REI.groebner_basis();GB
#compute groebner basis directly assuming equivalence relation
#correspond to trans2 objects=5: X1*(X2+X(2+1*5));X1*(X3+X(3+1*5));X1*(X4+X(4+1*5));X2*(X3+X(3+2*5));X2*(X4+X(4+2*5));X3*(X4+X(4+3*5));X7*(X8+X(8+1*5));X7*(X9+X(9+1*5));X8*(X9+X(9+2*5));X13*(X14+X(14+1*5)) then
#X2*(X1+X(2+1*5));X3*(X1+X(3+1*5));X4*(X1+X(4+1*5));X3*(X2+X(3+2*5));X4*(X2+X(4+2*5));X4*(X3+X(4+3*5));X8*(X7+X(8+1*5));X9*(X7+X(9+1*5));X9*(X8+X(9+2*5));X14*(X13+X(14+1*5)) that 20 polynomials
GB2 = []
PGB = []
for i in range(objects):
for j in range(i+1,objects):
for k in range(j+1,objects):
print(i);print(j);print(k)
print(R[i,j]*(R[i,k]+R[j,k]))
GB2.append(R[i,j]*(R[i,k]+R[j,k]))
print(PPR(ProbM(R[i,k])))
PGB.append(PPR(ProbM(R[i,k]))-PPR(ProbM(R[j,k])))
Prob(X1*X8 + X1*X2)

Defining X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20, X21, X22, X23, X24, X25, X26, X27, X28, X29, X30, X31, X32, X33, X34, X35 36 Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X35 = R[5,5] Defining x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31, x32, x33, x34, x35 t 1/3*t^3 - 1/2*t^2 + 7/6*t Polynomial Sequence with 61 Polynomials in 36 Variables 0 1 2 X1*X8 + X1*X2 Defining x2 x2 Defining x2 Defining x8 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 1 3 X1*X9 + X1*X3 Defining x3 x3 Defining x3 Defining x9 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 1 4 X1*X10 + X1*X4 Defining x4 x4 Defining x4 Defining x10 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 1 5 X1*X11 + X1*X5 Defining x5 x5 Defining x5 Defining x11 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 2 3 X2*X15 + X2*X3 Defining x3 x3 Defining x3 Defining x15 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 2 4 X2*X16 + X2*X4 Defining x4 x4 Defining x4 Defining x16 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 2 5 X2*X17 + X2*X5 Defining x5 x5 Defining x5 Defining x17 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 3 4 X3*X22 + X3*X4 Defining x4 x4 Defining x4 Defining x22 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 3 5 X3*X23 + X3*X5 Defining x5 x5 Defining x5 Defining x23 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 0 4 5 X4*X29 + X4*X5 Defining x5 x5 Defining x5 Defining x29 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 1 2 3 X8*X15 + X8*X9 Defining x9 x9 Defining x9 Defining x15 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 1 2 4 X8*X16 + X8*X10 Defining x10 x10 Defining x10 Defining x16 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 1 2 5 X8*X17 + X8*X11 Defining x11 x11 Defining x11 Defining x17 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 1 3 4 X9*X22 + X9*X10 Defining x10 x10 Defining x10 Defining x22 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 1 3 5 X9*X23 + X9*X11 Defining x11 x11 Defining x11 Defining x23 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 1 4 5 X10*X29 + X10*X11 Defining x11 x11 Defining x11 Defining x29 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 2 3 4 X15*X22 + X15*X16 Defining x16 x16 Defining x16 Defining x22 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 2 3 5 X15*X23 + X15*X17 Defining x17 x17 Defining x17 Defining x23 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 2 4 5 X16*X29 + X16*X17 Defining x17 x17 Defining x17 Defining x29 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8 3 4 5 X22*X29 + X22*X23 Defining x23 x23 Defining x23 Defining x29 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x2 + x1*x8
load('ProbRecurse.py')
#constraint
#mymat=matrix([[c.coefficient(g) for g in PPR.gens()] for c in constraint])
#PPR.gens()
#constraint=[ProbM(p)-1 for p in reflex]+symmtrc2+PGB
#print(constraint[1])
#print(constraint[1].coefficients())
#print(PPR(constraint[1]).coefficient(PPR.gens(1)))
myBP=X1*X8 + X1*X2
print(ProbM(X1*X8 + X1*X2))
print(Prob(X1*X8 + X1*X2))
print((Prob(X1*X8 + X1*X2)).monomials())
if type(myBP) == sage.rings.integer.Integer:
if myBP%2 == 0:
print(QQ(0))
else:
print(QQ(1))
myBP.variables()
mystr = str(myBP.variables()) # converting all uppers to lower and lowers to uppercase
newstr = ""
for l in mystr:
if l.isupper():
newstr = newstr + "1" # gotta mark them as upper or lower
else:
newstr = newstr + "0"
NewGens = ""
for ii in range(len(mystr)):
if newstr[ii] == "1":
NewGens = NewGens + mystr[ii].lower()
else:
NewGens = NewGens + mystr[ii].upper() # if boolean variable XcGGt1u then, prob var is xCggT1U etc.
NewGens = NewGens[1:len(NewGens) - 1] # generators for probabilty polynomial assuming boolean generators are independent
if len(myBP.variables()) == 1:
NewGens = NewGens.replace(',','')

CPR = PolynomialRing(QQ,len(myBP.variables()),NewGens.replace(', ',','),order='invlex') # needs order to be lex to work
CPR.inject_variables()
nn = myBP.nvariables() # number of generators in boolean polynomial
BPvars = myBP.variables()
BPmons = myBP.monomials()
atoms = 2**nn # number of atoms in corresponding boolean algebra
# should abort this with error, if too many variables in polynomial
Mvect = [0]*atoms
for jj in range(len(BPmons)):
mm = 0
for kk in range(len(BPmons[jj].variables())):
mm = mm + 2**BPvars.index(BPmons[jj].variables()[kk])
Mvect[mm] = 1
Mvect.reverse()
# Mvect order for instance is (X111,X110,X101,X100,X011,X010,X001,X000) has two nonequiv intepretations
# as monomials or atoms example was for three variables
Mvect = vector(Mvect)
MGinv = matrix(QQ,[1])
for ii in range(1,nn+1):
MGinv = block_matrix(QQ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False) # creating inverse of multigrade operator 'AND' matrix MG
# modulo 2 it's an involution however...
# MG and MGinv are upper triangular, MG is like Sierpinski triangle see
Avect = (MGinv*Mvect)%2
Pvect = MGinv*Avect # convert from boolean ring (mons) to boolean algebra (atoms) then to probability polynomial
# ((MGinv*Mvect)%2) is the atoms vector equivalent to the monnomials vector Mvect
# Pvect are the coefficients of prob poly in CPR
newpoly = CPR(1)
for ii in range(0,len(CPR.gens())):
newpoly = newpoly*(CPR(1) + CPR.gen(ii)) # to create all possible monomials in CPR
print(newpoly)
print(newpoly.monomials())
print(Pvect)
print(Pvect*vector(newpoly.monomials())) # inner product of the two vectors to create prob poly
print(ProbM(X1*X8 + X8*X2))
print(ProbM(X1*X2 + X8*X2))
print(Prob(X1*X2 + X8*X2))
print(ProbM(1+(X1*(1+X2))))
print(Prob(1+(X1*(1+X2))))

Defining x1, x2, x8 -2*x1*x2*x8 + x1*x8 + x1*x2 Defining x1, x2, x8 -2*x1*x2*x8 + x1*x8 + x1*x2 Defining x1, x2, x8 [x1*x2*x8, x1*x8, x1*x2] (X1, X2, X8) Defining x1, x2, x8 x1*x2*x8 + x2*x8 + x1*x8 + x8 + x1*x2 + x2 + x1 + 1 [x1*x2*x8, x2*x8, x1*x8, x8, x1*x2, x2, x1, 1] (-2, 0, 1, 0, 1, 0, 0, 0) -2*x1*x2*x8 + x1*x8 + x1*x2 Defining x1, x2, x8 -2*x1*x2*x8 + x2*x8 + x1*x8 Defining x1, x2, x8 -2*x1*x2*x8 + x2*x8 + x1*x2 Defining x1, x2, x8 -2*x1*x2*x8 + x2*x8 + x1*x2 Defining x1, x2 x1*x2 - x1 + 1 Defining x1, x2 x1*x2 - x1 + 1
KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])
GIVEN = lambda A, B: Prob(AND(A,B))/Prob(B)
convert = lambda f: reduce(lambda a,b:a+b,[reduce(lambda a,b:a*b,[(m.degree(g)+1+g) for g in BPR.gens()]) for m in f.monomials()]) #self inverse convert poly in gens to poly in atoms (missing var x implies, interpret as ~x)
# convert monomial terms A^i*B^j*C^k --> (A+1+i)*(B+1+j)*(C+1+k) in Boolean Ring i,j,k=0,1

rel_classes = 4
types = 3
objects = rel_classes*types; print('objects = ' + str(objects))
BPR = BooleanPolynomialRing(types*rel_classes*rel_classes,'X',order='degneglex')
BPR.inject_variables()
#print(str(objects^2) + ' Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X' + str(objects^2-1) + ' = R[' + str(objects-1) + ',' + str(objects-1) + ']')
R = matrix(BPR, objects, objects)
trans=[];exactlyone=[];clues = [];

for i in range(0, objects):    #reflexive
R[i,i] = BPR(1)

nn=-1
for i in range(0,types):    #symmetric
for j in range(i+1, types):
nn=nn+1
mm=0
for k in range(0,rel_classes):
for l in range(0,rel_classes):
R[i*rel_classes+k,rel_classes*j+l] = BPR.gens()[nn*16+mm]
R[rel_classes*j+l,i*rel_classes+k] = BPR.gens()[nn*16+mm]
mm=mm+1

for i in range(objects):     #transitive truth assertion built into these which assume symmtrc and reflex in place
for j in range(objects):
for k in range(objects):
trans.append(IF(R[i,j]*R[j,k],R[i,k]))

trans = list(set(trans))
trans.remove(1)

for i in range(0,types):
for j in range(i,types):
for k in range(0,rel_classes):
poly = UNIQUE([R[i*rel_classes+k,j*rel_classes+b] for b in range(0,rel_classes)])
exactlyone.append(poly)
poly = UNIQUE([R[j*rel_classes+b,i*rel_classes+k] for b in range(0,rel_classes)])
exactlyone.append(poly)

exactlyone=set(list(exactlyone))
exactlyone.remove(BPR(1))
exactlyone = list(exactlyone)

clues.append(UNIQUE([R[9,1],R[9,6]]))    # clue 1
clues.append(UNIQUE([AND(R[4,11],R[6,10]), AND(R[4,10],R[6,9]), AND(R[4,9],R[6,8])]))    # clue 2
clues.append(UNIQUE([AND(R[11,6],R[9,0]), AND(R[11,0],R[9,6])]))    # clue 3
clues.append(UNIQUE([AND(R[0,11],R[2,4]), AND(R[0,4],R[2,11])]))    # clue 4
clues.append(UNIQUE([AND(R[5,11],R[4,10]), AND(R[5,11],R[4,9]), AND(R[5,11],R[4,8]), AND(R[5,10],R[4,9]), AND(R[5,10],R[4,8]), AND(R[5,9],R[4,8])]))    # clue 5

prems = trans + exactlyone + clues
prems = [p+1 for p in prems]
REI=ideal(prems)
QR = BPR.quotient_ring(REI)
solnmat=matrix([[QR(R[i,j]) for i in range(objects)] for j in range(objects)])
sd=[j*rel_classes for j in range(1,types)]
solnmat.subdivide(sd,sd)
print(solnmat);print(' ')

objects = 12 Defining X0, X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, X16, X17, X18, X19, X20, X21, X22, X23, X24, X25, X26, X27, X28, X29, X30, X31, X32, X33, X34, X35, X36, X37, X38, X39, X40, X41, X42, X43, X44, X45, X46, X47 [1 0 0 0|0 1 0 0|0 0 0 1] [0 1 0 0|0 0 0 1|1 0 0 0] [0 0 1 0|1 0 0 0|0 0 1 0] [0 0 0 1|0 0 1 0|0 1 0 0] [-------+-------+-------] [0 0 1 0|1 0 0 0|0 0 1 0] [1 0 0 0|0 1 0 0|0 0 0 1] [0 0 0 1|0 0 1 0|0 1 0 0] [0 1 0 0|0 0 0 1|1 0 0 0] [-------+-------+-------] [0 1 0 0|0 0 0 1|1 0 0 0] [0 0 0 1|0 0 1 0|0 1 0 0] [0 0 1 0|1 0 0 0|0 0 1 0] [1 0 0 0|0 1 0 0|0 0 0 1]
KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])

rel_classes = 4
types = 3
objects = rel_classes*types; print('objects = ' + str(objects))
BPR = BooleanPolynomialRing(objects^2,'X',order='degrevlex')
print('types = ' + str(types))
print('rel_classes = ' + str(rel_classes))
print(str(objects^2) + ' Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X' + str(objects^2-1) + ' = R[' + str(objects-1) + ',' + str(objects-1) + ']')
#X = BPR.gens()
R = matrix(BPR, objects, objects, BPR.gens())
clues = []
for i in range(0,objects):    #symmetric
for j in range(0, objects):
poly = IF(R[i,j],R[j,i])
clues.append(poly)
for j in range(0, objects):    #transitive
for i in range(0, objects):
for k in range(0, objects):
poly = IF(AND(R[i,j],R[j,k]),R[i,k])
clues.append(poly)
clues.append(UNIQUE([R[9,1],R[9,6]]))    # clue 1
clues.append(UNIQUE([AND(R[4,11],R[6,10]), AND(R[4,10],R[6,9]), AND(R[4,9],R[6,8])]))    # clue 2
clues.append(UNIQUE([AND(R[11,6],R[9,0]), AND(R[11,0],R[9,6])]))    # clue 3
clues.append(UNIQUE([AND(R[0,11],R[2,4]), AND(R[0,4],R[2,11])]))    # clue 4
clues.append(UNIQUE([AND(R[5,11],R[4,10]), AND(R[5,11],R[4,9]), AND(R[5,11],R[4,8]), AND(R[5,10],R[4,9]), AND(R[5,10],R[4,8]), AND(R[5,9],R[4,8])]))    # clue 5
clues = [B + 1 for B in clues] #set to true
clues = list(set(clues)) #remove duplicates
if 0 in clues: clues.remove(0)
REI = ideal(clues); print('REI.gens() = ' + str(len(REI.gens())))
GB = REI.groebner_basis()
QR = BPR.quotient_ring(REI)
QRG = QR.gens()
quo = len(QRG) // objects
rem = len(QRG) % objects
if len(QRG) == objects^2:
A = matrix(BPR, objects, objects, 0)
for i in range(quo):
A[i,0:objects] = matrix(QRG[i*objects: (i + 1)*objects])
print('A 1st pass')
clues2 = []
for i in range(types-1):    # more hidden assumptions -- each object is related to exactly one oject in each other type
for j in range(i+1,types):
A2=A[i*rel_classes:(i+1)*rel_classes,j*rel_classes:(j+1)*rel_classes]
print('type '+ str(i) + ' to type ' + str(j));print(A2);print(' ')
for k in range(rel_classes):
poly = BPR(UNIQUE([A2[k,l] for l in range(rel_classes)]))
clues2.append(poly)
poly = BPR(UNIQUE([A2[l,k] for l in range(rel_classes)]))
clues2.append(poly)
clues2 = [B + 1 for B in clues2] #set to true
clues2 = list(set(clues2)) #remove duplicates
if 0 in clues2: clues2.remove(0)
stuff=[clues2[l].variables() for l in range(len(clues2))]
stuff2 = set(stuff[0])
for i in range(1,len(stuff)):
stuff2 = stuff2.union(set(stuff[i]))
print('number of variables left to be determined: ' + str(len(stuff2)))
clues = clues + clues2
REI = ideal(clues); print('REI.gens() = ' + str(len(REI.gens())))
GB = REI.groebner_basis()
QR = BPR.quotient_ring(REI)
QRG = QR.gens()
quo = len(QRG) // objects
rem = len(QRG) % objects
if len(QRG) == objects^2:
A3 = matrix(BPR, objects, objects, 0)
for i in range(quo):
A3[i,0:objects] = matrix(QRG[i*objects: (i + 1)*objects])
print('A 2nd pass')
for i in range(types-1):
for j in range(i+1,types):
A2=A3[i*rel_classes:(i+1)*rel_classes,j*rel_classes:(j+1)*rel_classes]
print('type '+ str(i) + ' to type ' + str(j));print(A2);print(' ')

objects = 12 types = 3 rel_classes = 4 144 Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X143 = R[11,11] REI.gens() = 1589 A 1st pass type 0 to type 1 [ X138 X124 + X126 + X128 0 X90 + X115 + X139] [ X49 X61 0 X85] [ X138 + 1 X124 + X125 + X126 + X129 + 1 X124 + X129 + 1 X88 + X90 + X115] [ X51 X63 X75 X87] type 0 to type 2 [ X138 + X140 X138 X126 + X129 + X142 X138 + 1] [ X97 X138 X121 X133] [X124 + X126 + X128 + X129 + 1 X124 + X129 + X138 + 1 X124 + X126 + X129 X138] [ X99 X111 X123 X135] type 1 to type 2 [X124 + X126 + X128 + X129 + X138 + 1 X124 + X129 + 1 X124 0] [ X101 X125 + X129 + X137 + 1 X125 X137] [ X116 + X138 X138 + 1 X126 X138] [ X103 X115 X127 X139] number of variables left to be determined: 31 REI.gens() = 1613 A 2nd pass type 0 to type 1 [0 1 0 0] [0 0 0 1] [1 0 0 0] [0 0 1 0] type 0 to type 2 [0 0 0 1] [1 0 0 0] [0 0 1 0] [0 1 0 0] type 1 to type 2 [0 0 1 0] [0 0 0 1] [0 1 0 0] [1 0 0 0]
KDELTA = lambda A, B: A.parent(A == B)
NOT = lambda A: A.parent(1) + A
XOR = lambda A, B: A + B
OR = lambda A, B: A + B + A*B
IF = lambda A, B: A.parent(1) + A + A*B
IFF = lambda A, B: A.parent(1) + A + B
AND = lambda A, B: A*B
NAND = lambda A, B: A.parent(1) + A*B
NOR = lambda A, B: A.parent(1) + A + B + A*B
FORALL = lambda mylist: reduce(AND, mylist)
EXISTS = lambda mylist: reduce(OR, mylist)
UNIQUE = lambda mylist: reduce(XOR,[FORALL([A + NOT(KDELTA(A,B)) for A in mylist]) for B in mylist])

rel_classes = 6
types = 4
objects = rel_classes*types; print('objects = ' + str(objects))
print('types = ' + str(types))
print('rel_classes = ' + str(rel_classes))
print(str(objects^2) + ' Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X' + str(objects^2-1) + ' = R[' + str(objects-1) + ',' + str(objects-1) + ']')
BPR = BooleanPolynomialRing(objects^2,'X',order='degneglex')
#X = BPR.gens()
R = matrix(BPR, objects, objects, BPR.gens())
clues = []
for i in range(0, objects):  #reflexive
poly = R[i,i]
clues.append(poly)
for i in range(0,objects):    #symmetric
for j in range(0, objects):
poly = IF(R[i,j],R[j,i])
clues.append(poly)
for j in range(0, objects):    #transitive
for i in range(0, objects):
for k in range(0, objects):
poly = IF(AND(R[i,j],R[j,k]),R[i,k])
clues.append(poly)
for i in range(0, types):    # distinct objects of the same type are not related
iobj = i*rel_classes
for j in range(0, rel_classes):
iobjj = iobj + j
for k in range(0, rel_classes):
if j!=k:
iobjk = iobj + k
poly = NOT(R[iobjj,iobjk])
clues.append(poly)
#for i in range(types):    # each object is related to exactly one oject of each other type
#    for j in range(rel_classes):
#        for k in range(types):
#            if k != i:
#                poly = UNIQUE([R[i*rel_classes + j, k*rel_classes + l] for l in range(rel_classes)])
#                clues.append(poly)
#for i in range(types-1):    # another way to do the same but placing either here exceeds size limitations for quotient ring and groebner basis algorithms
#    for j in range(rel_classes):
#        for k in range(i+1,types):
#            poly = UNIQUE([R[i*rel_classes + j, k*rel_classes + l] for l in range(rel_classes)])
#            clues.append(poly)
#            poly = UNIQUE([R[i*rel_classes + l, k*rel_classes + j] for l in range(rel_classes)])
#            clues.append(poly)
clues.append(UNIQUE([AND(R[11,22],R[16,1]), AND(R[11,1],R[16,22])]))    # clue 1
clues.append(UNIQUE([AND(R[1,18],R[15,19]), AND(R[1,18],R[15,20]), AND(R[1,18],R[15,21]), AND(R[1,18],R[15,22]), AND(R[1,18],R[15,23]),    # clue 2
AND(R[1,19],R[15,20]), AND(R[1,19],R[15,21]), AND(R[1,19],R[15,22]), AND(R[1,19],R[15,23]),
AND(R[1,20],R[15,21]), AND(R[1,20],R[15,22]), AND(R[1,20],R[15,23]),
AND(R[1,21],R[15,22]), AND(R[1,21],R[15,23]),
AND(R[1,22],R[15,23])]))
clues.append(UNIQUE([AND(R[13,23],R[9,21]), AND(R[13,22],R[9,20]), AND(R[13,21],R[9,19]), AND(R[13,20],R[9,18])]))    # clue 3
clues.append(R[5,11])    # clue 4
clues.append(NOT(R[4,23]))    # clue 5
clues.append(NOT(R[1,10]))    # clue 6.1
clues.append(NOT(R[1,14]))    # clue 6.2
clues.append(NOT(R[1,17]))    # clue 6.3
clues.append(NOT(R[1,9]))    # clue 6.4
clues.append(NOT(R[1,8]))    # clue 6.5
clues.append(NOT(R[10,14]))    # clue 6.6
clues.append(NOT(R[10,17]))    # clue 6.7
clues.append(NOT(R[14,9]))    # clue 6.8
clues.append(NOT(R[14,8]))    # clue 6.9
clues.append(NOT(R[17,9]))    # clue 6.10
clues.append(NOT(R[17,8]))    # clue 6.11
clues.append(R[21,6])    # clue 7
clues.append(UNIQUE([R[0,8], R[0,7]]))    # clue 8
clues.append(NOT(R[13,10]))    # clue 9
clues.append(NOT(R[3,14]))    # clue 10
clues.append(UNIQUE([AND(R[16,23],R[3,21]), AND(R[16,22],R[3,20]), AND(R[16,21],R[3,19]), AND(R[16,20],R[3,18])]))    # clue 11
clues = [B + 1 for B in clues] #set to true
clues = list(set(clues)) #remove duplicate
if 0 in clues: clues.remove(0)
REI = ideal(clues); print('REI.gens() = ' + str(len(REI.gens())))
GB = REI.groebner_basis()
QR = BPR.quotient_ring(REI)
QRG = QR.gens()
quo = len(QRG) // objects
rem = len(QRG) % objects
if len(QRG) == objects^2:
A = matrix(BPR, objects, objects, 0)
for i in range(quo):
A[i,0:objects] = matrix(QRG[i*objects: (i + 1)*objects])
print('A 1st pass')
clues2 = []
for i in range(types-1):    # more hidden assumptions -- each object is related to exactly one object in each other type
for j in range(i+1,types):
A2=A[i*rel_classes:(i+1)*rel_classes,j*rel_classes:(j+1)*rel_classes]
print('type ' + str(i) + ' to type ' + str(j));print(A2);print(' ')
for k in range(rel_classes):
poly = BPR(UNIQUE([A2[k,l] for l in range(rel_classes)]))
clues2.append(poly)
poly = BPR(UNIQUE([A2[l,k] for l in range(rel_classes)]))
clues2.append(poly)
clues2 = [B + 1 for B in clues2] #set to true
clues2 = list(set(clues2)) #remove duplicates
if 0 in clues2: clues2.remove(0)
stuff=[clues2[l].variables() for l in range(len(clues2))]
stuff2 = set(stuff[0])
for i in range(1,len(stuff)):
stuff2 = stuff2.union(set(stuff[i]))
print('number of variables left to be determined: ' + str(len(stuff2)))
clues = clues + clues2
REI = ideal(clues); print('REI.gens() = ' + str(len(REI.gens())))
GB = REI.groebner_basis()
QR = BPR.quotient_ring(REI)
QRG = QR.gens()
quo = len(QRG) // objects
rem = len(QRG) % objects
if len(QRG) == objects^2:
A3 = matrix(BPR, objects, objects, 0)
for i in range(quo):
A3[i,0:objects] = matrix(QRG[i*objects: (i + 1)*objects])
print('A 2nd pass')
for i in range(types-1):
for j in range(i+1,types):
A2=A3[i*rel_classes:(i+1)*rel_classes,j*rel_classes:(j+1)*rel_classes]
print('type ' + str(i) + ' to type ' + str(j));print(A2);print(' ')

objects = 24 types = 4 rel_classes = 6 576 Free Boolean Algebra generators: X0 = R[0,0], X1 = R[0,1], X2 = R[0,2],..., X575 = R[23,23] REI.gens() = 13413 A 1st pass type 0 to type 1 [ 0 X7 X7 + 1 0 0 0] [ X30 X31 0 0 0 0] [ X54 X55 X56 X57 X58 0] [ 0 X79 X80 0 X82 0] [ X102 X103 X104 X105 X106 0] [ 0 0 0 0 0 1] type 0 to type 2 [ X12 X13 X14 X15 0 X17] [ 0 0 0 0 1 0] [ X60 X61 X62 X63 0 X65] [ X84 0 0 0 0 X89] [X108 X109 X110 0 0 X113] [X132 X133 X134 X135 0 X137] type 0 to type 3 [ X18 0 X13 0 0 X23] [ 0 0 X30 + 1 X30 0 0] [ X66 X67 X68 X54 0 X71] [X30 + 1 X30 0 0 0 0] [ X114 X115 X116 X102 0 0] [ 0 0 0 0 1 0] type 1 to type 2 [ 0 X30 + 1 0 0 X30 0] [ X180 X181 X182 X183 X31 X185] [ X204 X205 0 X207 0 0] [ X228 0 0 0 0 0] [ X252 0 0 X255 0 0] [ X132 X133 X134 X135 0 X137] type 1 to type 3 [ 0 0 0 1 0 0] [ X186 X187 X181 + X31 0 0 X191] [ X210 X211 X205 0 0 X215] [X133 + X30 X30 + 1 X133 0 0 0] [ X258 X259 0 0 0 X263] [ 0 0 0 0 1 0] type 2 to type 3 [ X306 X307 X308 0 X132 X311] [ 0 0 X133 + X30 X30 + 1 X133 0] [ X354 0 0 0 X134 X359] [ 0 0 0 0 X135 X135 + 1] [ 0 0 X30 + 1 X30 0 0] [ X426 X427 0 0 X137 X431] number of variables left to be determined: 75 REI.gens() = 13463 A 2nd pass type 0 to type 1 [0 0 1 0 0 0] [1 0 0 0 0 0] [0 0 0 0 1 0] [0 1 0 0 0 0] [0 0 0 1 0 0] [0 0 0 0 0 1] type 0 to type 2 [0 1 0 0 0 0] [0 0 0 0 1 0] [0 0 0 1 0 0] [0 0 0 0 0 1] [1 0 0 0 0 0] [0 0 1 0 0 0] type 0 to type 3 [0 0 1 0 0 0] [0 0 0 1 0 0] [0 0 0 0 0 1] [0 1 0 0 0 0] [1 0 0 0 0 0] [0 0 0 0 1 0] type 1 to type 2 [0 0 0 0 1 0] [0 0 0 0 0 1] [0 1 0 0 0 0] [1 0 0 0 0 0] [0 0 0 1 0 0] [0 0 1 0 0 0] type 1 to type 3 [0 0 0 1 0 0] [0 1 0 0 0 0] [0 0 1 0 0 0] [1 0 0 0 0 0] [0 0 0 0 0 1] [0 0 0 0 1 0] type 2 to type 3 [1 0 0 0 0 0] [0 0 1 0 0 0] [0 0 0 0 1 0] [0 0 0 0 0 1] [0 0 0 1 0 0] [0 1 0 0 0 0]
print("percolation")
#list(range(0,4))
N=5

print(PP)
print()
for ii in range(1,N-1):
PPL.append(PP)
print(PP)
print()

for ii in range(1,4):
for jj in range(N):
for kk in range(N):
if jj != kk:
if NSP[jj,kk] == 0:
if (PPL[ii])[jj,kk] != 0:
SPL[jj,kk] = ii+1
NSP[jj,kk] = (PPL[ii])[jj,kk]

print("Number of Shortest Paths")
print(NSP)
print()
print("Shortest Path Length")
print(SPL)
print()
DB = []
nn = -1
for ii in range(N):
for jj in range(N):
if jj != ii:
for kk in range(N):
LL = []
if kk != ii and kk != jj:
nn = nn + 1
LL = [ii, jj, kk, NSP[ii,kk], SPL[ii,kk], NSP[ii, jj], SPL[ii, jj], NSP[jj, kk], SPL[jj, kk], NSP[ii, jj]*NSP[jj, kk], SPL[ii, jj] + SPL[jj, kk]]
DB.insert(nn,LL)
#print(LL)

for db in DB:
if db[10] == db[4]:
db.append(db[9])
else:
db.append(0)

print(matrix(DB))
print()
print(len(DB))


percolation [0 1 1 1 1] [1 0 1 0 0] [1 1 0 1 1] [1 0 1 0 1] [1 0 1 1 0] [4 1 3 2 2] [1 2 1 2 2] [3 1 4 2 2] [2 2 2 3 2] [2 2 2 2 3] [8 7 9 9 9] [7 2 7 4 4] [9 7 8 9 9] [9 4 9 6 7] [9 4 9 7 6] [34 17 33 26 26] [17 14 17 18 18] [33 17 34 26 26] [26 18 26 25 24] [26 18 26 24 25] Number of Shortest Paths [0 1 1 1 1] [1 0 1 2 2] [1 1 0 1 1] [1 2 1 0 1] [1 2 1 1 0] Shortest Path Length [0 1 1 1 1] [1 0 1 2 2] [1 1 0 1 1] [1 2 1 0 1] [1 2 1 1 0] [0 1 2 1 1 1 1 1 1 1 2 0] [0 1 3 1 1 1 1 2 2 2 3 0] [0 1 4 1 1 1 1 2 2 2 3 0] [0 2 1 1 1 1 1 1 1 1 2 0] [0 2 3 1 1 1 1 1 1 1 2 0] [0 2 4 1 1 1 1 1 1 1 2 0] [0 3 1 1 1 1 1 2 2 2 3 0] [0 3 2 1 1 1 1 1 1 1 2 0] [0 3 4 1 1 1 1 1 1 1 2 0] [0 4 1 1 1 1 1 2 2 2 3 0] [0 4 2 1 1 1 1 1 1 1 2 0] [0 4 3 1 1 1 1 1 1 1 2 0] [1 0 2 1 1 1 1 1 1 1 2 0] [1 0 3 2 2 1 1 1 1 1 2 1] [1 0 4 2 2 1 1 1 1 1 2 1] [1 2 0 1 1 1 1 1 1 1 2 0] [1 2 3 2 2 1 1 1 1 1 2 1] [1 2 4 2 2 1 1 1 1 1 2 1] [1 3 0 1 1 2 2 1 1 2 3 0] [1 3 2 1 1 2 2 1 1 2 3 0] [1 3 4 2 2 2 2 1 1 2 3 0] [1 4 0 1 1 2 2 1 1 2 3 0] [1 4 2 1 1 2 2 1 1 2 3 0] [1 4 3 2 2 2 2 1 1 2 3 0] [2 0 1 1 1 1 1 1 1 1 2 0] [2 0 3 1 1 1 1 1 1 1 2 0] [2 0 4 1 1 1 1 1 1 1 2 0] [2 1 0 1 1 1 1 1 1 1 2 0] [2 1 3 1 1 1 1 2 2 2 3 0] [2 1 4 1 1 1 1 2 2 2 3 0] [2 3 0 1 1 1 1 1 1 1 2 0] [2 3 1 1 1 1 1 2 2 2 3 0] [2 3 4 1 1 1 1 1 1 1 2 0] [2 4 0 1 1 1 1 1 1 1 2 0] [2 4 1 1 1 1 1 2 2 2 3 0] [2 4 3 1 1 1 1 1 1 1 2 0] [3 0 1 2 2 1 1 1 1 1 2 1] [3 0 2 1 1 1 1 1 1 1 2 0] [3 0 4 1 1 1 1 1 1 1 2 0] [3 1 0 1 1 2 2 1 1 2 3 0] [3 1 2 1 1 2 2 1 1 2 3 0] [3 1 4 1 1 2 2 2 2 4 4 0] [3 2 0 1 1 1 1 1 1 1 2 0] [3 2 1 2 2 1 1 1 1 1 2 1] [3 2 4 1 1 1 1 1 1 1 2 0] [3 4 0 1 1 1 1 1 1 1 2 0] [3 4 1 2 2 1 1 2 2 2 3 0] [3 4 2 1 1 1 1 1 1 1 2 0] [4 0 1 2 2 1 1 1 1 1 2 1] [4 0 2 1 1 1 1 1 1 1 2 0] [4 0 3 1 1 1 1 1 1 1 2 0] [4 1 0 1 1 2 2 1 1 2 3 0] [4 1 2 1 1 2 2 1 1 2 3 0] [4 1 3 1 1 2 2 2 2 4 4 0] [4 2 0 1 1 1 1 1 1 1 2 0] [4 2 1 2 2 1 1 1 1 1 2 1] [4 2 3 1 1 1 1 1 1 1 2 0] [4 3 0 1 1 1 1 1 1 1 2 0] [4 3 1 2 2 1 1 2 2 2 3 0] [4 3 2 1 1 1 1 1 1 1 2 0] 60
for db in DB:
if db[4] == 1:
db.append(0)
else:
db.append(1) #stub

print(matrix(DB))


[0 1 2 1 1 0] [0 1 3 1 1 0] [0 1 4 1 1 0] [0 2 1 1 1 0] [0 2 3 1 1 0] [0 2 4 1 1 0] [0 3 1 1 1 0] [0 3 2 1 1 0] [0 3 4 1 1 0] [0 4 1 1 1 0] [0 4 2 1 1 0] [0 4 3 1 1 0] [1 0 2 1 1 0] [1 0 3 2 2 1] [1 0 4 2 2 1] [1 2 0 1 1 0] [1 2 3 2 2 1] [1 2 4 2 2 1] [1 3 0 1 1 0] [1 3 2 1 1 0] [1 3 4 2 2 1] [1 4 0 1 1 0] [1 4 2 1 1 0] [1 4 3 2 2 1] [2 0 1 1 1 0] [2 0 3 1 1 0] [2 0 4 1 1 0] [2 1 0 1 1 0] [2 1 3 1 1 0] [2 1 4 1 1 0] [2 3 0 1 1 0] [2 3 1 1 1 0] [2 3 4 1 1 0] [2 4 0 1 1 0] [2 4 1 1 1 0] [2 4 3 1 1 0] [3 0 1 2 2 1] [3 0 2 1 1 0] [3 0 4 1 1 0] [3 1 0 1 1 0] [3 1 2 1 1 0] [3 1 4 1 1 0] [3 2 0 1 1 0] [3 2 1 2 2 1] [3 2 4 1 1 0] [3 4 0 1 1 0] [3 4 1 2 2 1] [3 4 2 1 1 0] [4 0 1 2 2 1] [4 0 2 1 1 0] [4 0 3 1 1 0] [4 1 0 1 1 0] [4 1 2 1 1 0] [4 1 3 1 1 0] [4 2 0 1 1 0] [4 2 1 2 2 1] [4 2 3 1 1 0] [4 3 0 1 1 0] [4 3 1 2 2 1] [4 3 2 1 1 0]