CoCalc -- ProbM.py

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