Double-star

1483 days ago by hogben

# This worksheet is generates data for # "Using Markov chains to determine expected propagation time for probabilistic zero forcing" # by Yu Chan, Emelie Curl, Jesse Geneson, Leslie Hogben, Kevin Liu, Issac Odegard, and Michael Ross # This worksheet is written by Leslie Hogben # It contains the code for computing the expected propagation time of # double stars and subdivided double stars. # The cell ept computes the expected propagation time of any graph if given its Markov matrix (using Theorem 2.2). # This worksheet is written by Leslie Hogben 
       
### The function ept computes expected propagation time given a Markov matrix M by applying Theorem 2.2 ### 
       
def ept(M): nstates=M.dimensions()[0] # matrix dimension jay=ones_matrix(nstates) eye=identity_matrix(nstates) # x.outer_product(y) does x*y^t for column vectors x,y Mi=(M-jay[0].outer_product(eye[nstates-1])-eye).inverse() ept=Mi[0,nstates-1]+1 return ept 
       
### The next cell defines the Markov matrix of the subdivided double star DS(t,t)_uu' # where the new vertex w is in the subdivded edge uu' # t is the # of leaves at each of the centers u and u' which are adjacent in DS(t,t) # the order is 2t+3 
       
def MDSsubd(t): # starting with a degee 2 vertex blue; u1 and u2 are the two centers d=t+1 # degree of HD vertices u1 and u2 # states are (0),(1) and (2,k1,k2) for k1,k2 = 0,...,t; (0) means neither u1 nor u2 blue # (1) means one of u2,u1 blue and other white; for book-keeping u2 blue # states (2,k1,k2) means u1,u2 blue and k1,k2 are the # of blue leaves m=zero_matrix(RR,(t+1)^2+2) # m is the matrix; use QQ for exact, RR for decimal # matrix row 0 is neither of u1,u2 blue = state (0) # matrix row 1 is 1 of u2,u1 blue (say u2) = state (1) # for both u1 and u2 blue, use lex order for states (2,k1,k2), i.e., # matrix row 2+k1(t+1)+k2 is state (2,k1,k2) # # forcing with both u1 and u2 white # row 0 = state (0): the degree 2 blue vertex can force 0, 1, or 2 of u1 and u2 m[0,0]=1/4 m[0,1]=1/2 m[0,2]=1/4 # forcing with u2 blue and u1 white # u2 blue implies u1 will become blue in next round, forced by the degree 2 vertex w pfu=2/d # probability of u2 forcing any one specific leaf neighbor pnfu=1-pfu # probability of u not forcing any one specific leaf neighbor for r in [0..t]: # r is the number of new blue leaves m[1,2+r]=binomial(t,r)*pfu^r*pnfu^(t-r) # starting with u1 and u2 both blue # state (2,k1,k2) is matrix index 2+k1(t+1)+k2 for k1,k2=0,...,t prm=zero_matrix(QQ,t+1) # prm is the leaf probability matrix; use QQ for exact, RR for decimal # prm[i,j] is the probability of going from i blue leaves to j blue leaf neighbprs of u2 (or u1) for i in [0..t-2]: # construct prm pfu=(i+2)/d # probability center forcing one specific leaf when i leaf neighbors are blue pnfu=1-pfu # probability center not forcing any one specific leaf neighbor when i are blue for j in [i..t]: x=t-i # x is # white vertices r=j-i # r is the number of new blue leaves prm[i,j]=binomial(x,r)*pfu^r*pnfu^(x-r) # prob r additional blue u-leaves prm[t-1,t]=1 prm[t,t]=1 # construct transition matrix m for k1 in [0..t]: for k2 in [0..t]: for r1 in [0..t-k1]: for r2 in [0..t-k2]: m[2+k1*(t+1)+k2,2+(k1+r1)*(t+1)+k2+r2]=prm[k1,k1+r1]*prm[k2,k2+r2] return m 
       
### The next cell defines the Markov matrix of the double star DS(t,t) # where t is the # of leaves at each of the adjacent centers u1 and u2 # starting with center u2 blue # the order of DS(t,t) is 2t+2 ### 
       
