Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download

Algebra to solve logic puzzles Convert boolean polynomials to probability

292 views
License: GPL3
ubuntu2204
Tweet
Share
Share
/ext/sage/10.6/local/var/lib/sage/venv-python3.12.5/lib/python3.12/site-packages/scikits/__init__.py:1: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html __import__("pkg_resources").declare_namespace(__name__) /ext/sage/10.6/local/var/lib/sage/venv-python3.12.5/lib/python3.12/site-packages/scikits/__init__.py:1: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('scikits')`. Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages __import__("pkg_resources").declare_namespace(__name__) 
load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
import numpy
import scipy.optimize
#from operator import add
#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:
    print("explicit contradiction in premises")
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')
load('Prob.py')
load('ProbM.py')
import numpy
import scipy.optimize
#from operator import add
#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:
    print("explicit contradiction in premises")
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')
load('Prob.py')
load('ProbM.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()
aa = 2**nn # number of atoms in corresponding boolean algebra
    # should abort this with error, if too many variables in polynomial
Mvect = [0]*aa
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
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
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))))
Error in lines 4-4 Traceback (most recent call last): File "/cocalc/lib/python3.11/site-packages/smc_sagews/sage_server.py", line 1250, in execute exec( File "", line 1, in <module> NameError: name 'X1' is not defined 
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
ADJ = zero_matrix(ZZ, N, N)
ADJ = matrix([[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]])

PP = ADJ
print(PP)
print()
PPL = [ADJ]
for ii in range(1,N-1):
    PP=PP*ADJ
    PPL.append(PP)
    print(PP)
    print()

NSP = copy(ADJ)
SPL = copy(ADJ)
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] 
Prob(X1*X8 + X1*X2);
ProbM(X1*X8 + X1*X2);
ProbM(X1*X8 + X1);
nn=3;
MGinv = matrix(QQ,[1]);
for ii in range(1,nn+1):
    MGinv = block_matrix(QQ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False)
print(MGinv);
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, x8 -x1*x8 + x1 [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 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])

GPR = BooleanPolynomialRing(3,'X',order='degneglex')
GPR.gens()
n=len(GPR.gens())
n
GPR.order
GPR.inject_variables()
#Prob(X1*X2 + X1*X0)
a1=ProbM2(X1*X2 + X1*X0)
len(a1)
#P1=ProbM(X1*X2 + X1)
#P1
#P1.coefficients()

mybp=(X0+1)*(X1+1)
mybp=FORALL([NOT(A) for A in GPR.gens()])
mybp
mybp.monomials()
mybp.variables()
#Pmybp=ProbM(mybp)
#Pmybp
#Pmybp.coefficients()

nn=len(GPR.gens());
MGinv = matrix(QQ,[1]);
for ii in range(1,nn+1):
    MGinv = block_matrix(QQ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False)
#print(MGinv);
#(x0+1)*(x1+1)*(x2+1)
CPR = PolynomialRing(QQ,len(GPR.gens()),'x',order='invlex')
CPR.inject_variables()
(FORALL([(1 + x) for x in CPR.gens()])).monomials()
#answer = ProbM2(X1*X2 + X1*X0)
(X0, X1, X2) 3 <built-in method order of sage.rings.polynomial.pbori.pbori.BooleanPolynomialRing object at 0x7f811f893040> Defining X0, X1, X2 Defining x0, x1, x2 2 X0*X1*X2 + X1*X2 + X0*X2 + X0*X1 + X2 + X1 + X0 + 1 [X0*X1*X2, X1*X2, X0*X2, X0*X1, X2, X1, X0, 1] (X0, X1, X2) Defining x0, x1, x2 [x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1] 
load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
load('ProbM2.py')
load('ProbSceme.py')
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])
GPR = BooleanPolynomialRing(3,'X',order='lp')
GPR.gens()
GPR.order
GPR.inject_variables()
Prob(X0*X2 + X1*X0+X0)
ProbM(X0*X2 + X1*X0+X0)
ProbM2(X0*X2 + X1*X0+X0)
#a1=ProbM2(X1*X2+X1)
#len(a1)
#a1[0]
#a1[1]
#mybp=(1+X0)*(1+X1)*(1+X2)#(X1*X2 + X1*X0)
MyBP=FORALL([(GPR(1) + x) for x in GPR.gens()])
MyBP
MyBP.parent()
MyBP.parent().gens()
len(MyBP.parent().gens())
MyBPMons=MyBP.monomials()
MyBP.monomials()
2**len(MyBP.parent().gens())
nn=len(GPR.gens());
MG=matrix(ZZ,[1])
MGinv = matrix(ZZ,[1]);
for ii in range(1,nn+1):
    MG = block_matrix(ZZ,2,2,[MG,MG,0,MG],subdivide=False)
    MGinv = block_matrix(ZZ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False)
print(MGinv);
#(x0+1)*(x1+1)*(x2+1)
CPR = PolynomialRing(QQ,len(GPR.gens()),'x',order='invlex')
CPR.inject_variables()
MyMons=(FORALL([(CPR(1) + x) for x in CPR.gens()])).monomials()
MyMons
#Mvect = vector([0,0,0,0,1,1,0,0])
#Avect = (MGinv*Mvect)%2
#Pvect = MGinv*Avect
#Pvect
Mvect = vector([0,1,0,0,1,0,0,0])
Avect = (MGinv*Mvect)%2
Pvect = MGinv*Avect
Mvect
Avect
Pvect
mystr='XcGGt1u'
ii=0
if mystr[ii]==mystr[ii].upper():
    mystr[ii].lower()
else:
    mystr[ii].upper()
mystr
#MG=MGinv%2
MG #*
(MyBP.monomials())[0]
MonVec= vector([x for x in MyBP.monomials()])
MonVec
MonVec2= vector([x for x in MyMons])
MonVec2
#MG*MonVec
MonVec*MG
(MonVec*MG)*MG
(X0, X1, X2) <built-in method order of sage.rings.polynomial.pbori.pbori.BooleanPolynomialRing object at 0x7fb60d09d9c0> Defining X0, X1, X2 Defining x0, x1, x2 2*x0*x1*x2 - x0*x2 - x0*x1 + x0 Defining x0, x1, x2 2*x0*x1*x2 - x0*x1 - x0*x2 + x0 Defining x0, x1, x2 ((2, 0, -1, 0, -1, 0, 1, 0), (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1)) X0*X1*X2 + X0*X1 + X0*X2 + X0 + X1*X2 + X1 + X2 + 1 Boolean PolynomialRing in X0, X1, X2 (X0, X1, X2) 3 [X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1] 8 [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 0 0 0 0 0 1] Defining x0, x1, x2 [x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1] (0, 1, 0, 0, 1, 0, 0, 0) (0, 1, 0, 0, 1, 0, 0, 0) (-2, 1, 0, 0, 1, 0, 0, 0) 'x' 'XcGGt1u' [1 1 1 1 1 1 1 1] [0 1 0 1 0 1 0 1] [0 0 1 1 0 0 1 1] [0 0 0 1 0 0 0 1] [0 0 0 0 1 1 1 1] [0 0 0 0 0 1 0 1] [0 0 0 0 0 0 1 1] [0 0 0 0 0 0 0 1] X0*X1*X2 (X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1) (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1) (X0*X1*X2, X0*X1*X2 + X0*X1, X0*X1*X2 + X0*X2, X0*X1*X2 + X0*X1 + X0*X2 + X0, X0*X1*X2 + X1*X2, X0*X1*X2 + X0*X1 + X1*X2 + X1, X0*X1*X2 + X0*X2 + X1*X2 + X2, X0*X1*X2 + X0*X1 + X0*X2 + X0 + X1*X2 + X1 + X2 + 1) (X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1) 
v=ProbSceme(3)
v
a=0
for ii in v:
    a=a+ii
a
#mg=MGinv%2
MG
v*MG
Mvect
(MGinv*Mvect)%2
Mvect = vector([0,0,0,0,1,1,0,0])
Mvect
(MGinv*Mvect)%2
(MG*Mvect)%2
ProbM(X1+X2*X0)
A=ProbM2(X1+X2*X0)
A
(0.216897236248995, 0.125181737169618, 0.123043301373761, 0.0748467738425043, 0.125643838286659, 0.0283746925135870, 0.122649911882767, 0.183362508682108) 1.00000000000000 [1 1 1 1 1 1 1 1] [0 1 0 1 0 1 0 1] [0 0 1 1 0 0 1 1] [0 0 0 1 0 0 0 1] [0 0 0 0 1 1 1 1] [0 0 0 0 0 1 0 1] [0 0 0 0 0 0 1 1] [0 0 0 0 0 0 0 1] (0.216897236248995, 0.342078973418613, 0.339940537622756, 0.539969048634878, 0.342541074535654, 0.496097504218859, 0.588234287792182, 1.00000000000000) (0, 1, 0, 0, 1, 0, 0, 0) (0, 1, 0, 0, 1, 0, 0, 0) (0, 0, 0, 0, 1, 1, 0, 0) (0, 1, 0, 0, 0, 1, 0, 0) (0, 1, 0, 0, 0, 1, 0, 0) Defining x0, x1, x2 -2*x0*x1*x2 + x0*x2 + x1 Defining x0, x1, x2 ((-2, 0, 1, 0, 0, 1, 0, 0), (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1)) 
myBP=X0*X2+ X1*X0+X0
myBP
BPvars = myBP.variables()
BPmons = myBP.monomials()
BPvars
BPmons
nn = 3
atoms = 2**nn # number of atoms in corresponding boolean algebra
    # should abort this with error, if too many variables in polynomial
Mvect = [0]*atoms
Mvect
len(BPmons)
for jj in range(len(BPmons)):
    mm = 0
    BPmons[jj].variables()
    for kk in range(len(BPmons[jj].variables())):
        #BPvars.index(BPmons[jj].variables()[kk])
        mm = mm + 2**BPvars.index(BPmons[jj].variables()[kk])
    'mm='
    mm
    Mvect[mm] = 1
Mvect
Mvect.reverse()
Mvect
Mvect = vector(Mvect)
MGinv = matrix(ZZ,[1])
for ii in range(1,nn+1):
    MGinv = block_matrix(ZZ,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
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
MGinv
Avect = (MGinv*Mvect)%2
Pvect = MGinv*Avect
Avect
Pvect
ProbM(myBP)
ProbM2(myBP)
Prob(myBP)
X0*X1 + X0*X2 + X0 (X0, X1, X2) [X0*X1, X0*X2, X0] [0, 0, 0, 0, 0, 0, 0, 0] 3 (X0, X1) 'mm=' 3 (X0, X2) 'mm=' 5 (X0,) 'mm=' 1 [0, 1, 0, 1, 0, 1, 0, 0] [0, 0, 1, 0, 1, 0, 1, 0] [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 0 0 0 0 0 1] (1, 0, 0, 0, 0, 0, 1, 0) (2, 0, -1, 0, -1, 0, 1, 0) Defining x0, x1, x2 2*x0*x1*x2 - x0*x1 - x0*x2 + x0 Defining x0, x1, x2 ((2, 0, -1, 0, -1, 0, 1, 0), (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1)) Defining x0, x1, x2 2*x0*x1*x2 - x0*x2 - x0*x1 + x0 
myBP=X0*X1+X0*X2+X0
myBP
if type(myBP) == sage.rings.integer.Integer:
    if myBP%2 == 0:
        myQP = QQ(0)
    else:
        myQP = QQ(1)
    #myBP.variables()
mystr = str(myBP.variables()) # converting all uppers to lower and lowers to uppercase
    #print(mystr)
NewGens=""
for ii in range(len(mystr)):
    if mystr[ii] == mystr[ii].upper():
        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(',','')
    #NewGens.replace(', ',',')
NewGens.replace(' ','')
    #print(mystr)
    #print(NewGens)
CPR = PolynomialRing(QQ,len(myBP.variables()),NewGens,order='invlex')
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(myBP.monomials())):
    mm = 0
    for kk in range(len(myBP.monomials()[jj].variables())):
        mm = mm + 2**myBP.variables().index(myBP.monomials()[jj].variables()[kk])
    Mvect[mm] = 1
Mvect.reverse()
mypoly=myBP.parent(1)
for ii in range(0,len(CPR.gens())):
    mypoly = mypoly*(myBP.parent(1) + myBP.parent().gen(ii))
#allmons=((1+X0)*(1+X1)*(1+X2)).monomials()
allmons = mypoly.monomials()
'Mvect'
Mvect
    #print Mvect
    # 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
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
Avect = (MGinv*Mvect)%2
'Avect'
Avect
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
myQP = Pvect*vector(newpoly.monomials()) # inner product of the two vectors to create prob poly
'Pvect'
Pvect
vector(newpoly.monomials())
myQP
Mvect2 = [0]*atoms
for y in myBP.monomials():
    for jj in range(len(allmons)):
        if y==allmons[jj]:
            Mvect2[jj]=1
Mvect2
myBP.parent().term_order()
X0*X1 + X0*X2 + X0 'x0,x1,x2' Defining x0, x1, x2 'Mvect' [0, 0, 1, 0, 1, 0, 1, 0] 'Avect' (1, 0, 0, 0, 0, 0, 1, 0) 'Pvect' (2, 0, -1, 0, -1, 0, 1, 0) (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1) 2*x0*x1*x2 - x0*x2 - x0*x1 + x0 [0, 1, 1, 1, 0, 0, 0, 0] lp term order 
'0'*(7-len(bin(0)[2:]))+bin(0)[2:]
bin(1)[2:]
bin(4)[2:]
n=0
n.nbits()
''
#n.binary()
#'0' * (3-n.nbits()) + n.binary()
nn=4
n=0
n.binary()
for ii in range(2**nn):
    if n==0:
        '0'*nn
    else:
        '0'*(nn-n.nbits())+n.binary()
    n=n+1
'0000000' '1' '100' 0 '' '0' '0000' '0001' '0010' '0011' '0100' '0101' '0110' '0111' '1000' '1001' '1010' '1011' '1100' '1101' '1110' '1111' 
from tabulate import tabulate
load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
load('ProbM2.py')
load('ProbScheme.py')
load('IndProbSchema.py')
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])
nn = 3 # of generators of Boolean Algebra
print('Creating a Boolean polynomial ring with '+str(nn)+ ' generators which will act as two state random variables.')
print(' ')
aa = 2**nn # number of atoms in the Boolean Algebra
BPR = BooleanPolynomialRing(nn,'X',order='lp')
print('The generators of the Boolean Polynomial Ring, BPR')
BPR.inject_variables()
BPR.term_order()
mypoly=BPR(1)
for ii in range(0,len(BPR.gens())):
    mypoly=mypoly*(BPR(1) + BPR.gen(ii))
print(' ')
CPR = PolynomialRing(QQ,nn,'x',order='lp') # probability polynomials over QQ
print('Generators for the probability polynomials over QQ')
CPR.inject_variables()
CPR.term_order()
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(' ')
IPS=IndProbSchema(nn)
print("Probabilities (heads) of the "+str(nn)+" generators taken as distinguishable (unfair) coin tosses. These are independent probabilities, and do not sum to 1.0")
prbheads = [round(x,6) for x in IPS[0]]
hdata = []
headsdata = ['prob. of heads']
headers=['generators']
cola=['left']
trials = 30
cointoss = matrix(ZZ,30,nn,0)
frequency = vector([0.0]*nn)
for ii in range(nn):
    headsdata.append(str(round(prbheads[ii],6)))
    headers.append('P(' + str(BPR.gens()[ii]) + ') = ' + str(CPR.gens()[ii]))
    cola.append('center')
hdata.append(headsdata)
for ii in range(trials):
    trialrow = ['coin toss '+ str(ii+1)]
    for jj in range(nn):
        if random() <= prbheads[jj]:
            cointoss[ii,jj] = 1
            trialrow.append('1')
        else:
            trialrow.append('0')
        frequency[jj] = frequency[jj] + cointoss[ii,jj]
    hdata.append(trialrow)
frequency = frequency/trials
trialrow = ['freq. of heads']
for jj in range(nn):
    trialrow.append(round(frequency[jj],6))
hdata.append(trialrow)
print(tabulate(hdata, headers=headers, tablefmt="grid", colalign=cola))
print(' ')
print('The number of atoms is '+str(aa)+' = 2**'+str(nn)+' generators. The generator probabilities induce the following probability scheme for the '+str(aa)+' atoms.')
print('Aprob')
[round(x,6) for x in IPS[1]]
print(' ')
mysum=0
for ii in range(aa):
    mysum=mysum+IPS[1][ii]
print("Total probability (just checking)")
round(mysum,6)
print(' ')
print('Table of atoms and monomials. With the exception of the 0th atom and monomial the rest are not equal.')
AMdata = []
for ii in range(aa):
    rowAM = []
    mystring = '0'*(nn-len(bin(ii)[2:]))+bin(ii)[2:]
    mystring2 = ''
    for jj in range(nn):
        mystring2 = mystring2 + '*(' + str(BPR.gen(jj)) + '+' + mystring[jj] +')'
    rowAM.append(str(ii))
    rowAM.append(mystring)
    rowAM.append(mystring2[1:])
    rowAM.append(str(mypoly.monomials()[ii]))
    AMdata.append(rowAM)
    #print(str(ii)+'  '+ mystring +'  ' + mystring2[1:] + '  ' + str(mypoly.monomials()[ii]))
print(tabulate(AMdata, headers=('#','bin','atoms', 'monomials'), tablefmt="grid", colalign=('center','center','center','left')))
print(' ')
###############################################################################################
# should abort this with error, if too many variables in polynomial
myBP = X0*X1 + X0*X2 + X0 + 1
###############################################################################################
print('An example Boolean polynomial:')
print(' ')
print('myBP')
myBP
print(' ')
Mvect = [0]*aa
for y in myBP.monomials():
    for jj in range(len(mypoly.monomials())):
        if y==mypoly.monomials()[jj]:
            Mvect[jj]=1
MGAinv = matrix(ZZ,[1])
for ii in range(1,nn+1):
    MGAinv = block_matrix(ZZ,2,2,[MGAinv,-MGAinv,0,MGAinv],subdivide=False) # creating inverse of multigrade operator 'AND' matrix MGA
        # modulo 2 it's an involution however...
        # MGA and MGAinv are upper triangular, MGA is like Sierpinski triangle see
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
        # over ZZ it is totally unimodular
MGA = MGAinv%2
Mvect=vector(Mvect)
print('Monomial vector')
print('MonVec')
print(str(vector(mypoly.monomials())))
print(' ')
print('Polynomial vector')
print('Mvect = (MGA*Avect)mod 2, Mvect is 0-1 vector with a 1 corresponding to each monomial in the Bollean polynomial, 0 otherwise')
Mvect
print(' ')
Avect = MGA*Mvect%2
print('Atom vector')
print('Avect = (MGA*Mvect)mod 2')
Avect
print(' ')
print("Multigrade And matrix, it's totally unimodular over ZZ, and it's an involution(self inverse) mod 2")
print('MGA =')
MGA
print(' ')
print('Inverse of Multigrade And matrix in ZZ, QQ, or RR')
print('MGAinv')
MGAinv
Pvect = MGAinv*Avect
print(' ')
Mprob=IPS[1]*MGA
print('Monomial probabilities(which do not sum to 1.0)')
print('but the last term is the sum of all atom probabilities = 1.0')
print('Mprob = Aprob*MGA')
[round(x,6) for x in Mprob]
print(' ')
print('Generators for the probability polynomial or formula over RR')
print(' ')
print('Pvect = MGAinv*Avect, the coefficients of the probability polynomial or formula')
[round(x,6) for x in Pvect]
print(' ')
print('Monoms, monomials over RR')
Monoms = vector(newpoly.monomials())
print(' ')
print('Probability polynomial when Boolean Ring generators are independent two state random variables')
print('ProbPoly = Pvect*Monoms')
Pvect*Monoms
print(' ')
ProbExpr = ''
for ii in range(aa):
    if Pvect[ii] > 1: mystr = ' + '+str(Pvect[ii])+'*'
    if Pvect[ii] == 1: mystr = str(' + ')
    if Pvect[ii] == -1: mystr = str(' - ')
    if Pvect[ii] < -1: mystr = ' - '+str(abs(Pvect[ii]))+'*'
    if Pvect[ii] != 0: ProbExpr = ProbExpr + mystr+'P('+ str(mypoly.monomials()[ii])+')'
if ProbExpr[1] == '+' : ProbExpr = ProbExpr[3:]
if ProbExpr[1] == '-' : ProbExpr = '-' + ProbExpr[3:]
print('Probability formula whether or not the generators are independent,')
print('P('+str(myBP)+') = ' + ProbExpr)
print(' ')
print("In this example with independent generators the numerical probability (using the dot product of vectors), P("+str(myBP)+") = ")
print('Pvect*Mprob =')
Pvect*Mprob
print('or Avect*Aprob =')
IPS[1]*Avect
print(' ')
print('An example where the generators are not independent is constructed by creating a probability schema for the atoms arbitrarily')
PS = ProbScheme(nn)
print('Probability schema for the atoms, PS =')
[round(x,6) for x in PS]
print(' ')
MP=PS*MGA
print('Monomial probabilities')
print('MP = PS*MGA =')
[round(x,6) for x in MP]
print(' ')
print('In this example with non-independent generators the numerical probability, P('+str(myBP)+') = ')
print('Pvect*MP =')
Pvect*MP
print('Avect*PS =')
Avect*PS
Creating a Boolean polynomial ring with 3 generators which will act as two state random variables. The generators of the Boolean Polynomial Ring, BPR Defining X0, X1, X2 lp term order Generators for the probability polynomials over QQ Defining x0, x1, x2 lp term order Probabilities (heads) of the 3 generators taken as distinguishable (unfair) coin tosses. These are independent probabilities, and do not sum to 1.0 +----------------+--------------+--------------+--------------+ | generators | P(X0) = x0 | P(X1) = x1 | P(X2) = x2 | +================+==============+==============+==============+ | prob. of heads | 0.970432 | 0.994518 | 0.119172 | +----------------+--------------+--------------+--------------+ | coin toss 1 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 2 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 3 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 4 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 5 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 6 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 7 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 8 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 9 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 10 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 11 | 1 | 1 | 1 | +----------------+--------------+--------------+--------------+ | coin toss 12 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 13 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 14 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 15 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 16 | 1 | 1 | 1 | +----------------+--------------+--------------+--------------+ | coin toss 17 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 18 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 19 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 20 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 21 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 22 | 1 | 1 | 1 | +----------------+--------------+--------------+--------------+ | coin toss 23 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 24 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 25 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 26 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 27 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 28 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 29 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 30 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | freq. of heads | 1 | 1 | 0.1 | +----------------+--------------+--------------+--------------+ The number of atoms is 8 = 2**3 generators. The generator probabilities induce the following probability scheme for the 8 atoms. Aprob [0.115014, 0.850098, 0.000634, 0.004686, 0.003504, 0.025901, 1.9e-05, 0.000143] Total probability (just checking) 1.0 Table of atoms and monomials. With the exception of the 0th atom and monomial the rest are not equal. +-----+-------+----------------------+-------------+ | # | bin | atoms | monomials | +=====+=======+======================+=============+ | 0 | 000 | (X0+0)*(X1+0)*(X2+0) | X0*X1*X2 | +-----+-------+----------------------+-------------+ | 1 | 001 | (X0+0)*(X1+0)*(X2+1) | X0*X1 | +-----+-------+----------------------+-------------+ | 2 | 010 | (X0+0)*(X1+1)*(X2+0) | X0*X2 | +-----+-------+----------------------+-------------+ | 3 | 011 | (X0+0)*(X1+1)*(X2+1) | X0 | +-----+-------+----------------------+-------------+ | 4 | 100 | (X0+1)*(X1+0)*(X2+0) | X1*X2 | +-----+-------+----------------------+-------------+ | 5 | 101 | (X0+1)*(X1+0)*(X2+1) | X1 | +-----+-------+----------------------+-------------+ | 6 | 110 | (X0+1)*(X1+1)*(X2+0) | X2 | +-----+-------+----------------------+-------------+ | 7 | 111 | (X0+1)*(X1+1)*(X2+1) | 1 | +-----+-------+----------------------+-------------+ An example Boolean polynomial: myBP X0*X1 + X0*X2 + X0 + 1 Monomial vector MonVec (X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1) Polynomial vector Mvect = (MGA*Avect)mod 2, Mvect is 0-1 vector with a 1 corresponding to each monomial in the Bollean polynomial, 0 otherwise (0, 1, 1, 1, 0, 0, 0, 1) Atom vector Avect = (MGA*Mvect)mod 2 (0, 1, 1, 0, 1, 1, 1, 1) Multigrade And matrix, it's totally unimodular over ZZ, and it's an involution(self inverse) mod 2 MGA = [1 1 1 1 1 1 1 1] [0 1 0 1 0 1 0 1] [0 0 1 1 0 0 1 1] [0 0 0 1 0 0 0 1] [0 0 0 0 1 1 1 1] [0 0 0 0 0 1 0 1] [0 0 0 0 0 0 1 1] [0 0 0 0 0 0 0 1] Inverse of Multigrade And matrix in ZZ, QQ, or RR MGAinv [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 0 0 0 0 0 1] Monomial probabilities(which do not sum to 1.0) but the last term is the sum of all atom probabilities = 1.0 Mprob = Aprob*MGA [0.115014, 0.965112, 0.115648, 0.970432, 0.118519, 0.994518, 0.119172, 1.0] Generators for the probability polynomial or formula over RR Pvect = MGAinv*Avect, the coefficients of the probability polynomial or formula [-2.0, 1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 1.0] Monoms, monomials over RR Probability polynomial when Boolean Ring generators are independent two state random variables ProbPoly = Pvect*Monoms -2*x0*x1*x2 + x0*x1 + x0*x2 - x0 + 1 Probability formula whether or not the generators are independent, P(X0*X1 + X0*X2 + X0 + 1) = -2*P(X0*X1*X2) + P(X0*X1) + P(X0*X2) - P(X0) + P(1) In this example with independent generators the numerical probability (using the dot product of vectors), P(X0*X1 + X0*X2 + X0 + 1) = Pvect*Mprob = 0.880299732641442 or Avect*Aprob = 0.880299732641442 An example where the generators are not independent is constructed by creating a probability schema for the atoms arbitrarily Probability schema for the atoms, PS = [0.036326, 0.026352, 0.219059, 0.037976, 0.113768, 0.101268, 0.225121, 0.240132] Monomial probabilities MP = PS*MGA = [0.036326, 0.062677, 0.255384, 0.319712, 0.150093, 0.277713, 0.594273, 1.0] In this example with non-independent generators the numerical probability, P(X0*X1 + X0*X2 + X0 + 1) = Pvect*MP = 0.925698661517172 Avect*PS = 0.925698661517172 
from tabulate import tabulate
#
nn=7
TRM=matrix(RR,nn,nn)
for ii in [1..nn]:
    for jj in [1..nn]:
        TRM[ii-1,jj-1]=int(round(abs(normalvariate(0,1000))))

NET = TRM-TRM.transpose()
for ii in [1..nn]:
    for jj in [1..nn]:
        if NET[ii-1,jj-1] < 0.0: NET[ii-1,jj-1] = 0.0

b=[0.0]*nn
x=[0.0]*(nn-1)
for kk in [1..nn]:
    for ii in [1..nn]:
        b[kk-1] = b[kk-1]-TRM[ii-1,kk-1]+TRM[kk-1,ii-1]

m = 2*nn
for ii in [1..nn-1]:
    for jj in [1..nn-1]:
        x[ii-1] = b[jj-1] + x[ii-1]
        if ii == jj: x[ii-1] = b[jj-1] + x[ii-1]
    x[ii-1] = x[ii-1]/m

LMD = matrix(RR,nn,nn)
for ii in [2..nn]:
    for jj in [1..ii-1]:
        for kk in [1..nn-1]:
            if ii==kk: LMD[ii-1,jj-1] = LMD[ii-1,jj-1] + x[kk-1]
            if jj==kk: LMD[ii-1,jj-1] = LMD[ii-1,jj-1] - x[kk-1]
        LMD[jj-1,ii-1] = -LMD[ii-1,jj-1]

print(' ')
NETL = LMD-LMD.transpose()
for ii in range(nn):
    for jj in range(nn):
        if NETL[ii,jj] < 0.0: NETL[ii,jj] = 0.0
CYC = TRM-LMD
CYCS = (CYC+CYC.transpose())/2 # symmetric cycle
CYCK = (CYC-CYC.transpose())/2 # skew-symmetric cycle
#rep=matrix(ZZ,nn,1,[1]*nn)
recvTRM = matrix(RR,nn,1,0)
sendTRM = matrix(RR,nn,1,0)
recvNET = matrix(RR,nn,1,0)
sendNET = matrix(RR,nn,1,0)
recvLMD = matrix(RR,nn,1,0)
sendLMD = matrix(RR,nn,1,0)
recvNETL = matrix(RR,nn,1,0)
sendNETL = matrix(RR,nn,1,0)
recvCYCS = matrix(RR,nn,1,0)
sendCYCS = matrix(RR,nn,1,0)
recvCYCK = matrix(RR,nn,1,0)
sendCYCK = matrix(RR,nn,1,0)
for ii in range(nn):
    recvTRM[ii,0] = sum(TRM[ii])
    recvNET[ii,0] = sum(NET[ii])
    recvLMD[ii,0] = sum(LMD[ii])
    recvNETL[ii,0] = sum(NETL[ii])
    recvCYCS[ii,0] = sum(CYCS[ii])
    recvCYCK[ii,0] = sum(CYCK[ii])
    sendTRM[ii,0] = sum(TRM.transpose()[ii])
    sendNET[ii,0] = sum(NET.transpose()[ii])
    sendLMD[ii,0] = sum(LMD.transpose()[ii])
    sendNETL[ii,0] = sum(NETL.transpose()[ii])
    sendCYCS[ii,0] = sum(CYCS.transpose()[ii])
    sendCYCK[ii,0] = sum(CYCK.transpose()[ii])

Nodes = [' ']
for ii in range(nn):
    Nodes.append('N'+str(ii))
Nodes.append(' ')
Nodes.append('R')
Nodes.append(' ')
Nodes.append('R-S')
TRMdata=[]
NETdata=[]
LMDdata=[]
NETLdata=[]
CYCSdata = []
CYCKdata = []
for ii in range(nn):
    rowT = [Nodes[ii+1]]
    rowN = [Nodes[ii+1]]
    rowL = [Nodes[ii+1]]
    rowNL = [Nodes[ii+1]]
    rowCS = [Nodes[ii+1]]
    rowCK = [Nodes[ii+1]]
    for jj in range(nn):
        rowT.append('{0:.2f}'.format(TRM[ii,jj]))
        rowN.append('{0:.2f}'.format(NET[ii,jj]))
        rowL.append('{0:.2f}'.format(LMD[ii,jj]))
        rowNL.append('{0:.2f}'.format(NETL[ii,jj]))
        rowCS.append('{0:.2f}'.format(CYCS[ii,jj]))
        rowCK.append('{0:.2f}'.format(CYCK[ii,jj]))
    rowT.append(' ')
    rowT.append('{0:.2f}'.format(recvTRM[ii,0]))
    rowT.append(' ')
    rowT.append('{0:.2f}'.format((recvTRM-sendTRM)[ii,0]))

    rowN.append(' ')
    rowN.append('{0:.2f}'.format(recvNET[ii,0]))
    rowN.append(' ')
    rowN.append('{0:.2f}'.format((recvNET-sendNET)[ii,0]))

    rowL.append(' ')
    rowL.append('{0:.2f}'.format(recvLMD[ii,0]))
    rowL.append(' ')
    rowL.append('{0:.2f}'.format((recvLMD-sendLMD)[ii,0]))

    rowNL.append(' ')
    rowNL.append('{0:.2f}'.format(recvNETL[ii,0]))
    rowNL.append(' ')
    rowNL.append('{0:.2f}'.format((recvNETL-sendNETL)[ii,0]))

    rowCS.append(' ')
    rowCS.append('{0:.2f}'.format(recvCYCS[ii,0]))
    rowCS.append(' ')
    rowCS.append('{0:.2f}'.format((recvCYCS-sendCYCS)[ii,0]))

    rowCK.append(' ')
    rowCK.append('{0:.2f}'.format(recvCYCK[ii,0]))
    rowCK.append(' ')
    rowCK.append('{0:.2f}'.format((recvCYCK-sendCYCK)[ii,0]))

    TRMdata.append(rowT)
    NETdata.append(rowN)
    LMDdata.append(rowL)
    NETLdata.append(rowNL)
    CYCSdata.append(rowCS)
    CYCKdata.append(rowCK)

TRMdata.append([' ']*(nn+5))
TRMdata.append([' ']*(nn+5))
TRMdata[nn+1][0] = 'S'
TRMdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvTRM)))
TRMdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvTRM-sendTRM)))
NETdata.append([' ']*(nn+5))
NETdata.append([' ']*(nn+5))
NETdata[nn+1][0] = 'S'
NETdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvNET)))
NETdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvNET-sendNET)))
LMDdata.append([' ']*(nn+5))
LMDdata.append([' ']*(nn+5))
LMDdata[nn+1][0] = 'S'
LMDdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvLMD)))
LMDdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvLMD-sendLMD)))
NETLdata.append([' ']*(nn+5))
NETLdata.append([' ']*(nn+5))
NETLdata[nn+1][0] = 'S'
NETLdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvNETL)))
NETLdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvNETL-sendNETL)))
CYCSdata.append([' ']*(nn+5))
CYCSdata.append([' ']*(nn+5))
CYCSdata[nn+1][0] = 'S'
CYCSdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvCYCS)))
CYCSdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvCYCS-sendCYCS)))
CYCKdata.append([' ']*(nn+5))
CYCKdata.append([' ']*(nn+5))
CYCKdata[nn+1][0] = 'S'
CYCKdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvCYCK)))
CYCKdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvCYCK-sendCYCK)))
SKEW = TRM-TRM.transpose()
SKEWdata = []
for ii in range(nn):
    TRMdata[nn+1][ii+1]='{0:.2f}'.format(sendTRM[ii,0])
    NETdata[nn+1][ii+1]='{0:.2f}'.format(sendNET[ii,0])
    LMDdata[nn+1][ii+1]='{0:.2f}'.format(sendLMD[ii,0])
    NETLdata[nn+1][ii+1]='{0:.2f}'.format(sendNETL[ii,0])
    CYCSdata[nn+1][ii+1]='{0:.2f}'.format(sendCYCS[ii,0])
    CYCKdata[nn+1][ii+1]='{0:.2f}'.format(sendCYCK[ii,0])
    rowSK = [Nodes[ii+1]]
    for jj in range(nn):
        rowSK.append('{0:.2f}'.format(SKEW[ii][jj]))
    SKEWdata.append(rowSK)

cola = ["right" for x in Nodes]
cola[0] = "center"
Nodes[0]='TRM'
mystring = 'The square matrix, TRM, of data in the following table represents the flow of goods and services in a complete directed multigraph of '+str(nn)+' nodes as valued in a single currency. The nodes of the graph might represent the nations of a small world. The data can be interpreted in either of two ways. Here we choose the interpretation that the value, X, in the cell corresponding to row 0 and column 1 is the amount received by node 0 from node 1. The value, Y, corresponding to row 1 and column 0 is the amount received by node 1 from node 0. One of these values will generally be greater than(though possibly equal to) the other. The net effect is a positive flow of |X-Y| in one direction or the other. A negative entry represents a flow in the direction opposite to a given edge in the graph. The column, R, is the sum along each row, and represents the total amount received from the network for each node. The row, S, represents the amount sent into the network for each node. The column, R-S, represents the net inflow or outflow of money for each node during a given time step. The data for this table was generated with a random variable from a half-normal distribution, i.e. |Z| where Z ~ N(0,sigma^2). TRM is my abbreviation for transaction matrix.'
print(mystring)
print(tabulate(TRMdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
cola2 = ["right" for ii in range(nn)]
cola2[0] = 'center'
Nodes[0] = 'SKEW'
mystring = "To see the net effect of money transfers from node to node, we first form a skew-symmetric matrix(each entry is the negative of it's mirror reflection through the main diagonal) from TRM given by SKEW = TRM - TRM.transpose()."
print(mystring)
print(tabulate(SKEWdata, headers=Nodes, tablefmt="grid", colalign=cola2))
print(' ')
Nodes[0]='NET'
mystring = "Next, we set the negative entries of SKEW to 0, and keep the rest to obtain the matrix, NET. For instance if you send me $10, and I send you $5, then the net effect is you sent me $5, and I sent you $0, we're just expressing everything in non-negative numbers. We could apply x*H(x) where H is the Heaviside step function to each entry of SKEW to obtain NET. Note that the R-S column vector for NET has not changed from that of TRM. We also note that for any square matrix, the sum over the R vector equals the sum over the S vector which is just the sum of all entries in the matrix, and therefore the R-S vector always sums to 0."
print(mystring)
print(tabulate(NETdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='LMD'
mystring = 'Closer analysis using some graph theory and linear algebra allows us to express the net effect of all transactions in TRM efficiently as the projection of TRM into the cocycle space of the complete directed multigraph on N vertices and 2*N*(N-1) edges. The dimension of the cocycle space is N-1. Projecting TRM thusly we obtain the skew-symmetric matrix, LMD. I called it this for no good reason having inherited it from a Visual Basic version I wrote a number of years ago. Note here also, the R-S vector remains unchanged from TRM. As a further note the numbers here are rounded to two decimal places whereas they have many more significant digits so the numbers may not add up exactly, if checking.'
print(mystring)
print(tabulate(LMDdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='NETL'
mystring = 'To express the transactions in LMD as non-negative numbers only form SKEW = LMD - LMD.transpose, and set the negatives to 0. We obtain the matrix NETL. Again R-S is unchanged.'
print(mystring)
print(tabulate(NETLdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='CYCS'
mystring = "Any matrix where the R-S vector is the 0-vector(not merely summing to 0, but every entry is 0) lies in the cycle space of the complete directed graph. We can project TRM into the cycle space to obtain a matrix, CYC. This can be obtained by simply subtracting LMD from TRM, thus CYC = TRM - LMD. CYC can be broken down further into symmetric and skew-symmetric matrices, CYCS and CYCK respectively where CYCS = (CYC + CYC.transpose())/2 and CYCK = (CYC - CYC.transpose())/2. We have TRM = LMD + CYCS + CYCK."
print(mystring)
print(tabulate(CYCSdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='CYCK'
mystring = 'Die Quantenmechanik ist sehr achtung-gebietend. Aber eine innere Stimme sagt mir, daß das doch nicht der wahre Jakob ist. Die Theorie liefert viel, aber dem Geheimnis des Alten bringt sie uns kaum näher. Jedenfalls bin ich überzeugt, daß der nicht würfelt.'
mystring = "The space of all N×N square matrices' dimension is N^2. The cocycle space is dimension N-1, so that the cycle space is dimension N^2 - N + 1. The cycle space is further broken down into the symmetric cycles and the skew-symmetric cycles. Every symmetric matrix is a cycle, and the dimension of the symmetric matrices is N*(N+1)/2. The dimension of the skew-symmetric matrices is N*(N-1)/2 leaving the dimension of the skew-symmetric cycles (N-1)*(N-2)/2."
print(mystring)
print(tabulate(CYCKdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
The square matrix, TRM, of data in the following table represents the flow of goods and services in a complete directed multigraph of 7 nodes as valued in a single currency. The nodes of the graph might represent the nations of a small world. The data can be interpreted in either of two ways. Here we choose the interpretation that the value, X, in the cell corresponding to row 0 and column 1 is the amount received by node 0 from node 1. The value, Y, corresponding to row 1 and column 0 is the amount received by node 1 from node 0. One of these values will generally be greater than(though possibly equal to) the other. The net effect is a positive flow of |X-Y| in one direction or the other. A negative entry represents a flow in the direction opposite to a given edge in the graph. The column, R, is the sum along each row, and represents the total amount received from the network for each node. The row, S, represents the amount sent into the network for each node. The column, R-S, represents the net inflow or outflow of money for each node during a given time step. The data for this table was generated with a random variable from a half-normal distribution, i.e. |Z| where Z ~ N(0,sigma^2). TRM is my abbreviation for transaction matrix. +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | TRM | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +=======+=========+==========+=========+=========+=========+=========+=========+=====+==========+=====+==========+ | N0 | 1606.00 | 877.00 | 343.00 | 764.00 | 731.00 | 495.00 | 1351.00 | | 6167.00 | | -859.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N1 | 459.00 | 1154.00 | 678.00 | 817.00 | 719.00 | 367.00 | 2389.00 | | 6583.00 | | -4081.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N2 | 1160.00 | 1838.00 | 539.00 | 295.00 | 672.00 | 582.00 | 193.00 | | 5279.00 | | -1926.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N3 | 1100.00 | 2073.00 | 888.00 | 1065.00 | 330.00 | 1092.00 | 542.00 | | 7090.00 | | 1607.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N4 | 588.00 | 1241.00 | 1516.00 | 1399.00 | 468.00 | 9.00 | 436.00 | | 5657.00 | | 197.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N5 | 1628.00 | 951.00 | 1535.00 | 114.00 | 1333.00 | 666.00 | 17.00 | | 6244.00 | | 2711.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N6 | 485.00 | 2530.00 | 1706.00 | 1029.00 | 1207.00 | 322.00 | 1872.00 | | 9151.00 | | 2351.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | | | | | | | | | | | | | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | S | 7026.00 | 10664.00 | 7205.00 | 5483.00 | 5460.00 | 3533.00 | 6800.00 | | 46171.00 | | 0.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ To see the net effect of money transfers from node to node, we first form a skew-symmetric matrix(each entry is the negative of it's mirror reflection through the main diagonal) from TRM given by SKEW = TRM - TRM.transpose(). +--------+------+------+-------+-------+-------+-------+-------+ | SKEW | N0 | N1 | N2 | N3 | N4 | N5 | N6 | +========+======+======+=======+=======+=======+=======+=======+ | N0 | 0 | 418 | -817 | -336 | 143 | -1133 | 866 | +--------+------+------+-------+-------+-------+-------+-------+ | N1 | -418 | 0 | -1160 | -1256 | -522 | -584 | -141 | +--------+------+------+-------+-------+-------+-------+-------+ | N2 | 817 | 1160 | 0 | -593 | -844 | -953 | -1513 | +--------+------+------+-------+-------+-------+-------+-------+ | N3 | 336 | 1256 | 593 | 0 | -1069 | 978 | -487 | +--------+------+------+-------+-------+-------+-------+-------+ | N4 | -143 | 522 | 844 | 1069 | 0 | -1324 | -771 | +--------+------+------+-------+-------+-------+-------+-------+ | N5 | 1133 | 584 | 953 | -978 | 1324 | 0 | -305 | +--------+------+------+-------+-------+-------+-------+-------+ | N6 | -866 | 141 | 1513 | 487 | 771 | 305 | 0 | +--------+------+------+-------+-------+-------+-------+-------+ Next, we set the negative entries of SKEW to 0, and keep the rest to obtain the matrix, NET. For instance if you send me $10, and I send you $5, then the net effect is you sent me $5, and I sent you $0, we're just expressing everything in non-negative numbers. We could apply x*H(x) where H is the Heaviside step function to each entry of SKEW to obtain NET. Note that the R-S column vector for NET has not changed from that of TRM. We also note that for any square matrix, the sum over the R vector equals the sum over the S vector which is just the sum of all entries in the matrix, and therefore the R-S vector always sums to 0. +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | NET | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +=======+=========+=========+=========+=========+=========+=========+========+=====+==========+=====+==========+ | N0 | 0.00 | 418.00 | 0.00 | 0.00 | 143.00 | 0.00 | 866.00 | | 1427.00 | | -859.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 0.00 | | -4081.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N2 | 817.00 | 1160.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 1977.00 | | -1926.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N3 | 336.00 | 1256.00 | 593.00 | 0.00 | 0.00 | 978.00 | 0.00 | | 3163.00 | | 1607.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N4 | 0.00 | 522.00 | 844.00 | 1069.00 | 0.00 | 0.00 | 0.00 | | 2435.00 | | 197.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N5 | 1133.00 | 584.00 | 953.00 | 0.00 | 1324.00 | 0.00 | 0.00 | | 3994.00 | | 2711.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N6 | 0.00 | 141.00 | 1513.00 | 487.00 | 771.00 | 305.00 | 0.00 | | 3217.00 | | 2351.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | | | | | | | | | | | | | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | S | 2286.00 | 4081.00 | 3903.00 | 1556.00 | 2238.00 | 1283.00 | 866.00 | | 16213.00 | | 0.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ Closer analysis using some graph theory and linear algebra allows us to express the net effect of all transactions in TRM efficiently as the projection of TRM into the cocycle space of the complete directed multigraph on N vertices and 2*N*(N-1) edges. The dimension of the cocycle space is N-1. Projecting TRM thusly we obtain the skew-symmetric matrix, LMD. I called it this for no good reason having inherited it from a Visual Basic version I wrote a number of years ago. Note here also, the R-S vector remains unchanged from TRM. As a further note the numbers here are rounded to two decimal places whereas they have many more significant digits so the numbers may not add up exactly, if checking. +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | LMD | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +=======+=========+=========+=========+=========+=========+==========+==========+=====+==========+=====+==========+ | N0 | 0.00 | 230.14 | 76.21 | -176.14 | -75.43 | -255.00 | -229.29 | | -429.50 | | -859.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N1 | -230.14 | 0.00 | -153.93 | -406.29 | -305.57 | -485.14 | -459.43 | | -2040.50 | | -4081.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N2 | -76.21 | 153.93 | 0.00 | -252.36 | -151.64 | -331.21 | -305.50 | | -963.00 | | -1926.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N3 | 176.14 | 406.29 | 252.36 | 0.00 | 100.71 | -78.86 | -53.14 | | 803.50 | | 1607.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N4 | 75.43 | 305.57 | 151.64 | -100.71 | 0.00 | -179.57 | -153.86 | | 98.50 | | 197.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N5 | 255.00 | 485.14 | 331.21 | 78.86 | 179.57 | 0.00 | 25.71 | | 1355.50 | | 2711.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N6 | 229.29 | 459.43 | 305.50 | 53.14 | 153.86 | -25.71 | 0.00 | | 1175.50 | | 2351.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | | | | | | | | | | | | | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | S | 429.50 | 2040.50 | 963.00 | -803.50 | -98.50 | -1355.50 | -1175.50 | | 0.00 | | 0.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ To express the transactions in LMD as non-negative numbers only form SKEW = LMD - LMD.transpose, and set the negatives to 0. We obtain the matrix NETL. Again R-S is unchanged. +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | NETL | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +========+=========+=========+=========+========+========+======+=======+=====+=========+=====+==========+ | N0 | 0.00 | 460.29 | 152.43 | 0.00 | 0.00 | 0.00 | 0.00 | | 612.71 | | -859.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 0.00 | | -4081.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N2 | 0.00 | 307.86 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 307.86 | | -1926.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N3 | 352.29 | 812.57 | 504.71 | 0.00 | 201.43 | 0.00 | 0.00 | | 1871.00 | | 1607.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N4 | 150.86 | 611.14 | 303.29 | 0.00 | 0.00 | 0.00 | 0.00 | | 1065.29 | | 197.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N5 | 510.00 | 970.29 | 662.43 | 157.71 | 359.14 | 0.00 | 51.43 | | 2711.00 | | 2711.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N6 | 458.57 | 918.86 | 611.00 | 106.29 | 307.71 | 0.00 | 0.00 | | 2402.43 | | 2351.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | | | | | | | | | | | | | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | S | 1471.71 | 4081.00 | 2233.86 | 264.00 | 868.29 | 0.00 | 51.43 | | 8970.29 | | 0.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ Any matrix where the R-S vector is the 0-vector(not merely summing to 0, but every entry is 0) lies in the cycle space of the complete directed graph. We can project TRM into the cycle space to obtain a matrix, CYC. This can be obtained by simply subtracting LMD from TRM, thus CYC = TRM - LMD. CYC can be broken down further into symmetric and skew-symmetric matrices, CYCS and CYCK respectively where CYCS = (CYC + CYC.transpose())/2 and CYCK = (CYC - CYC.transpose())/2. We have TRM = LMD + CYCS + CYCK. +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | CYCS | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +========+=========+=========+=========+=========+=========+=========+=========+=====+==========+=====+=======+ | N0 | 1606.00 | 668.00 | 751.50 | 932.00 | 659.50 | 1061.50 | 918.00 | | 6596.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N1 | 668.00 | 1154.00 | 1258.00 | 1445.00 | 980.00 | 659.00 | 2459.50 | | 8623.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N2 | 751.50 | 1258.00 | 539.00 | 591.50 | 1094.00 | 1058.50 | 949.50 | | 6242.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N3 | 932.00 | 1445.00 | 591.50 | 1065.00 | 864.50 | 603.00 | 785.50 | | 6286.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N4 | 659.50 | 980.00 | 1094.00 | 864.50 | 468.00 | 671.00 | 821.50 | | 5558.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N5 | 1061.50 | 659.00 | 1058.50 | 603.00 | 671.00 | 666.00 | 169.50 | | 4888.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N6 | 918.00 | 2459.50 | 949.50 | 785.50 | 821.50 | 169.50 | 1872.00 | | 7975.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | | | | | | | | | | | | | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | S | 6596.50 | 8623.50 | 6242.00 | 6286.50 | 5558.50 | 4888.50 | 7975.50 | | 46171.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ The space of all N×N square matrices' dimension is N^2. The cocycle space is dimension N-1, so that the cycle space is dimension N^2 - N + 1. The cycle space is further broken down into the symmetric cycles and the skew-symmetric cycles. Every symmetric matrix is a cycle, and the dimension of the symmetric matrices is N*(N+1)/2. The dimension of the skew-symmetric matrices is N*(N-1)/2 leaving the dimension of the skew-symmetric cycles (N-1)*(N-2)/2. +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | CYCK | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +========+=========+=========+=========+=========+=========+=========+=========+=====+=======+=====+=======+ | N0 | 0.00 | -21.14 | -484.71 | 8.14 | 146.93 | -311.50 | 662.29 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N1 | 21.14 | 0.00 | -426.07 | -221.71 | 44.57 | 193.14 | 388.93 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N2 | 484.71 | 426.07 | 0.00 | -44.14 | -270.36 | -145.29 | -451.00 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N3 | -8.14 | 221.71 | 44.14 | 0.00 | -635.21 | 567.86 | -190.36 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N4 | -146.93 | -44.57 | 270.36 | 635.21 | 0.00 | -482.43 | -231.64 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N5 | 311.50 | -193.14 | 145.29 | -567.86 | 482.43 | 0.00 | -178.21 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N6 | -662.29 | -388.93 | 451.00 | 190.36 | 231.64 | 178.21 | 0.00 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | | | | | | | | | | | | | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | S | 0.00 | 0.00 | 0.00 | 0.00 | -0.00 | 0.00 | 0.00 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ 
CPR = PolynomialRing(QQ,nn,'x',order=myBP.parent().term_order())
CPR.inject_variables()
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
newpoly
Defining x0, x1, x2 x0*x1*x2 + x0*x1 + x0*x2 + x0 + x1*x2 + x1 + x2 + 1 
from bitarray import bitarray
b = bitarray()
b.append(1)
b
b.extend([1,0])
b
len(b)
rows = 10
cols = 22
B=[]
B.append(rows)
B.append(cols)
B.append(bitarray(rows*cols))
B.append(bitarray(rows*cols))
B
B[2][B[1]*5+19]
B[3][B[1]*5+19]=1
B[3][B[1]*5+19]
B[3]
bitarray('1') bitarray('110') 3 [10, 22, bitarray('0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000'), bitarray('0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000')] 0 1 bitarray('0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000') 
import csv
import sys
from bitarray import bitarray
from tabulate import tabulate
load('RVal.py')
load('SVal.py')
load('Xchange.py')
load('XchangeA.py')
load('MinShip.py')
load('printPHI.py')
rows = 10
cols = 22
PHI=[]
PHI.append(rows)
PHI.append(cols)
PHI.append(bitarray(rows*cols))
PHI.append(bitarray(rows*cols))
f = open('Book1.csv', 'rt')
try:
    reader = csv.reader(f)
    ii = 0
    for row in reader:
        jj = 0
        for x in row:
            SVal(PHI,ii,jj,int(x))
            jj = jj + 1
        ii = ii + 1
finally:
    f.close()
PHIorig = deepcopy(PHI)
print(' ')
#supply and demand vector
Rmat = matrix(ZZ,rows,rows,0)
for ii in range(rows):
    for jj in range(rows):
        Rmat[ii,jj] = round(normalvariate(0,1000))
SDvect = [0]*rows
for ii in range(rows):
    for jj in range(rows):
        SDvect[ii] = SDvect[ii] + Rmat[ii, jj] - Rmat[jj, ii]
mysum = 0
for ii in range(rows):
    mysum = mysum + SDvect[ii]
SDvect
mysum
#print(' ')
cost = [1]*PHI[1]
corner = [0]
#SDvect = [-6,-14,-18,-6,-46,-14,24,-9,30,59]
colH = [1..PHI[1]]
rowH = [PHI[1]+1..PHI[1]+PHI[0]]
costorig = cost[:]
SDvectorig = SDvect[:]
colHorig = colH[:]
rowHorig = rowH[:]
cornerorig = corner[:]

print(' ')
PHI
print(' ')
printPHI(PHI,cost,SDvect,corner,rowH,colH,'PHI')

#pivotrow = 0
#pivotcol = 0
#XchangeA(PHI, pivotrow, pivotcol, cost, SDvect, corner, rowH, colH)
print(' ')
MS = MinShip(PHI, cost, SDvect, corner, rowH, colH)

print('Simplex method complete after '+str(MS[0])+' Xchange operations')
PHI
print(' ')
printPHI(PHI,cost,SDvect,corner,rowH,colH,'PHI')
print(' ')
print('Shipping vector =')
[round(x,2) for x in MS[1]]

TransMat = matrix(RR,rows,rows,0.0)
for jj in range(cols):
    for ii in range(rows):
        if RVal(PHIorig,ii,jj) > 0: aa = ii
        if RVal(PHIorig,ii,jj) < 0: bb = ii
    TransMat[aa,bb] = MS[1][jj]

print(' ')
TransMatdata = []
for ii in range(rows):
    row = []
    for jj in range(rows+1):
        row.append('7')
    TransMatdata.append(row)
for ii in range(rows):
    TransMatdata[ii][0] = str(rowHorig[ii])
    for jj in range(rows):
        TransMatdata[ii][jj+1] = round(TransMat[ii][jj],2)

hd = [str(x) for x in rowHorig]
cola = ['right']*rows
print(tabulate(TransMatdata, headers=hd, tablefmt="grid", colalign=cola))
RSvect = [0]*rows
for ii in range(rows):
    for jj in range(rows):
        RSvect[ii] = RSvect[ii] + TransMat[ii, jj] - TransMat[jj, ii]
mysum = 0
for ii in range(rows):
    mysum = mysum + RSvect[ii]
print(' ')
print('Receive minus Send vector')
[round(x,2) for x in RSvect]
mysum
#printPHI(PHIorig,costorig,SDvectorig,cornerorig,rowHorig,colHorig,'PHI')
[-8628, 1767, 3650, -5510, -73, -3795, -929, 7069, 1690, 4759] 0 [10, 22, bitarray('0000000000000000000000100000000000000000000001010000000000000000000010101000000000000000000001010100000000000000000000100100000000000000000000101000100000000000000000010001000000000000000000100011000000000000000001000011'), bitarray('1110000000000000000000100111000000000000000001010011100000000000000010101001100000000000000001010101111100000000000000100100001100000000000000101000101000000000000000010001011000000000000000100011010000000000000001000011')] +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | PHI | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | SD | +=======+=====+=====+=====+=====+=====+=====+=====+=====+=====+======+======+======+======+======+======+======+======+======+======+======+======+======+=======+ | 23 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -8628 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 24 | -1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1767 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 25 | 0 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3650 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 26 | 0 | 0 | -1 | 0 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -5510 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 27 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | -1 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | -73 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 28 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | -3795 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 29 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | 0 | 0 | -1 | 0 | 1 | 0 | 0 | 0 | -929 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 30 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | -1 | 0 | 1 | 1 | 0 | 7069 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 31 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | -1 | -1 | 0 | 1 | 1690 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 32 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | -1 | -1 | 4759 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | cost | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ Simplex method complete after 21 Xchange operations [10, 22, bitarray('0000000000000000000000000000010000000000000000011000001100000000000100000000000000000000000010000011100010000000000000000000000010000000001001000010100000000000010000010101110000000000000000000000111110010011000111000000'), bitarray('0000000100000100011000000000110100000010000011011000001100000000000101100000100000000000111010000011100010000000001100000110000010000000001101000010100110000000010000010101110111100100110001110000111110010011000111000000')] +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | PHI | 23 | 24 | 25 | 4 | 5 | 26 | 12 | 17 | 27 | 28 | 6 | 3 | 7 | 29 | 30 | 31 | 9 | 13 | 11 | 20 | 21 | 22 | SD | +=======+======+======+======+=====+=====+======+======+======+======+======+=====+=====+=====+======+======+======+=====+======+======+======+======+======+=======+ | 19 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 929 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 18 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | -1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 3795 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 2 | 1 | 1 | 0 | -1 | -1 | 0 | 0 | 0 | 0 | 0 | -1 | -1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 6861 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 1 | 0 | -1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1767 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 8 | 1 | 1 | 1 | 0 | -1 | 0 | 0 | 0 | 0 | 0 | -1 | -1 | -1 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | 0 | 3211 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 10 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 5510 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 14 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 1 | 0 | -1 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 0 | 3274 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 15 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | -1 | -1 | -1 | 0 | 1 | 761 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 16 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | -1 | -1 | 4759 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 32 | -1 | -1 | -1 | 0 | 0 | -1 | 0 | 0 | -1 | -1 | 0 | 0 | 0 | -1 | -1 | -1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | cost | 3 | 2 | 2 | 1 | 1 | 2 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 30867 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ Shipping vector = [1767.0, 6861.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3211.0, 0.0, 5510.0, 0.0, 0.0, 0.0, 3274.0, 761.0, 4759.0, 0.0, 3795.0, 929.0, 0.0, 0.0, 0.0] +----+------+------+------+------+------+------+------+------+------+------+ | | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | +====+======+======+======+======+======+======+======+======+======+======+ | 23 | 0 | 1767 | 6861 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 24 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 25 | 0 | 0 | 0 | 0 | 3211 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 26 | 0 | 0 | 0 | 0 | 5510 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 27 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3274 | 761 | 4759 | +----+------+------+------+------+------+------+------+------+------+------+ | 28 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3795 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 29 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 929 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 30 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 31 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 32 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ Receive minus Send vector [8628.0, -1767.0, -3650.0, 5510.0, 73.0, 3795.0, 929.0, -7069.0, -1690.0, -4759.0] 0.000000000000000 











load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
load('ProbM2.py')
import numpy
import scipy.optimize
/ext/sage/10.6/local/var/lib/sage/venv-python3.12.5/lib/python3.12/site-packages/scikits/__init__.py:1: DeprecationWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html __import__("pkg_resources").declare_namespace(__name__) /ext/sage/10.6/local/var/lib/sage/venv-python3.12.5/lib/python3.12/site-packages/scikits/__init__.py:1: DeprecationWarning: Deprecated call to `pkg_resources.declare_namespace('scikits')`. Implementing implicit namespace packages (as specified in PEP 420) is preferred to `pkg_resources.declare_namespace`. See https://setuptools.pypa.io/en/latest/references/keywords.html#keyword-namespace-packages __import__("pkg_resources").declare_namespace(__name__) 
load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
import numpy
import scipy.optimize
#from operator import add
#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:
    print("explicit contradiction in premises")
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')
load('Prob.py')
load('ProbM.py')
import numpy
import scipy.optimize
#from operator import add
#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:
    print("explicit contradiction in premises")
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')
load('Prob.py')
load('ProbM.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()
aa = 2**nn # number of atoms in corresponding boolean algebra
    # should abort this with error, if too many variables in polynomial
Mvect = [0]*aa
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
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
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))))
Error in lines 4-4 Traceback (most recent call last): File "/cocalc/lib/python3.11/site-packages/smc_sagews/sage_server.py", line 1250, in execute exec( File "", line 1, in <module> NameError: name 'X1' is not defined 
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
ADJ = zero_matrix(ZZ, N, N)
ADJ = matrix([[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]])

PP = ADJ
print(PP)
print()
PPL = [ADJ]
for ii in range(1,N-1):
    PP=PP*ADJ
    PPL.append(PP)
    print(PP)
    print()

NSP = copy(ADJ)
SPL = copy(ADJ)
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] 
Prob(X1*X8 + X1*X2);
ProbM(X1*X8 + X1*X2);
ProbM(X1*X8 + X1);
nn=3;
MGinv = matrix(QQ,[1]);
for ii in range(1,nn+1):
    MGinv = block_matrix(QQ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False)
print(MGinv);
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, x8 -x1*x8 + x1 [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 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])

GPR = BooleanPolynomialRing(3,'X',order='degneglex')
GPR.gens()
n=len(GPR.gens())
n
GPR.order
GPR.inject_variables()
#Prob(X1*X2 + X1*X0)
a1=ProbM2(X1*X2 + X1*X0)
len(a1)
#P1=ProbM(X1*X2 + X1)
#P1
#P1.coefficients()

mybp=(X0+1)*(X1+1)
mybp=FORALL([NOT(A) for A in GPR.gens()])
mybp
mybp.monomials()
mybp.variables()
#Pmybp=ProbM(mybp)
#Pmybp
#Pmybp.coefficients()

nn=len(GPR.gens());
MGinv = matrix(QQ,[1]);
for ii in range(1,nn+1):
    MGinv = block_matrix(QQ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False)
#print(MGinv);
#(x0+1)*(x1+1)*(x2+1)
CPR = PolynomialRing(QQ,len(GPR.gens()),'x',order='invlex')
CPR.inject_variables()
(FORALL([(1 + x) for x in CPR.gens()])).monomials()
#answer = ProbM2(X1*X2 + X1*X0)
(X0, X1, X2) 3 <built-in method order of sage.rings.polynomial.pbori.pbori.BooleanPolynomialRing object at 0x7f811f893040> Defining X0, X1, X2 Defining x0, x1, x2 2 X0*X1*X2 + X1*X2 + X0*X2 + X0*X1 + X2 + X1 + X0 + 1 [X0*X1*X2, X1*X2, X0*X2, X0*X1, X2, X1, X0, 1] (X0, X1, X2) Defining x0, x1, x2 [x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1] 
load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
load('ProbM2.py')
load('ProbSceme.py')
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])
GPR = BooleanPolynomialRing(3,'X',order='lp')
GPR.gens()
GPR.order
GPR.inject_variables()
Prob(X0*X2 + X1*X0+X0)
ProbM(X0*X2 + X1*X0+X0)
ProbM2(X0*X2 + X1*X0+X0)
#a1=ProbM2(X1*X2+X1)
#len(a1)
#a1[0]
#a1[1]
#mybp=(1+X0)*(1+X1)*(1+X2)#(X1*X2 + X1*X0)
MyBP=FORALL([(GPR(1) + x) for x in GPR.gens()])
MyBP
MyBP.parent()
MyBP.parent().gens()
len(MyBP.parent().gens())
MyBPMons=MyBP.monomials()
MyBP.monomials()
2**len(MyBP.parent().gens())
nn=len(GPR.gens());
MG=matrix(ZZ,[1])
MGinv = matrix(ZZ,[1]);
for ii in range(1,nn+1):
    MG = block_matrix(ZZ,2,2,[MG,MG,0,MG],subdivide=False)
    MGinv = block_matrix(ZZ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False)
print(MGinv);
#(x0+1)*(x1+1)*(x2+1)
CPR = PolynomialRing(QQ,len(GPR.gens()),'x',order='invlex')
CPR.inject_variables()
MyMons=(FORALL([(CPR(1) + x) for x in CPR.gens()])).monomials()
MyMons
#Mvect = vector([0,0,0,0,1,1,0,0])
#Avect = (MGinv*Mvect)%2
#Pvect = MGinv*Avect
#Pvect
Mvect = vector([0,1,0,0,1,0,0,0])
Avect = (MGinv*Mvect)%2
Pvect = MGinv*Avect
Mvect
Avect
Pvect
mystr='XcGGt1u'
ii=0
if mystr[ii]==mystr[ii].upper():
    mystr[ii].lower()
else:
    mystr[ii].upper()
mystr
#MG=MGinv%2
MG #*
(MyBP.monomials())[0]
MonVec= vector([x for x in MyBP.monomials()])
MonVec
MonVec2= vector([x for x in MyMons])
MonVec2
#MG*MonVec
MonVec*MG
(MonVec*MG)*MG
(X0, X1, X2) <built-in method order of sage.rings.polynomial.pbori.pbori.BooleanPolynomialRing object at 0x7fb60d09d9c0> Defining X0, X1, X2 Defining x0, x1, x2 2*x0*x1*x2 - x0*x2 - x0*x1 + x0 Defining x0, x1, x2 2*x0*x1*x2 - x0*x1 - x0*x2 + x0 Defining x0, x1, x2 ((2, 0, -1, 0, -1, 0, 1, 0), (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1)) X0*X1*X2 + X0*X1 + X0*X2 + X0 + X1*X2 + X1 + X2 + 1 Boolean PolynomialRing in X0, X1, X2 (X0, X1, X2) 3 [X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1] 8 [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 0 0 0 0 0 1] Defining x0, x1, x2 [x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1] (0, 1, 0, 0, 1, 0, 0, 0) (0, 1, 0, 0, 1, 0, 0, 0) (-2, 1, 0, 0, 1, 0, 0, 0) 'x' 'XcGGt1u' [1 1 1 1 1 1 1 1] [0 1 0 1 0 1 0 1] [0 0 1 1 0 0 1 1] [0 0 0 1 0 0 0 1] [0 0 0 0 1 1 1 1] [0 0 0 0 0 1 0 1] [0 0 0 0 0 0 1 1] [0 0 0 0 0 0 0 1] X0*X1*X2 (X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1) (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1) (X0*X1*X2, X0*X1*X2 + X0*X1, X0*X1*X2 + X0*X2, X0*X1*X2 + X0*X1 + X0*X2 + X0, X0*X1*X2 + X1*X2, X0*X1*X2 + X0*X1 + X1*X2 + X1, X0*X1*X2 + X0*X2 + X1*X2 + X2, X0*X1*X2 + X0*X1 + X0*X2 + X0 + X1*X2 + X1 + X2 + 1) (X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1) 
v=ProbSceme(3)
v
a=0
for ii in v:
    a=a+ii
a
#mg=MGinv%2
MG
v*MG
Mvect
(MGinv*Mvect)%2
Mvect = vector([0,0,0,0,1,1,0,0])
Mvect
(MGinv*Mvect)%2
(MG*Mvect)%2
ProbM(X1+X2*X0)
A=ProbM2(X1+X2*X0)
A
(0.216897236248995, 0.125181737169618, 0.123043301373761, 0.0748467738425043, 0.125643838286659, 0.0283746925135870, 0.122649911882767, 0.183362508682108) 1.00000000000000 [1 1 1 1 1 1 1 1] [0 1 0 1 0 1 0 1] [0 0 1 1 0 0 1 1] [0 0 0 1 0 0 0 1] [0 0 0 0 1 1 1 1] [0 0 0 0 0 1 0 1] [0 0 0 0 0 0 1 1] [0 0 0 0 0 0 0 1] (0.216897236248995, 0.342078973418613, 0.339940537622756, 0.539969048634878, 0.342541074535654, 0.496097504218859, 0.588234287792182, 1.00000000000000) (0, 1, 0, 0, 1, 0, 0, 0) (0, 1, 0, 0, 1, 0, 0, 0) (0, 0, 0, 0, 1, 1, 0, 0) (0, 1, 0, 0, 0, 1, 0, 0) (0, 1, 0, 0, 0, 1, 0, 0) Defining x0, x1, x2 -2*x0*x1*x2 + x0*x2 + x1 Defining x0, x1, x2 ((-2, 0, 1, 0, 0, 1, 0, 0), (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1)) 
myBP=X0*X2+ X1*X0+X0
myBP
BPvars = myBP.variables()
BPmons = myBP.monomials()
BPvars
BPmons
nn = 3
atoms = 2**nn # number of atoms in corresponding boolean algebra
    # should abort this with error, if too many variables in polynomial
Mvect = [0]*atoms
Mvect
len(BPmons)
for jj in range(len(BPmons)):
    mm = 0
    BPmons[jj].variables()
    for kk in range(len(BPmons[jj].variables())):
        #BPvars.index(BPmons[jj].variables()[kk])
        mm = mm + 2**BPvars.index(BPmons[jj].variables()[kk])
    'mm='
    mm
    Mvect[mm] = 1
Mvect
Mvect.reverse()
Mvect
Mvect = vector(Mvect)
MGinv = matrix(ZZ,[1])
for ii in range(1,nn+1):
    MGinv = block_matrix(ZZ,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
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
MGinv
Avect = (MGinv*Mvect)%2
Pvect = MGinv*Avect
Avect
Pvect
ProbM(myBP)
ProbM2(myBP)
Prob(myBP)
X0*X1 + X0*X2 + X0 (X0, X1, X2) [X0*X1, X0*X2, X0] [0, 0, 0, 0, 0, 0, 0, 0] 3 (X0, X1) 'mm=' 3 (X0, X2) 'mm=' 5 (X0,) 'mm=' 1 [0, 1, 0, 1, 0, 1, 0, 0] [0, 0, 1, 0, 1, 0, 1, 0] [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 0 0 0 0 0 1] (1, 0, 0, 0, 0, 0, 1, 0) (2, 0, -1, 0, -1, 0, 1, 0) Defining x0, x1, x2 2*x0*x1*x2 - x0*x1 - x0*x2 + x0 Defining x0, x1, x2 ((2, 0, -1, 0, -1, 0, 1, 0), (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1)) Defining x0, x1, x2 2*x0*x1*x2 - x0*x2 - x0*x1 + x0 
myBP=X0*X1+X0*X2+X0
myBP
if type(myBP) == sage.rings.integer.Integer:
    if myBP%2 == 0:
        myQP = QQ(0)
    else:
        myQP = QQ(1)
    #myBP.variables()
mystr = str(myBP.variables()) # converting all uppers to lower and lowers to uppercase
    #print(mystr)
NewGens=""
for ii in range(len(mystr)):
    if mystr[ii] == mystr[ii].upper():
        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(',','')
    #NewGens.replace(', ',',')
NewGens.replace(' ','')
    #print(mystr)
    #print(NewGens)
CPR = PolynomialRing(QQ,len(myBP.variables()),NewGens,order='invlex')
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(myBP.monomials())):
    mm = 0
    for kk in range(len(myBP.monomials()[jj].variables())):
        mm = mm + 2**myBP.variables().index(myBP.monomials()[jj].variables()[kk])
    Mvect[mm] = 1
Mvect.reverse()
mypoly=myBP.parent(1)
for ii in range(0,len(CPR.gens())):
    mypoly = mypoly*(myBP.parent(1) + myBP.parent().gen(ii))
#allmons=((1+X0)*(1+X1)*(1+X2)).monomials()
allmons = mypoly.monomials()
'Mvect'
Mvect
    #print Mvect
    # 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
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
Avect = (MGinv*Mvect)%2
'Avect'
Avect
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
myQP = Pvect*vector(newpoly.monomials()) # inner product of the two vectors to create prob poly
'Pvect'
Pvect
vector(newpoly.monomials())
myQP
Mvect2 = [0]*atoms
for y in myBP.monomials():
    for jj in range(len(allmons)):
        if y==allmons[jj]:
            Mvect2[jj]=1
Mvect2
myBP.parent().term_order()
X0*X1 + X0*X2 + X0 'x0,x1,x2' Defining x0, x1, x2 'Mvect' [0, 0, 1, 0, 1, 0, 1, 0] 'Avect' (1, 0, 0, 0, 0, 0, 1, 0) 'Pvect' (2, 0, -1, 0, -1, 0, 1, 0) (x0*x1*x2, x1*x2, x0*x2, x2, x0*x1, x1, x0, 1) 2*x0*x1*x2 - x0*x2 - x0*x1 + x0 [0, 1, 1, 1, 0, 0, 0, 0] lp term order 
'0'*(7-len(bin(0)[2:]))+bin(0)[2:]
bin(1)[2:]
bin(4)[2:]
n=0
n.nbits()
''
#n.binary()
#'0' * (3-n.nbits()) + n.binary()
nn=4
n=0
n.binary()
for ii in range(2**nn):
    if n==0:
        '0'*nn
    else:
        '0'*(nn-n.nbits())+n.binary()
    n=n+1
'0000000' '1' '100' 0 '' '0' '0000' '0001' '0010' '0011' '0100' '0101' '0110' '0111' '1000' '1001' '1010' '1011' '1100' '1101' '1110' '1111' 
from tabulate import tabulate
load('ProbRecurse.py')
load('Prob.py')
load('ProbM.py')
load('ProbM2.py')
load('ProbScheme.py')
load('IndProbSchema.py')
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])
nn = 3 # of generators of Boolean Algebra
print('Creating a Boolean polynomial ring with '+str(nn)+ ' generators which will act as two state random variables.')
print(' ')
aa = 2**nn # number of atoms in the Boolean Algebra
BPR = BooleanPolynomialRing(nn,'X',order='lp')
print('The generators of the Boolean Polynomial Ring, BPR')
BPR.inject_variables()
BPR.term_order()
mypoly=BPR(1)
for ii in range(0,len(BPR.gens())):
    mypoly=mypoly*(BPR(1) + BPR.gen(ii))
print(' ')
CPR = PolynomialRing(QQ,nn,'x',order='lp') # probability polynomials over QQ
print('Generators for the probability polynomials over QQ')
CPR.inject_variables()
CPR.term_order()
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(' ')
IPS=IndProbSchema(nn)
print("Probabilities (heads) of the "+str(nn)+" generators taken as distinguishable (unfair) coin tosses. These are independent probabilities, and do not sum to 1.0")
prbheads = [round(x,6) for x in IPS[0]]
hdata = []
headsdata = ['prob. of heads']
headers=['generators']
cola=['left']
trials = 30
cointoss = matrix(ZZ,30,nn,0)
frequency = vector([0.0]*nn)
for ii in range(nn):
    headsdata.append(str(round(prbheads[ii],6)))
    headers.append('P(' + str(BPR.gens()[ii]) + ') = ' + str(CPR.gens()[ii]))
    cola.append('center')
hdata.append(headsdata)
for ii in range(trials):
    trialrow = ['coin toss '+ str(ii+1)]
    for jj in range(nn):
        if random() <= prbheads[jj]:
            cointoss[ii,jj] = 1
            trialrow.append('1')
        else:
            trialrow.append('0')
        frequency[jj] = frequency[jj] + cointoss[ii,jj]
    hdata.append(trialrow)
frequency = frequency/trials
trialrow = ['freq. of heads']
for jj in range(nn):
    trialrow.append(round(frequency[jj],6))
hdata.append(trialrow)
print(tabulate(hdata, headers=headers, tablefmt="grid", colalign=cola))
print(' ')
print('The number of atoms is '+str(aa)+' = 2**'+str(nn)+' generators. The generator probabilities induce the following probability scheme for the '+str(aa)+' atoms.')
print('Aprob')
[round(x,6) for x in IPS[1]]
print(' ')
mysum=0
for ii in range(aa):
    mysum=mysum+IPS[1][ii]
print("Total probability (just checking)")
round(mysum,6)
print(' ')
print('Table of atoms and monomials. With the exception of the 0th atom and monomial the rest are not equal.')
AMdata = []
for ii in range(aa):
    rowAM = []
    mystring = '0'*(nn-len(bin(ii)[2:]))+bin(ii)[2:]
    mystring2 = ''
    for jj in range(nn):
        mystring2 = mystring2 + '*(' + str(BPR.gen(jj)) + '+' + mystring[jj] +')'
    rowAM.append(str(ii))
    rowAM.append(mystring)
    rowAM.append(mystring2[1:])
    rowAM.append(str(mypoly.monomials()[ii]))
    AMdata.append(rowAM)
    #print(str(ii)+'  '+ mystring +'  ' + mystring2[1:] + '  ' + str(mypoly.monomials()[ii]))
print(tabulate(AMdata, headers=('#','bin','atoms', 'monomials'), tablefmt="grid", colalign=('center','center','center','left')))
print(' ')
###############################################################################################
# should abort this with error, if too many variables in polynomial
myBP = X0*X1 + X0*X2 + X0 + 1
###############################################################################################
print('An example Boolean polynomial:')
print(' ')
print('myBP')
myBP
print(' ')
Mvect = [0]*aa
for y in myBP.monomials():
    for jj in range(len(mypoly.monomials())):
        if y==mypoly.monomials()[jj]:
            Mvect[jj]=1
MGAinv = matrix(ZZ,[1])
for ii in range(1,nn+1):
    MGAinv = block_matrix(ZZ,2,2,[MGAinv,-MGAinv,0,MGAinv],subdivide=False) # creating inverse of multigrade operator 'AND' matrix MGA
        # modulo 2 it's an involution however...
        # MGA and MGAinv are upper triangular, MGA is like Sierpinski triangle see
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
        # over ZZ it is totally unimodular
MGA = MGAinv%2
Mvect=vector(Mvect)
print('Monomial vector')
print('MonVec')
print(str(vector(mypoly.monomials())))
print(' ')
print('Polynomial vector')
print('Mvect = (MGA*Avect)mod 2, Mvect is 0-1 vector with a 1 corresponding to each monomial in the Bollean polynomial, 0 otherwise')
Mvect
print(' ')
Avect = MGA*Mvect%2
print('Atom vector')
print('Avect = (MGA*Mvect)mod 2')
Avect
print(' ')
print("Multigrade And matrix, it's totally unimodular over ZZ, and it's an involution(self inverse) mod 2")
print('MGA =')
MGA
print(' ')
print('Inverse of Multigrade And matrix in ZZ, QQ, or RR')
print('MGAinv')
MGAinv
Pvect = MGAinv*Avect
print(' ')
Mprob=IPS[1]*MGA
print('Monomial probabilities(which do not sum to 1.0)')
print('but the last term is the sum of all atom probabilities = 1.0')
print('Mprob = Aprob*MGA')
[round(x,6) for x in Mprob]
print(' ')
print('Generators for the probability polynomial or formula over RR')
print(' ')
print('Pvect = MGAinv*Avect, the coefficients of the probability polynomial or formula')
[round(x,6) for x in Pvect]
print(' ')
print('Monoms, monomials over RR')
Monoms = vector(newpoly.monomials())
print(' ')
print('Probability polynomial when Boolean Ring generators are independent two state random variables')
print('ProbPoly = Pvect*Monoms')
Pvect*Monoms
print(' ')
ProbExpr = ''
for ii in range(aa):
    if Pvect[ii] > 1: mystr = ' + '+str(Pvect[ii])+'*'
    if Pvect[ii] == 1: mystr = str(' + ')
    if Pvect[ii] == -1: mystr = str(' - ')
    if Pvect[ii] < -1: mystr = ' - '+str(abs(Pvect[ii]))+'*'
    if Pvect[ii] != 0: ProbExpr = ProbExpr + mystr+'P('+ str(mypoly.monomials()[ii])+')'
if ProbExpr[1] == '+' : ProbExpr = ProbExpr[3:]
if ProbExpr[1] == '-' : ProbExpr = '-' + ProbExpr[3:]
print('Probability formula whether or not the generators are independent,')
print('P('+str(myBP)+') = ' + ProbExpr)
print(' ')
print("In this example with independent generators the numerical probability (using the dot product of vectors), P("+str(myBP)+") = ")
print('Pvect*Mprob =')
Pvect*Mprob
print('or Avect*Aprob =')
IPS[1]*Avect
print(' ')
print('An example where the generators are not independent is constructed by creating a probability schema for the atoms arbitrarily')
PS = ProbScheme(nn)
print('Probability schema for the atoms, PS =')
[round(x,6) for x in PS]
print(' ')
MP=PS*MGA
print('Monomial probabilities')
print('MP = PS*MGA =')
[round(x,6) for x in MP]
print(' ')
print('In this example with non-independent generators the numerical probability, P('+str(myBP)+') = ')
print('Pvect*MP =')
Pvect*MP
print('Avect*PS =')
Avect*PS
Creating a Boolean polynomial ring with 3 generators which will act as two state random variables. The generators of the Boolean Polynomial Ring, BPR Defining X0, X1, X2 lp term order Generators for the probability polynomials over QQ Defining x0, x1, x2 lp term order Probabilities (heads) of the 3 generators taken as distinguishable (unfair) coin tosses. These are independent probabilities, and do not sum to 1.0 +----------------+--------------+--------------+--------------+ | generators | P(X0) = x0 | P(X1) = x1 | P(X2) = x2 | +================+==============+==============+==============+ | prob. of heads | 0.970432 | 0.994518 | 0.119172 | +----------------+--------------+--------------+--------------+ | coin toss 1 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 2 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 3 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 4 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 5 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 6 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 7 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 8 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 9 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 10 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 11 | 1 | 1 | 1 | +----------------+--------------+--------------+--------------+ | coin toss 12 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 13 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 14 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 15 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 16 | 1 | 1 | 1 | +----------------+--------------+--------------+--------------+ | coin toss 17 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 18 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 19 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 20 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 21 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 22 | 1 | 1 | 1 | +----------------+--------------+--------------+--------------+ | coin toss 23 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 24 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 25 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 26 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 27 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 28 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 29 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | coin toss 30 | 1 | 1 | 0 | +----------------+--------------+--------------+--------------+ | freq. of heads | 1 | 1 | 0.1 | +----------------+--------------+--------------+--------------+ The number of atoms is 8 = 2**3 generators. The generator probabilities induce the following probability scheme for the 8 atoms. Aprob [0.115014, 0.850098, 0.000634, 0.004686, 0.003504, 0.025901, 1.9e-05, 0.000143] Total probability (just checking) 1.0 Table of atoms and monomials. With the exception of the 0th atom and monomial the rest are not equal. +-----+-------+----------------------+-------------+ | # | bin | atoms | monomials | +=====+=======+======================+=============+ | 0 | 000 | (X0+0)*(X1+0)*(X2+0) | X0*X1*X2 | +-----+-------+----------------------+-------------+ | 1 | 001 | (X0+0)*(X1+0)*(X2+1) | X0*X1 | +-----+-------+----------------------+-------------+ | 2 | 010 | (X0+0)*(X1+1)*(X2+0) | X0*X2 | +-----+-------+----------------------+-------------+ | 3 | 011 | (X0+0)*(X1+1)*(X2+1) | X0 | +-----+-------+----------------------+-------------+ | 4 | 100 | (X0+1)*(X1+0)*(X2+0) | X1*X2 | +-----+-------+----------------------+-------------+ | 5 | 101 | (X0+1)*(X1+0)*(X2+1) | X1 | +-----+-------+----------------------+-------------+ | 6 | 110 | (X0+1)*(X1+1)*(X2+0) | X2 | +-----+-------+----------------------+-------------+ | 7 | 111 | (X0+1)*(X1+1)*(X2+1) | 1 | +-----+-------+----------------------+-------------+ An example Boolean polynomial: myBP X0*X1 + X0*X2 + X0 + 1 Monomial vector MonVec (X0*X1*X2, X0*X1, X0*X2, X0, X1*X2, X1, X2, 1) Polynomial vector Mvect = (MGA*Avect)mod 2, Mvect is 0-1 vector with a 1 corresponding to each monomial in the Bollean polynomial, 0 otherwise (0, 1, 1, 1, 0, 0, 0, 1) Atom vector Avect = (MGA*Mvect)mod 2 (0, 1, 1, 0, 1, 1, 1, 1) Multigrade And matrix, it's totally unimodular over ZZ, and it's an involution(self inverse) mod 2 MGA = [1 1 1 1 1 1 1 1] [0 1 0 1 0 1 0 1] [0 0 1 1 0 0 1 1] [0 0 0 1 0 0 0 1] [0 0 0 0 1 1 1 1] [0 0 0 0 0 1 0 1] [0 0 0 0 0 0 1 1] [0 0 0 0 0 0 0 1] Inverse of Multigrade And matrix in ZZ, QQ, or RR MGAinv [ 1 -1 -1 1 -1 1 1 -1] [ 0 1 0 -1 0 -1 0 1] [ 0 0 1 -1 0 0 -1 1] [ 0 0 0 1 0 0 0 -1] [ 0 0 0 0 1 -1 -1 1] [ 0 0 0 0 0 1 0 -1] [ 0 0 0 0 0 0 1 -1] [ 0 0 0 0 0 0 0 1] Monomial probabilities(which do not sum to 1.0) but the last term is the sum of all atom probabilities = 1.0 Mprob = Aprob*MGA [0.115014, 0.965112, 0.115648, 0.970432, 0.118519, 0.994518, 0.119172, 1.0] Generators for the probability polynomial or formula over RR Pvect = MGAinv*Avect, the coefficients of the probability polynomial or formula [-2.0, 1.0, 1.0, -1.0, 0.0, 0.0, 0.0, 1.0] Monoms, monomials over RR Probability polynomial when Boolean Ring generators are independent two state random variables ProbPoly = Pvect*Monoms -2*x0*x1*x2 + x0*x1 + x0*x2 - x0 + 1 Probability formula whether or not the generators are independent, P(X0*X1 + X0*X2 + X0 + 1) = -2*P(X0*X1*X2) + P(X0*X1) + P(X0*X2) - P(X0) + P(1) In this example with independent generators the numerical probability (using the dot product of vectors), P(X0*X1 + X0*X2 + X0 + 1) = Pvect*Mprob = 0.880299732641442 or Avect*Aprob = 0.880299732641442 An example where the generators are not independent is constructed by creating a probability schema for the atoms arbitrarily Probability schema for the atoms, PS = [0.036326, 0.026352, 0.219059, 0.037976, 0.113768, 0.101268, 0.225121, 0.240132] Monomial probabilities MP = PS*MGA = [0.036326, 0.062677, 0.255384, 0.319712, 0.150093, 0.277713, 0.594273, 1.0] In this example with non-independent generators the numerical probability, P(X0*X1 + X0*X2 + X0 + 1) = Pvect*MP = 0.925698661517172 Avect*PS = 0.925698661517172 
from tabulate import tabulate
#
nn=7
TRM=matrix(RR,nn,nn)
for ii in [1..nn]:
    for jj in [1..nn]:
        TRM[ii-1,jj-1]=int(round(abs(normalvariate(0,1000))))

NET = TRM-TRM.transpose()
for ii in [1..nn]:
    for jj in [1..nn]:
        if NET[ii-1,jj-1] < 0.0: NET[ii-1,jj-1] = 0.0

b=[0.0]*nn
x=[0.0]*(nn-1)
for kk in [1..nn]:
    for ii in [1..nn]:
        b[kk-1] = b[kk-1]-TRM[ii-1,kk-1]+TRM[kk-1,ii-1]

m = 2*nn
for ii in [1..nn-1]:
    for jj in [1..nn-1]:
        x[ii-1] = b[jj-1] + x[ii-1]
        if ii == jj: x[ii-1] = b[jj-1] + x[ii-1]
    x[ii-1] = x[ii-1]/m

LMD = matrix(RR,nn,nn)
for ii in [2..nn]:
    for jj in [1..ii-1]:
        for kk in [1..nn-1]:
            if ii==kk: LMD[ii-1,jj-1] = LMD[ii-1,jj-1] + x[kk-1]
            if jj==kk: LMD[ii-1,jj-1] = LMD[ii-1,jj-1] - x[kk-1]
        LMD[jj-1,ii-1] = -LMD[ii-1,jj-1]

print(' ')
NETL = LMD-LMD.transpose()
for ii in range(nn):
    for jj in range(nn):
        if NETL[ii,jj] < 0.0: NETL[ii,jj] = 0.0
CYC = TRM-LMD
CYCS = (CYC+CYC.transpose())/2 # symmetric cycle
CYCK = (CYC-CYC.transpose())/2 # skew-symmetric cycle
#rep=matrix(ZZ,nn,1,[1]*nn)
recvTRM = matrix(RR,nn,1,0)
sendTRM = matrix(RR,nn,1,0)
recvNET = matrix(RR,nn,1,0)
sendNET = matrix(RR,nn,1,0)
recvLMD = matrix(RR,nn,1,0)
sendLMD = matrix(RR,nn,1,0)
recvNETL = matrix(RR,nn,1,0)
sendNETL = matrix(RR,nn,1,0)
recvCYCS = matrix(RR,nn,1,0)
sendCYCS = matrix(RR,nn,1,0)
recvCYCK = matrix(RR,nn,1,0)
sendCYCK = matrix(RR,nn,1,0)
for ii in range(nn):
    recvTRM[ii,0] = sum(TRM[ii])
    recvNET[ii,0] = sum(NET[ii])
    recvLMD[ii,0] = sum(LMD[ii])
    recvNETL[ii,0] = sum(NETL[ii])
    recvCYCS[ii,0] = sum(CYCS[ii])
    recvCYCK[ii,0] = sum(CYCK[ii])
    sendTRM[ii,0] = sum(TRM.transpose()[ii])
    sendNET[ii,0] = sum(NET.transpose()[ii])
    sendLMD[ii,0] = sum(LMD.transpose()[ii])
    sendNETL[ii,0] = sum(NETL.transpose()[ii])
    sendCYCS[ii,0] = sum(CYCS.transpose()[ii])
    sendCYCK[ii,0] = sum(CYCK.transpose()[ii])

Nodes = [' ']
for ii in range(nn):
    Nodes.append('N'+str(ii))
Nodes.append(' ')
Nodes.append('R')
Nodes.append(' ')
Nodes.append('R-S')
TRMdata=[]
NETdata=[]
LMDdata=[]
NETLdata=[]
CYCSdata = []
CYCKdata = []
for ii in range(nn):
    rowT = [Nodes[ii+1]]
    rowN = [Nodes[ii+1]]
    rowL = [Nodes[ii+1]]
    rowNL = [Nodes[ii+1]]
    rowCS = [Nodes[ii+1]]
    rowCK = [Nodes[ii+1]]
    for jj in range(nn):
        rowT.append('{0:.2f}'.format(TRM[ii,jj]))
        rowN.append('{0:.2f}'.format(NET[ii,jj]))
        rowL.append('{0:.2f}'.format(LMD[ii,jj]))
        rowNL.append('{0:.2f}'.format(NETL[ii,jj]))
        rowCS.append('{0:.2f}'.format(CYCS[ii,jj]))
        rowCK.append('{0:.2f}'.format(CYCK[ii,jj]))
    rowT.append(' ')
    rowT.append('{0:.2f}'.format(recvTRM[ii,0]))
    rowT.append(' ')
    rowT.append('{0:.2f}'.format((recvTRM-sendTRM)[ii,0]))

    rowN.append(' ')
    rowN.append('{0:.2f}'.format(recvNET[ii,0]))
    rowN.append(' ')
    rowN.append('{0:.2f}'.format((recvNET-sendNET)[ii,0]))

    rowL.append(' ')
    rowL.append('{0:.2f}'.format(recvLMD[ii,0]))
    rowL.append(' ')
    rowL.append('{0:.2f}'.format((recvLMD-sendLMD)[ii,0]))

    rowNL.append(' ')
    rowNL.append('{0:.2f}'.format(recvNETL[ii,0]))
    rowNL.append(' ')
    rowNL.append('{0:.2f}'.format((recvNETL-sendNETL)[ii,0]))

    rowCS.append(' ')
    rowCS.append('{0:.2f}'.format(recvCYCS[ii,0]))
    rowCS.append(' ')
    rowCS.append('{0:.2f}'.format((recvCYCS-sendCYCS)[ii,0]))

    rowCK.append(' ')
    rowCK.append('{0:.2f}'.format(recvCYCK[ii,0]))
    rowCK.append(' ')
    rowCK.append('{0:.2f}'.format((recvCYCK-sendCYCK)[ii,0]))

    TRMdata.append(rowT)
    NETdata.append(rowN)
    LMDdata.append(rowL)
    NETLdata.append(rowNL)
    CYCSdata.append(rowCS)
    CYCKdata.append(rowCK)

TRMdata.append([' ']*(nn+5))
TRMdata.append([' ']*(nn+5))
TRMdata[nn+1][0] = 'S'
TRMdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvTRM)))
TRMdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvTRM-sendTRM)))
NETdata.append([' ']*(nn+5))
NETdata.append([' ']*(nn+5))
NETdata[nn+1][0] = 'S'
NETdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvNET)))
NETdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvNET-sendNET)))
LMDdata.append([' ']*(nn+5))
LMDdata.append([' ']*(nn+5))
LMDdata[nn+1][0] = 'S'
LMDdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvLMD)))
LMDdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvLMD-sendLMD)))
NETLdata.append([' ']*(nn+5))
NETLdata.append([' ']*(nn+5))
NETLdata[nn+1][0] = 'S'
NETLdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvNETL)))
NETLdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvNETL-sendNETL)))
CYCSdata.append([' ']*(nn+5))
CYCSdata.append([' ']*(nn+5))
CYCSdata[nn+1][0] = 'S'
CYCSdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvCYCS)))
CYCSdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvCYCS-sendCYCS)))
CYCKdata.append([' ']*(nn+5))
CYCKdata.append([' ']*(nn+5))
CYCKdata[nn+1][0] = 'S'
CYCKdata[nn+1][nn+2] = '{0:.2f}'.format(sum(vector(recvCYCK)))
CYCKdata[nn+1][nn+4] = '{0:.2f}'.format(sum(vector(recvCYCK-sendCYCK)))
SKEW = TRM-TRM.transpose()
SKEWdata = []
for ii in range(nn):
    TRMdata[nn+1][ii+1]='{0:.2f}'.format(sendTRM[ii,0])
    NETdata[nn+1][ii+1]='{0:.2f}'.format(sendNET[ii,0])
    LMDdata[nn+1][ii+1]='{0:.2f}'.format(sendLMD[ii,0])
    NETLdata[nn+1][ii+1]='{0:.2f}'.format(sendNETL[ii,0])
    CYCSdata[nn+1][ii+1]='{0:.2f}'.format(sendCYCS[ii,0])
    CYCKdata[nn+1][ii+1]='{0:.2f}'.format(sendCYCK[ii,0])
    rowSK = [Nodes[ii+1]]
    for jj in range(nn):
        rowSK.append('{0:.2f}'.format(SKEW[ii][jj]))
    SKEWdata.append(rowSK)

cola = ["right" for x in Nodes]
cola[0] = "center"
Nodes[0]='TRM'
mystring = 'The square matrix, TRM, of data in the following table represents the flow of goods and services in a complete directed multigraph of '+str(nn)+' nodes as valued in a single currency. The nodes of the graph might represent the nations of a small world. The data can be interpreted in either of two ways. Here we choose the interpretation that the value, X, in the cell corresponding to row 0 and column 1 is the amount received by node 0 from node 1. The value, Y, corresponding to row 1 and column 0 is the amount received by node 1 from node 0. One of these values will generally be greater than(though possibly equal to) the other. The net effect is a positive flow of |X-Y| in one direction or the other. A negative entry represents a flow in the direction opposite to a given edge in the graph. The column, R, is the sum along each row, and represents the total amount received from the network for each node. The row, S, represents the amount sent into the network for each node. The column, R-S, represents the net inflow or outflow of money for each node during a given time step. The data for this table was generated with a random variable from a half-normal distribution, i.e. |Z| where Z ~ N(0,sigma^2). TRM is my abbreviation for transaction matrix.'
print(mystring)
print(tabulate(TRMdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
cola2 = ["right" for ii in range(nn)]
cola2[0] = 'center'
Nodes[0] = 'SKEW'
mystring = "To see the net effect of money transfers from node to node, we first form a skew-symmetric matrix(each entry is the negative of it's mirror reflection through the main diagonal) from TRM given by SKEW = TRM - TRM.transpose()."
print(mystring)
print(tabulate(SKEWdata, headers=Nodes, tablefmt="grid", colalign=cola2))
print(' ')
Nodes[0]='NET'
mystring = "Next, we set the negative entries of SKEW to 0, and keep the rest to obtain the matrix, NET. For instance if you send me $10, and I send you $5, then the net effect is you sent me $5, and I sent you $0, we're just expressing everything in non-negative numbers. We could apply x*H(x) where H is the Heaviside step function to each entry of SKEW to obtain NET. Note that the R-S column vector for NET has not changed from that of TRM. We also note that for any square matrix, the sum over the R vector equals the sum over the S vector which is just the sum of all entries in the matrix, and therefore the R-S vector always sums to 0."
print(mystring)
print(tabulate(NETdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='LMD'
mystring = 'Closer analysis using some graph theory and linear algebra allows us to express the net effect of all transactions in TRM efficiently as the projection of TRM into the cocycle space of the complete directed multigraph on N vertices and 2*N*(N-1) edges. The dimension of the cocycle space is N-1. Projecting TRM thusly we obtain the skew-symmetric matrix, LMD. I called it this for no good reason having inherited it from a Visual Basic version I wrote a number of years ago. Note here also, the R-S vector remains unchanged from TRM. As a further note the numbers here are rounded to two decimal places whereas they have many more significant digits so the numbers may not add up exactly, if checking.'
print(mystring)
print(tabulate(LMDdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='NETL'
mystring = 'To express the transactions in LMD as non-negative numbers only form SKEW = LMD - LMD.transpose, and set the negatives to 0. We obtain the matrix NETL. Again R-S is unchanged.'
print(mystring)
print(tabulate(NETLdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='CYCS'
mystring = "Any matrix where the R-S vector is the 0-vector(not merely summing to 0, but every entry is 0) lies in the cycle space of the complete directed graph. We can project TRM into the cycle space to obtain a matrix, CYC. This can be obtained by simply subtracting LMD from TRM, thus CYC = TRM - LMD. CYC can be broken down further into symmetric and skew-symmetric matrices, CYCS and CYCK respectively where CYCS = (CYC + CYC.transpose())/2 and CYCK = (CYC - CYC.transpose())/2. We have TRM = LMD + CYCS + CYCK."
print(mystring)
print(tabulate(CYCSdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
Nodes[0]='CYCK'
mystring = 'Die Quantenmechanik ist sehr achtung-gebietend. Aber eine innere Stimme sagt mir, daß das doch nicht der wahre Jakob ist. Die Theorie liefert viel, aber dem Geheimnis des Alten bringt sie uns kaum näher. Jedenfalls bin ich überzeugt, daß der nicht würfelt.'
mystring = "The space of all N×N square matrices' dimension is N^2. The cocycle space is dimension N-1, so that the cycle space is dimension N^2 - N + 1. The cycle space is further broken down into the symmetric cycles and the skew-symmetric cycles. Every symmetric matrix is a cycle, and the dimension of the symmetric matrices is N*(N+1)/2. The dimension of the skew-symmetric matrices is N*(N-1)/2 leaving the dimension of the skew-symmetric cycles (N-1)*(N-2)/2."
print(mystring)
print(tabulate(CYCKdata, headers=Nodes, tablefmt="grid", colalign=cola))
print(' ')
The square matrix, TRM, of data in the following table represents the flow of goods and services in a complete directed multigraph of 7 nodes as valued in a single currency. The nodes of the graph might represent the nations of a small world. The data can be interpreted in either of two ways. Here we choose the interpretation that the value, X, in the cell corresponding to row 0 and column 1 is the amount received by node 0 from node 1. The value, Y, corresponding to row 1 and column 0 is the amount received by node 1 from node 0. One of these values will generally be greater than(though possibly equal to) the other. The net effect is a positive flow of |X-Y| in one direction or the other. A negative entry represents a flow in the direction opposite to a given edge in the graph. The column, R, is the sum along each row, and represents the total amount received from the network for each node. The row, S, represents the amount sent into the network for each node. The column, R-S, represents the net inflow or outflow of money for each node during a given time step. The data for this table was generated with a random variable from a half-normal distribution, i.e. |Z| where Z ~ N(0,sigma^2). TRM is my abbreviation for transaction matrix. +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | TRM | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +=======+=========+==========+=========+=========+=========+=========+=========+=====+==========+=====+==========+ | N0 | 1606.00 | 877.00 | 343.00 | 764.00 | 731.00 | 495.00 | 1351.00 | | 6167.00 | | -859.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N1 | 459.00 | 1154.00 | 678.00 | 817.00 | 719.00 | 367.00 | 2389.00 | | 6583.00 | | -4081.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N2 | 1160.00 | 1838.00 | 539.00 | 295.00 | 672.00 | 582.00 | 193.00 | | 5279.00 | | -1926.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N3 | 1100.00 | 2073.00 | 888.00 | 1065.00 | 330.00 | 1092.00 | 542.00 | | 7090.00 | | 1607.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N4 | 588.00 | 1241.00 | 1516.00 | 1399.00 | 468.00 | 9.00 | 436.00 | | 5657.00 | | 197.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N5 | 1628.00 | 951.00 | 1535.00 | 114.00 | 1333.00 | 666.00 | 17.00 | | 6244.00 | | 2711.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | N6 | 485.00 | 2530.00 | 1706.00 | 1029.00 | 1207.00 | 322.00 | 1872.00 | | 9151.00 | | 2351.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | | | | | | | | | | | | | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ | S | 7026.00 | 10664.00 | 7205.00 | 5483.00 | 5460.00 | 3533.00 | 6800.00 | | 46171.00 | | 0.00 | +-------+---------+----------+---------+---------+---------+---------+---------+-----+----------+-----+----------+ To see the net effect of money transfers from node to node, we first form a skew-symmetric matrix(each entry is the negative of it's mirror reflection through the main diagonal) from TRM given by SKEW = TRM - TRM.transpose(). +--------+------+------+-------+-------+-------+-------+-------+ | SKEW | N0 | N1 | N2 | N3 | N4 | N5 | N6 | +========+======+======+=======+=======+=======+=======+=======+ | N0 | 0 | 418 | -817 | -336 | 143 | -1133 | 866 | +--------+------+------+-------+-------+-------+-------+-------+ | N1 | -418 | 0 | -1160 | -1256 | -522 | -584 | -141 | +--------+------+------+-------+-------+-------+-------+-------+ | N2 | 817 | 1160 | 0 | -593 | -844 | -953 | -1513 | +--------+------+------+-------+-------+-------+-------+-------+ | N3 | 336 | 1256 | 593 | 0 | -1069 | 978 | -487 | +--------+------+------+-------+-------+-------+-------+-------+ | N4 | -143 | 522 | 844 | 1069 | 0 | -1324 | -771 | +--------+------+------+-------+-------+-------+-------+-------+ | N5 | 1133 | 584 | 953 | -978 | 1324 | 0 | -305 | +--------+------+------+-------+-------+-------+-------+-------+ | N6 | -866 | 141 | 1513 | 487 | 771 | 305 | 0 | +--------+------+------+-------+-------+-------+-------+-------+ Next, we set the negative entries of SKEW to 0, and keep the rest to obtain the matrix, NET. For instance if you send me $10, and I send you $5, then the net effect is you sent me $5, and I sent you $0, we're just expressing everything in non-negative numbers. We could apply x*H(x) where H is the Heaviside step function to each entry of SKEW to obtain NET. Note that the R-S column vector for NET has not changed from that of TRM. We also note that for any square matrix, the sum over the R vector equals the sum over the S vector which is just the sum of all entries in the matrix, and therefore the R-S vector always sums to 0. +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | NET | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +=======+=========+=========+=========+=========+=========+=========+========+=====+==========+=====+==========+ | N0 | 0.00 | 418.00 | 0.00 | 0.00 | 143.00 | 0.00 | 866.00 | | 1427.00 | | -859.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 0.00 | | -4081.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N2 | 817.00 | 1160.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 1977.00 | | -1926.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N3 | 336.00 | 1256.00 | 593.00 | 0.00 | 0.00 | 978.00 | 0.00 | | 3163.00 | | 1607.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N4 | 0.00 | 522.00 | 844.00 | 1069.00 | 0.00 | 0.00 | 0.00 | | 2435.00 | | 197.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N5 | 1133.00 | 584.00 | 953.00 | 0.00 | 1324.00 | 0.00 | 0.00 | | 3994.00 | | 2711.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | N6 | 0.00 | 141.00 | 1513.00 | 487.00 | 771.00 | 305.00 | 0.00 | | 3217.00 | | 2351.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | | | | | | | | | | | | | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ | S | 2286.00 | 4081.00 | 3903.00 | 1556.00 | 2238.00 | 1283.00 | 866.00 | | 16213.00 | | 0.00 | +-------+---------+---------+---------+---------+---------+---------+--------+-----+----------+-----+----------+ Closer analysis using some graph theory and linear algebra allows us to express the net effect of all transactions in TRM efficiently as the projection of TRM into the cocycle space of the complete directed multigraph on N vertices and 2*N*(N-1) edges. The dimension of the cocycle space is N-1. Projecting TRM thusly we obtain the skew-symmetric matrix, LMD. I called it this for no good reason having inherited it from a Visual Basic version I wrote a number of years ago. Note here also, the R-S vector remains unchanged from TRM. As a further note the numbers here are rounded to two decimal places whereas they have many more significant digits so the numbers may not add up exactly, if checking. +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | LMD | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +=======+=========+=========+=========+=========+=========+==========+==========+=====+==========+=====+==========+ | N0 | 0.00 | 230.14 | 76.21 | -176.14 | -75.43 | -255.00 | -229.29 | | -429.50 | | -859.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N1 | -230.14 | 0.00 | -153.93 | -406.29 | -305.57 | -485.14 | -459.43 | | -2040.50 | | -4081.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N2 | -76.21 | 153.93 | 0.00 | -252.36 | -151.64 | -331.21 | -305.50 | | -963.00 | | -1926.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N3 | 176.14 | 406.29 | 252.36 | 0.00 | 100.71 | -78.86 | -53.14 | | 803.50 | | 1607.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N4 | 75.43 | 305.57 | 151.64 | -100.71 | 0.00 | -179.57 | -153.86 | | 98.50 | | 197.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N5 | 255.00 | 485.14 | 331.21 | 78.86 | 179.57 | 0.00 | 25.71 | | 1355.50 | | 2711.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | N6 | 229.29 | 459.43 | 305.50 | 53.14 | 153.86 | -25.71 | 0.00 | | 1175.50 | | 2351.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | | | | | | | | | | | | | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ | S | 429.50 | 2040.50 | 963.00 | -803.50 | -98.50 | -1355.50 | -1175.50 | | 0.00 | | 0.00 | +-------+---------+---------+---------+---------+---------+----------+----------+-----+----------+-----+----------+ To express the transactions in LMD as non-negative numbers only form SKEW = LMD - LMD.transpose, and set the negatives to 0. We obtain the matrix NETL. Again R-S is unchanged. +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | NETL | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +========+=========+=========+=========+========+========+======+=======+=====+=========+=====+==========+ | N0 | 0.00 | 460.29 | 152.43 | 0.00 | 0.00 | 0.00 | 0.00 | | 612.71 | | -859.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 0.00 | | -4081.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N2 | 0.00 | 307.86 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | | 307.86 | | -1926.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N3 | 352.29 | 812.57 | 504.71 | 0.00 | 201.43 | 0.00 | 0.00 | | 1871.00 | | 1607.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N4 | 150.86 | 611.14 | 303.29 | 0.00 | 0.00 | 0.00 | 0.00 | | 1065.29 | | 197.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N5 | 510.00 | 970.29 | 662.43 | 157.71 | 359.14 | 0.00 | 51.43 | | 2711.00 | | 2711.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | N6 | 458.57 | 918.86 | 611.00 | 106.29 | 307.71 | 0.00 | 0.00 | | 2402.43 | | 2351.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | | | | | | | | | | | | | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ | S | 1471.71 | 4081.00 | 2233.86 | 264.00 | 868.29 | 0.00 | 51.43 | | 8970.29 | | 0.00 | +--------+---------+---------+---------+--------+--------+------+-------+-----+---------+-----+----------+ Any matrix where the R-S vector is the 0-vector(not merely summing to 0, but every entry is 0) lies in the cycle space of the complete directed graph. We can project TRM into the cycle space to obtain a matrix, CYC. This can be obtained by simply subtracting LMD from TRM, thus CYC = TRM - LMD. CYC can be broken down further into symmetric and skew-symmetric matrices, CYCS and CYCK respectively where CYCS = (CYC + CYC.transpose())/2 and CYCK = (CYC - CYC.transpose())/2. We have TRM = LMD + CYCS + CYCK. +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | CYCS | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +========+=========+=========+=========+=========+=========+=========+=========+=====+==========+=====+=======+ | N0 | 1606.00 | 668.00 | 751.50 | 932.00 | 659.50 | 1061.50 | 918.00 | | 6596.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N1 | 668.00 | 1154.00 | 1258.00 | 1445.00 | 980.00 | 659.00 | 2459.50 | | 8623.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N2 | 751.50 | 1258.00 | 539.00 | 591.50 | 1094.00 | 1058.50 | 949.50 | | 6242.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N3 | 932.00 | 1445.00 | 591.50 | 1065.00 | 864.50 | 603.00 | 785.50 | | 6286.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N4 | 659.50 | 980.00 | 1094.00 | 864.50 | 468.00 | 671.00 | 821.50 | | 5558.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N5 | 1061.50 | 659.00 | 1058.50 | 603.00 | 671.00 | 666.00 | 169.50 | | 4888.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | N6 | 918.00 | 2459.50 | 949.50 | 785.50 | 821.50 | 169.50 | 1872.00 | | 7975.50 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | | | | | | | | | | | | | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ | S | 6596.50 | 8623.50 | 6242.00 | 6286.50 | 5558.50 | 4888.50 | 7975.50 | | 46171.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+----------+-----+-------+ The space of all N×N square matrices' dimension is N^2. The cocycle space is dimension N-1, so that the cycle space is dimension N^2 - N + 1. The cycle space is further broken down into the symmetric cycles and the skew-symmetric cycles. Every symmetric matrix is a cycle, and the dimension of the symmetric matrices is N*(N+1)/2. The dimension of the skew-symmetric matrices is N*(N-1)/2 leaving the dimension of the skew-symmetric cycles (N-1)*(N-2)/2. +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | CYCK | N0 | N1 | N2 | N3 | N4 | N5 | N6 | | R | | R-S | +========+=========+=========+=========+=========+=========+=========+=========+=====+=======+=====+=======+ | N0 | 0.00 | -21.14 | -484.71 | 8.14 | 146.93 | -311.50 | 662.29 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N1 | 21.14 | 0.00 | -426.07 | -221.71 | 44.57 | 193.14 | 388.93 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N2 | 484.71 | 426.07 | 0.00 | -44.14 | -270.36 | -145.29 | -451.00 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N3 | -8.14 | 221.71 | 44.14 | 0.00 | -635.21 | 567.86 | -190.36 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N4 | -146.93 | -44.57 | 270.36 | 635.21 | 0.00 | -482.43 | -231.64 | | 0.00 | | 0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N5 | 311.50 | -193.14 | 145.29 | -567.86 | 482.43 | 0.00 | -178.21 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | N6 | -662.29 | -388.93 | 451.00 | 190.36 | 231.64 | 178.21 | 0.00 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | | | | | | | | | | | | | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ | S | 0.00 | 0.00 | 0.00 | 0.00 | -0.00 | 0.00 | 0.00 | | -0.00 | | -0.00 | +--------+---------+---------+---------+---------+---------+---------+---------+-----+-------+-----+-------+ 
CPR = PolynomialRing(QQ,nn,'x',order=myBP.parent().term_order())
CPR.inject_variables()
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
newpoly
Defining x0, x1, x2 x0*x1*x2 + x0*x1 + x0*x2 + x0 + x1*x2 + x1 + x2 + 1 
from bitarray import bitarray
b = bitarray()
b.append(1)
b
b.extend([1,0])
b
len(b)
rows = 10
cols = 22
B=[]
B.append(rows)
B.append(cols)
B.append(bitarray(rows*cols))
B.append(bitarray(rows*cols))
B
B[2][B[1]*5+19]
B[3][B[1]*5+19]=1
B[3][B[1]*5+19]
B[3]
bitarray('1') bitarray('110') 3 [10, 22, bitarray('0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000'), bitarray('0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000')] 0 1 bitarray('0000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000') 
import csv
import sys
from bitarray import bitarray
from tabulate import tabulate
load('RVal.py')
load('SVal.py')
load('Xchange.py')
load('XchangeA.py')
load('MinShip.py')
load('printPHI.py')
rows = 10
cols = 22
PHI=[]
PHI.append(rows)
PHI.append(cols)
PHI.append(bitarray(rows*cols))
PHI.append(bitarray(rows*cols))
f = open('Book1.csv', 'rt')
try:
    reader = csv.reader(f)
    ii = 0
    for row in reader:
        jj = 0
        for x in row:
            SVal(PHI,ii,jj,int(x))
            jj = jj + 1
        ii = ii + 1
finally:
    f.close()
PHIorig = deepcopy(PHI)
print(' ')
#supply and demand vector
Rmat = matrix(ZZ,rows,rows,0)
for ii in range(rows):
    for jj in range(rows):
        Rmat[ii,jj] = round(normalvariate(0,1000))
SDvect = [0]*rows
for ii in range(rows):
    for jj in range(rows):
        SDvect[ii] = SDvect[ii] + Rmat[ii, jj] - Rmat[jj, ii]
mysum = 0
for ii in range(rows):
    mysum = mysum + SDvect[ii]
SDvect
mysum
#print(' ')
cost = [1]*PHI[1]
corner = [0]
#SDvect = [-6,-14,-18,-6,-46,-14,24,-9,30,59]
colH = [1..PHI[1]]
rowH = [PHI[1]+1..PHI[1]+PHI[0]]
costorig = cost[:]
SDvectorig = SDvect[:]
colHorig = colH[:]
rowHorig = rowH[:]
cornerorig = corner[:]

print(' ')
PHI
print(' ')
printPHI(PHI,cost,SDvect,corner,rowH,colH,'PHI')

#pivotrow = 0
#pivotcol = 0
#XchangeA(PHI, pivotrow, pivotcol, cost, SDvect, corner, rowH, colH)
print(' ')
MS = MinShip(PHI, cost, SDvect, corner, rowH, colH)

print('Simplex method complete after '+str(MS[0])+' Xchange operations')
PHI
print(' ')
printPHI(PHI,cost,SDvect,corner,rowH,colH,'PHI')
print(' ')
print('Shipping vector =')
[round(x,2) for x in MS[1]]

TransMat = matrix(RR,rows,rows,0.0)
for jj in range(cols):
    for ii in range(rows):
        if RVal(PHIorig,ii,jj) > 0: aa = ii
        if RVal(PHIorig,ii,jj) < 0: bb = ii
    TransMat[aa,bb] = MS[1][jj]

print(' ')
TransMatdata = []
for ii in range(rows):
    row = []
    for jj in range(rows+1):
        row.append('7')
    TransMatdata.append(row)
for ii in range(rows):
    TransMatdata[ii][0] = str(rowHorig[ii])
    for jj in range(rows):
        TransMatdata[ii][jj+1] = round(TransMat[ii][jj],2)

hd = [str(x) for x in rowHorig]
cola = ['right']*rows
print(tabulate(TransMatdata, headers=hd, tablefmt="grid", colalign=cola))
RSvect = [0]*rows
for ii in range(rows):
    for jj in range(rows):
        RSvect[ii] = RSvect[ii] + TransMat[ii, jj] - TransMat[jj, ii]
mysum = 0
for ii in range(rows):
    mysum = mysum + RSvect[ii]
print(' ')
print('Receive minus Send vector')
[round(x,2) for x in RSvect]
mysum
#printPHI(PHIorig,costorig,SDvectorig,cornerorig,rowHorig,colHorig,'PHI')
[-8628, 1767, 3650, -5510, -73, -3795, -929, 7069, 1690, 4759] 0 [10, 22, bitarray('0000000000000000000000100000000000000000000001010000000000000000000010101000000000000000000001010100000000000000000000100100000000000000000000101000100000000000000000010001000000000000000000100011000000000000000001000011'), bitarray('1110000000000000000000100111000000000000000001010011100000000000000010101001100000000000000001010101111100000000000000100100001100000000000000101000101000000000000000010001011000000000000000100011010000000000000001000011')] +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | PHI | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | 17 | 18 | 19 | 20 | 21 | 22 | SD | +=======+=====+=====+=====+=====+=====+=====+=====+=====+=====+======+======+======+======+======+======+======+======+======+======+======+======+======+=======+ | 23 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -8628 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 24 | -1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1767 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 25 | 0 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3650 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 26 | 0 | 0 | -1 | 0 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -5510 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 27 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | -1 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | -73 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 28 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | -3795 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 29 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | 0 | 0 | -1 | 0 | 1 | 0 | 0 | 0 | -929 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 30 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | -1 | 0 | 1 | 1 | 0 | 7069 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 31 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | -1 | -1 | 0 | 1 | 1690 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | 32 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | -1 | -1 | 4759 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ | cost | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 0 | +-------+-----+-----+-----+-----+-----+-----+-----+-----+-----+------+------+------+------+------+------+------+------+------+------+------+------+------+-------+ Simplex method complete after 21 Xchange operations [10, 22, bitarray('0000000000000000000000000000010000000000000000011000001100000000000100000000000000000000000010000011100010000000000000000000000010000000001001000010100000000000010000010101110000000000000000000000111110010011000111000000'), bitarray('0000000100000100011000000000110100000010000011011000001100000000000101100000100000000000111010000011100010000000001100000110000010000000001101000010100110000000010000010101110111100100110001110000111110010011000111000000')] +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | PHI | 23 | 24 | 25 | 4 | 5 | 26 | 12 | 17 | 27 | 28 | 6 | 3 | 7 | 29 | 30 | 31 | 9 | 13 | 11 | 20 | 21 | 22 | SD | +=======+======+======+======+=====+=====+======+======+======+======+======+=====+=====+=====+======+======+======+=====+======+======+======+======+======+=======+ | 19 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 929 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 18 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | -1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 3795 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 2 | 1 | 1 | 0 | -1 | -1 | 0 | 0 | 0 | 0 | 0 | -1 | -1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 6861 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 1 | 0 | -1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1767 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 8 | 1 | 1 | 1 | 0 | -1 | 0 | 0 | 0 | 0 | 0 | -1 | -1 | -1 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | 0 | 3211 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 10 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 5510 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 14 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 1 | 0 | -1 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | 0 | 1 | 1 | 0 | 3274 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 15 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | 0 | 0 | 0 | 0 | -1 | 0 | -1 | 0 | -1 | -1 | -1 | 0 | 1 | 761 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 16 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 1 | 1 | 0 | 0 | 0 | 0 | -1 | -1 | 4759 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | 32 | -1 | -1 | -1 | 0 | 0 | -1 | 0 | 0 | -1 | -1 | 0 | 0 | 0 | -1 | -1 | -1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ | cost | 3 | 2 | 2 | 1 | 1 | 2 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 1 | 0 | 1 | 1 | 1 | 30867 | +-------+------+------+------+-----+-----+------+------+------+------+------+-----+-----+-----+------+------+------+-----+------+------+------+------+------+-------+ Shipping vector = [1767.0, 6861.0, 0.0, 0.0, 0.0, 0.0, 0.0, 3211.0, 0.0, 5510.0, 0.0, 0.0, 0.0, 3274.0, 761.0, 4759.0, 0.0, 3795.0, 929.0, 0.0, 0.0, 0.0] +----+------+------+------+------+------+------+------+------+------+------+ | | 23 | 24 | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | +====+======+======+======+======+======+======+======+======+======+======+ | 23 | 0 | 1767 | 6861 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 24 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 25 | 0 | 0 | 0 | 0 | 3211 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 26 | 0 | 0 | 0 | 0 | 5510 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 27 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3274 | 761 | 4759 | +----+------+------+------+------+------+------+------+------+------+------+ | 28 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 3795 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 29 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 929 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 30 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 31 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ | 32 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | +----+------+------+------+------+------+------+------+------+------+------+ Receive minus Send vector [8628.0, -1767.0, -3650.0, 5510.0, 73.0, 3795.0, 929.0, -7069.0, -1690.0, -4759.0] 0.000000000000000