1
def ProbM(myBP): # myBP is a boolean polynomial being converted to probabilty expression
2
    # myBP.parent().term_order() must be 'lp'
3
    if type(myBP) == sage.rings.integer.Integer:
4
        if myBP%2 == 0:
5
            return QQ(0)
6
        else:
7
            return QQ(1)
8
    mystr = str(myBP.variables()) # converting all uppers to lower and lowers to uppercase I SEE A BUG HERE!!!!!! mixed upper and lower var names sort differently
9
    NewGens=""
10
    for ii in range(len(mystr)):
11
        if mystr[ii] == mystr[ii].upper():
12
            NewGens = NewGens + mystr[ii].lower()
13
        else:
14
            NewGens = NewGens + mystr[ii].upper() # if boolean variable XcGGt1u then, prob var is xCggT1U etc. Mixed case generators are a problem - they sort differently
15
    NewGens = NewGens[1:len(NewGens) - 1] # generators for probabilty polynomial assuming boolean generators are independent
16
    if len(myBP.variables()) == 1:
17
        NewGens = NewGens.replace(',','')
18
    NewGens.replace(' ','')
19
    CPR = PolynomialRing(QQ,len(myBP.variables()),NewGens,order=myBP.parent().term_order()) #myBP.parent().term_order()
20
    CPR.inject_variables()
21
    nn = myBP.nvariables() # number of generators in boolean polynomial
22
    atoms = 2**nn # number of atoms in corresponding boolean algebra
23
    # should abort this with error, if too many variables in polynomial
24
    mypoly=myBP.parent(1)
25
    for ii in range(0,len(myBP.parent().gens())):
26
        mypoly = mypoly*(myBP.parent(1) + myBP.parent().gen(ii))
27
    Mvect = [0]*atoms
28
    for y in myBP.monomials():
29
        for jj in range(len(mypoly.monomials())):
30
            if y==mypoly.monomials()[jj]:
31
                Mvect[jj]=1
32
    #Mvect = [0]*atoms # This loop only works when CPR term_order is 'invlex'
33
    #for jj in range(len(myBP.monomials())):
34
    #    mm = 0
35
    #    for kk in range(len(myBP.monomials()[jj].variables())):
36
    #        mm = mm + 2**myBP.variables().index(myBP.monomials()[jj].variables()[kk])
37
    #    Mvect[mm] = 1
38
    #Mvect.reverse()
39
    # Mvect order for instance is (X111,X110,X101,X100,X011,X010,X001,X000) has two nonequiv intepretations
40
    # as monomials or atoms example was for three variables
41
    Mvect = vector(Mvect)
42
    MGinv = matrix(ZZ,[1])
43
    for ii in range(1,nn+1):
44
        MGinv = block_matrix(ZZ,2,2,[MGinv,-MGinv,0,MGinv],subdivide=False) # creating inverse of multigrade operator 'AND' matrix MG
45
        # modulo 2 it's an involution however...
46
        # MG and MGinv are upper triangular, MG is like Sierpinski triangle see  
47
        # https://commons.wikimedia.org/wiki/File:Multigrade_operator_AND.svg
48
    Avect = (MGinv*Mvect)%2
49
    Pvect = MGinv*Avect # convert from boolean ring (mons) to boolean algebra (atoms) then to probability polynomial
50
    # ((MGinv*Mvect)%2) is the atoms vector equivalent to the monnomials vector Mvect
51
    # Pvect are the coefficients of prob poly in CPR
52
    # CPR.term_order() = myBP.parent().term_order()
53
    newpoly = CPR(1)
54
    for ii in range(0,len(CPR.gens())):
55
         newpoly = newpoly*(CPR(1) + CPR.gen(ii)) # to create all possible monomials in CPR
56
    return Pvect*vector(newpoly.monomials()) # inner product of the two vectors to create prob poly
57