def MDSc(t): # starting with center u2 blue; u1 and u2 are the two centers d=t+1 # degree of centers vertices u1 and u2 # states are (1,k2) and (2,k1,k2) for k1,k2 = 0,...,t # (1,k2) means u2 blue and u1 white, k2 is # blue leaves of u2. # states (2,k1,k2) means u1,u2 blue and k1,k2 are the # of blue leaves of each m=zero_matrix(RR,1+t+(t+1)^2) # m is the matrix; use QQ for exact, RR for decimal # # matrix row 0...t: row k2 is u2 blue and k2 blue leaves = state (1,k2) # forcing with u2 blue and u1 white # disinguish whether u1 gets forced for k2 in [0..t-1]: # k2 is the number of new blue leaves starting this round pfu=(k2+1)/d # probability of u2 forcing any one specific leaf neighbor pnfu=1-pfu # probability of u not forcing any one specific leaf neighbor # when u1 is not forced (k2 blue => index k2) for r in [0..t-k2]: # r is the number of new blue leaves m[k2,k2+r]=binomial(t-k2,r)*pfu^r*pnfu^(t-k2-r+1) # when u1 is forced (k2 blue => index k2+t+1) for r in [0..t-k2]: # r is the number of new blue leaves m[k2,1+t+k2+r]=binomial(t-k2,r)*pfu^(r+1)*pnfu^(t-k2-r) m[t,2*t+1]=1 # starting with u1 and u2 both blue # for both u1 and u2 blue, use lex order for states (2,k1,k2), i.e., # matrix row 1+t+k1(t+1)+k2 is state (2,k1,k2) prm=zero_matrix(QQ,t+1) # prm is the leaf probability matrix; use QQ for exact, RR for decimal # prm[i,j] is the probability of going from i blue leaves to j blue leaf neighbprs of u2 (or u1) for i in [0..t-2]: # construct prm pfu=(i+2)/d # probability center forcing one specific leaf when i leaf neighbors are blue pnfu=1-pfu # probability center not forcing any one specific leaf neighbor when i are blue for j in [i..t]: x=t-i # x is # white vertices r=j-i # r is the number of new blue leaves prm[i,j]=binomial(x,r)*pfu^r*pnfu^(x-r) # prob r additional blue u-leaves prm[t-1,t]=1 prm[t,t]=1 # construct transition matrix m for k1 in [0..t]: for k2 in [0..t]: for r1 in [0..t-k1]: for r2 in [0..t-k2]: m[1+t+k1*(t+1)+k2,1+t+(k1+r1)*(t+1)+k2+r2]=prm[k1,k1+r1]*prm[k2,k2+r2] return m 
       
### The next cell defines the Markov matrix of the double star DS(t,t) # where t is the # of leaves at each of the adjacent centers u1 and u2 # starting with a blue leaf neighbor of center u2 ### 
       
def MDSl(t): # starting with one leaf neighbor u2 blue; u1 and u2 are the two centers d=t+1 # degree of HD vertices u1 and u2 # states are (0) = no centers blue, (1,k2) and (2,k1,k2) for k1=1,...,t and k2 = 0,...,t # (1,k2) means u2 blue and u1 white, k2 is # blue leaves of u2. # states (2,k1,k2) means u1,u2 blue and k1,k2 are the # of blue leaves of each m=zero_matrix(QQ,1+t+(t+1)*t) # m is the matrix; use QQ for exact, RR for decimal # # matrix row 0: bkue leaf forces to state (1,1) = index 1 m[0,1]=1 # matrix row 1...t: row k2 is u2 blue and k2 blue leaves = state (1,k2) # forcing with u2 blue and u1 white # disinguish whether u1 gets forced for k2 in [1..t-1]: # k2 is the number of new blue leaves starting this round pfu=(k2+1)/d # probability of u2 forcing any one specific leaf neighbor pnfu=1-pfu # probability of u not forcing any one specific leaf neighbor # when u2 is not forced (k2 blue => index k2) for r in [0..t-k2]: # r is the number of new blue leaves m[k2,k2+r]=binomial(t-k2,r)*pfu^r*pnfu^(t-k2-r+1) # when u2 is forced (k2 blue => index k2+t) for r in [0..t-k2]: # r is the number of new blue leaves m[k2,1+t+k2+r]=binomial(t-k2,r)*pfu^(r+1)*pnfu^(t-k2-r) m[t,2*t+1]=1 # starting with u1 and u2 both blue # for both u1 and u2 blue, use lex order for states (2,k1,k2), i.e., # matrix row t+k1(t+1)+k2 is state (2,k1,k2) prm=zero_matrix(QQ,t+1) # prm is the leaf probability matrix; use QQ for exact, RR for decimal # prm[i,j] is the probability of going from i blue leaves to j blue leaf neighbprs of u2 (or u1) for i in [0..t-2]: # construct prm pfu=(i+2)/d # probability center forcing one specific leaf when i leaf neighbors are blue pnfu=1-pfu # probability center not forcing any one specific leaf neighbor when i are blue for j in [i..t]: x=t-i # x is # white vertices r=j-i # r is the number of new blue leaves prm[i,j]=binomial(x,r)*pfu^r*pnfu^(x-r) # prob r additional blue u-leaves prm[t-1,t]=1 prm[t,t]=1 # construct transition matrix m for k1 in [0..t]: for k2 in [1..t]: for r1 in [0..t-k1]: for r2 in [0..t-k2]: m[t+k1*t+k2,t+(k1+r1)*t+k2+r2]=prm[k1,k1+r1]*prm[k2,k2+r2] return m 
       
def MDS1(t): # starting with one leaf neighbor u2 blue; u1 and u2 are the two centers d=t+1 # degree of HD vertices u1 and u2 # states are (0) = no centers blue, (1,k2) and (2,k1,k2) for k1=1,...,t and k2 = 0,...,t # (1,k2) means u2 blue and u1 white, k2 is # blue leaves of u2. # states (2,k1,k2) means u1,u2 blue and k1,k2 are the # of blue leaves of each m=zero_matrix(RR,1+t+(t+1)*t) # m is the matrix; use QQ for exact, RR for decimal # # matrix row 0: bkue leaf forces to state (1,1) = index 1 m[0,1]=1 # matrix row 1...t: row k2 is u2 blue and k2 blue leaves = state (1,k2) # forcing with u2 blue and u1 white # disinguish whether u1 gets forced for k2 in [1..t-1]: # k2 is the number of new blue leaves starting this round pfu=(k2+1)/d # probability of u2 forcing any one specific leaf neighbor pnfu=1-pfu # probability of u not forcing any one specific leaf neighbor # when u2 is not forced (k2 blue => index k2) for r in [0..t-k2]: # r is the number of new blue leaves m[k2,k2+r]=binomial(t-k2,r)*pfu^r*pnfu^(t-k2-r+1) # when u2 is forced (k2 blue => index k2+t) for r in [0..t-k2]: # r is the number of new blue leaves m[k2,t+k2+r]=binomial(t-k2,r)*pfu^(r+1)*pnfu^(t-k2-r) m[t,2*t]=1 # starting with u1 and u2 both blue # for both u1 and u2 blue, use lex order for states (2,k1,k2), i.e., # matrix row t+k1(t+1)+k2 is state (2,k1,k2) prm=zero_matrix(QQ,t+1) # prm is the leaf probability matrix; use QQ for exact, RR for decimal # prm[i,j] is the probability of going from i blue leaves to j blue leaf neighbprs of u2 (or u1) for i in [0..t-2]: # construct prm pfu=(i+2)/d # probability center forcing one specific leaf when i leaf neighbors are blue pnfu=1-pfu # probability center not forcing any one specific leaf neighbor when i are blue for j in [i..t]: x=t-i # x is # white vertices r=j-i # r is the number of new blue leaves prm[i,j]=binomial(x,r)*pfu^r*pnfu^(x-r) # prob r additional blue u-leaves prm[t-1,t]=1 prm[t,t]=1 # construct transition matrix m for k1 in [0..t]: for k2 in [1..t]: for r1 in [0..t-k1]: for r2 in [0..t-k2]: m[t+k1*t+k2,t+(k1+r1)*t+k2+r2]=prm[k1,k1+r1]*prm[k2,k2+r2] return m 
       
for t in [2..12]: print t, ept(MDSc(t)).n(digits=6), ept(MDS1(t)).n(digits=6), ept(MDSsubd(t)).n(digits=6) 
       
2 3.79605 4.00000 3.75231
3 4.71003 4.85852 4.36807
4 5.42428 5.53190 4.84690
5 6.02081 6.09666 5.24553
6 6.53274 6.58471 5.58244
7 6.97868 7.01109 5.87223
8 7.37310 7.38881 6.12644
9 7.72679 7.72817 6.35263
10 8.04739 8.03639 6.55591
11 8.34040 8.31854 6.74013
12 8.61003 8.57850 6.90842
2 3.79605 4.00000 3.75231
3 4.71003 4.85852 4.36807
4 5.42428 5.53190 4.84690
5 6.02081 6.09666 5.24553
6 6.53274 6.58471 5.58244
7 6.97868 7.01109 5.87223
8 7.37310 7.38881 6.12644
9 7.72679 7.72817 6.35263
10 8.04739 8.03639 6.55591
11 8.34040 8.31854 6.74013
12 8.61003 8.57850 6.90842
for t in [13..15]: print t, ept(MDSc(t)).n(digits=6), ept(MDS1(t)).n(digits=6), ept(MDSsubd(t)).n(digits=6) 
       
13 8.85959 8.81937 7.06324
14 9.09181 9.04372 7.20656
15 9.30892 9.25368 7.33989
13 8.85959 8.81937 7.06324
14 9.09181 9.04372 7.20656
15 9.30892 9.25368 7.33989