PropagationTime_Throttling

954 days ago by hogben

# Sage code and examples for variants of propagation time and throttling 
       
#################################### # This Sage worksheet is written by Leslie Hogben based on code written by numerous other people. # The focus of this worksheet is on propagation time and throttling. # Examples are from the book, # "Inverse Problems and Zero Forcing for Graphs" # by Leslie Hogben, Jehian C.H. Lin, and Bryan L. Shader # Please report errors to hogben@aimath.org and jephianlin@gmail.com # # This worksheet begins with examples, which are ready to read. # In order to run the examples for yourself, you first need to enter the function cells. # To do that, go to the bottom of this file (or search for FUNCTIONS-HERE) and enter each cell that begins # F-. #################################### 
       
#################################### # Section 10.2 # Propagation time of a set S, minimum propagation time of a set of k elements, # and propagation time of a graph are illustrated. # These are computed by functions ptz(G,S), ptk(G), and pt(G). # Examples illustrate possible effects of edge subdivision. #################################### 
       
# Example PT.2-1: # The center of W_7 is in every efficient standard zero forcing set. # A set of 3 consecutive degree 3 vertices is an efficient for W_8 that does not contain the center. 
       
# Wheel W_7 g=graphs.WheelGraph(7) show(g) print 'Z(G) =',Z(g) 
       
Z(G) = 3
Z(G) = 3
# propagation time of a set B = {0,1,2} ptz(g,[0,1,2]) 
       
2
2
ptz(g,[1,2,3]) 
       
3
3
# Propagation time of W_7. print 'pt(G) =', pt(g) 
       
pt(G) = 2
pt(G) = 2
# Wheel W_8 g=graphs.WheelGraph(8) show(g) z=Z(g) print 'Z(G)=',z print 'pt(G)=', pt(g) 
       
Z(G)= 3
pt(G)= 3
Z(G)= 3
pt(G)= 3
ptz(g,[1,2,3]) 
       
3
3
# Example PT.2-2: Book graph B_3. 
       
r=3 star=graphs.CompleteBipartiteGraph(1,r) p2=graphs.PathGraph(2) g=star.cartesian_product(p2) g.relabel() show(g) z=Z(g) print 'Z =', z, 'pt =',ptk(g,z) 
       
Z = 3 pt = 3
Z = 3 pt = 3
S=ZFsets(g,r) for s in S: print s, 'pt =',ptz(g,s) 
       
(0, 2, 4) pt = 3
(0, 2, 5) pt = 4
(0, 2, 6) pt = 3
(0, 2, 7) pt = 4
(0, 3, 4) pt = 4
(0, 3, 6) pt = 4
(0, 4, 6) pt = 3
(0, 4, 7) pt = 4
(0, 5, 6) pt = 4
(1, 2, 5) pt = 4
(1, 2, 7) pt = 4
(1, 3, 4) pt = 4
(1, 3, 5) pt = 3
(1, 3, 6) pt = 4
(1, 3, 7) pt = 3
(1, 4, 7) pt = 4
(1, 5, 6) pt = 4
(1, 5, 7) pt = 3
(2, 3, 4) pt = 3
(2, 3, 5) pt = 3
(2, 3, 6) pt = 3
(2, 3, 7) pt = 3
(2, 4, 5) pt = 3
(2, 6, 7) pt = 3
(3, 4, 5) pt = 3
(3, 6, 7) pt = 3
(4, 5, 6) pt = 3
(4, 5, 7) pt = 3
(4, 6, 7) pt = 3
(5, 6, 7) pt = 3
(0, 2, 4) pt = 3
(0, 2, 5) pt = 4
(0, 2, 6) pt = 3
(0, 2, 7) pt = 4
(0, 3, 4) pt = 4
(0, 3, 6) pt = 4
(0, 4, 6) pt = 3
(0, 4, 7) pt = 4
(0, 5, 6) pt = 4
(1, 2, 5) pt = 4
(1, 2, 7) pt = 4
(1, 3, 4) pt = 4
(1, 3, 5) pt = 3
(1, 3, 6) pt = 4
(1, 3, 7) pt = 3
(1, 4, 7) pt = 4
(1, 5, 6) pt = 4
(1, 5, 7) pt = 3
(2, 3, 4) pt = 3
(2, 3, 5) pt = 3
(2, 3, 6) pt = 3
(2, 3, 7) pt = 3
(2, 4, 5) pt = 3
(2, 6, 7) pt = 3
(3, 4, 5) pt = 3
(3, 6, 7) pt = 3
(4, 5, 6) pt = 3
(4, 5, 7) pt = 3
(4, 6, 7) pt = 3
(5, 6, 7) pt = 3
# Example PT.2-3: Edge subdivision can lower propagation time substantially. 
       
# construct the graph G(k) k=5 g=graphs.PathGraph(2*k) g.relabel([1..2*k]) p2=graphs.PathGraph(3) p2.relabel([2*k+1..2*k+3]) g=g.union(p2) g.add_edge(k,2*k+2) show(g) 
       
z=Z(g) print 'Z(G)=',z print 'pt(G)=', ptk(g,z) 
       
Z(G)= 2
pt(G)= 9
Z(G)= 2
pt(G)= 9
g.subdivide_edge(5,12,1) show(g) z=Z(g) print 'Z(G)=',z print 'pt(G)=', pt(g) 
       
Z(G)= 3
pt(G)= 6
Z(G)= 3
pt(G)= 6
# Example PT.2-4: Edge subdivision can raise propagation time substantially. 
       
# construct the graph H(k) k=7 g=graphs.PathGraph(k) g.relabel([1..k]) p2=graphs.PathGraph(k) p2.relabel([k+1..2*k]) g=g.union(p2) g.add_edges([(2,k+2),(3,k+3)]) show(g) 
       
z=Z(g) print 'Z(G)=',z print 'pt(G)=', ptk(g,z) 
       
Z(G)= 2
pt(G)= 6
Z(G)= 2
pt(G)= 6
g.subdivide_edge(3,10,1) show(g) 
       
z=Z(g) print 'Z(G)=',z print 'pt(G)=', ptk(g,z) 
       
Z(G)= 2
pt(G)= 10
Z(G)= 2
pt(G)= 10
# Example PT.2-5: All minimum zero forcing sets of H(k). 
       
# construct the graph H(k) k=7 h=graphs.PathGraph(k) h.relabel([1..k]) p2=graphs.PathGraph(k) p2.relabel([k+1..2*k]) h=h.union(p2) h.add_edges([(2,k+2),(3,k+3)]) show(h) z=Z(h) print 'Z(G)=',z 
       
Z(G)= 2
Z(G)= 2
zsets=ZFsets(h,z) zsets 
       
[(1, 7), (1, 8), (7, 14), (8, 14)]
[(1, 7), (1, 8), (7, 14), (8, 14)]
B2=zsets[0] B1=zsets[1] print 'B_1 =', B1, 'pt(H(7),B_1) =',ptz(h,B1) print 'B_2 =', B2, 'pt(H(7),B_2) =',ptz(h,B2) 
       
B_1 = (1, 8) pt(H(7),B_1) = 6
B_2 = (1, 7) pt(H(7),B_2) = 9
B_1 = (1, 8) pt(H(7),B_1) = 6
B_2 = (1, 7) pt(H(7),B_2) = 9
#################################### # Section 10.3 # Propagation time of a set S, minimum propagation time of a set of k elements, # and propagation time of a graph are illustrated. # These are computed by functions pt_plus(G,S), ptpk(G,k), and ptp(G). # Examples illustrate possible effect of contraction. ################################### 
       
# Example PT.3-1: The wheel does not have a minimum PSD forcing set with full forcing. 
       
# Wheel W_8 g=graphs.WheelGraph(8) show(g) z=Zplus(g) print 'Z(G)=',z print 'ptp(G)=', ptpk(g,z) 
       
Z(G)= 3
ptp(G)= 2
Z(G)= 3
ptp(G)= 2
pt_plus(g,[0,1,4]) 
       
2
2
pt_plus(g,[1,4,6]) 
       
-1
-1
pt_plus(g,[1,2,3]) 
       
3
3
# pt_plus = =1 means not a PSD forcing set 
       
# Example PT.3-2: Contracting an edge can lower PSD propagation time substantially. 
       
# construct the graph L(k) k=8 x=graphs.CycleGraph(2*k) x.add_edges([(0,k),(0,2),(0,2*k-2)]) x.relabel([1..2*k]) show(x) z=Zplus(x) print 'Zplus(L(k))=',z print 'ptp(L(k))=', ptp(x) 
       
Zplus(L(k))= 2
ptp(L(k))= 7
Zplus(L(k))= 2
ptp(L(k))= 7
x.contract_edge(1,k+1) show(x) z=Zplus(x) print 'Zplus(L(k)_e)=',z print 'ptp(L(k)_e)=', ptp(x) 
       
Zplus(L(k)_e)= 3
ptp(L(k)_e)= 2
Zplus(L(k)_e)= 3
ptp(L(k)_e)= 2
#################################### # Section 10.4 # Skew propagation time of a set B in a graph G and skew propagation time of G are illustrated. # These are computed by functions prop_time_unlooped(G,B) and pts(G). # Use of the simultaneous leaf stripping algorithm is illustrated. # Examples are presented illustrating effect of deletion of an independent twin and # showing the skew propagation interval is not full. ################################### 
       
# Example PT.4-1: pt_-(P_{4k+1})= k < 2k = ecc(P_{4k+1}) 
       
k=5 g=graphs.PathGraph(4*k+1) zm=Zskew(g) print 'Z_-(G) =', zm, 'pt_-(G) =',pts(g),'ecc(G) =',g.eccentricity(2*k) 
       
Z_-(G) = 1 pt_-(G) = 5 ecc(G) = 10
Z_-(G) = 1 pt_-(G) = 5 ecc(G) = 10
# Example PT.4-2: Simultaneous leaf stripping. 
       
g=Graph({2:[1,4,5,6],4:[3,5,6],5:[6,18],7:[6,8,10],8:[9],10:[8,11,12,13],14:[12,18],15:[13,14,16],17:[18,16]}) show(g) print ptsk(g,0) leafstripsimult(g) 
       
6
6
6
6
# Example PT.4-3: Independent twins v and w in G and Zskew(G-v)=0. # This is the graph constructed from Example PT.4-2 by twinning vertex 18. 
       
g=Graph({2:[1,4,5,6],4:[3,5,6],5:[6,18,19],7:[6,8,10],8:[9],10:[8,11,12,13],14:[12,18,19],15:[13,14,16],17:[18,16,19]}) show(g) z=Zskew(g) print z, ptsk(g,z) 
       
1 6
1 6
# Example PT.4-4: the skew propagation time interval need not be not full. 
       
g=Graph({1:[2,5,6],2:[7,3],4:[3,5]}) show(g) for v in [1..7]: print v, prop_time_unlooped(g,[v]) 
       
1 -1
2 -1
3 2
4 -1
5 2
6 4
7 4
1 -1
2 -1
3 2
4 -1
5 2
6 4
7 4
#################################### # Section 10.5 # Power propagation time # Power propagation time of a set S, minimum power propagation time of a set of k elements, # and power propagation time of a graph are illustrated. # These are computed by functions pptS(G,S), pptk(G,k), and ppt(G). ################################### 
       
# Example PT.5-1: ppt(T_1)=3. 
       
T1=Graph({2:[1,3],3:[2,4,5],5:[6,7],7:[8,9],10:[9,11,12]}) show(T1) 
       
PowerDom(T1) 
       
2
2
pptS(T1,[3,10]) 
       
3
3
ppt(T1) 
       
3
3
#################################### # Section 11.2 # Standard throttling number is computed by the functions # th(G) and thk(G); the latter retuns both th(G) and least k that achieves it. # There are examples for th(G-v), th(G-e), th(G/e). ################################### 
       
# Example TH.2-1: th(T_1). 
       
T1=Graph({2:[1,3],3:[2,4,5],5:[6,7],7:[8,9],10:[9,11,12]}) show(T1) th(T1) 
       
7
7
# Example TH.2-2: Throttling number of the wheel W_7. 
       
g=graphs.WheelGraph(7) show(g) zfs=zero_forcing_set_bruteforce(g) z=len(zfs) print 'zfs=', zfs print 'Z(G)=',z print 'pt(G)=', pt(g) t=ptk(g,z) thg=z+t for k in [z..z+t]: ptkg=ptk(g,k) print k, ptkg, 'th(G,k)=',k+ptkg if k+ptkg<thg: thg=k+ptkg print 'th(G)=',thg 
       
zfs= {0, 1, 2}
Z(G)= 3
pt(G)= 2
3 2 th(G,k)= 5
4 2 th(G,k)= 6
5 1 th(G,k)= 6
th(G)= 5
zfs= {0, 1, 2}
Z(G)= 3
pt(G)= 2
3 2 th(G,k)= 5
4 2 th(G,k)= 6
5 1 th(G,k)= 6
th(G)= 5
# All optimal throttling sets have 3 vertices. Only the minimum zero forcing sets with the dominating vertex have propagation time 2. 
       
# Example TH.2-3: th(G-v)= th(G)+1. 
       
k3=graphs.CompleteGraph(3) p3=graphs.PathGraph(3) g=p3.cartesian_product(k3) g.relabel() show(g) print 'pt(G)=',pt(g),'th(G)=',th(g) 
       
pt(G)= 2 th(G)= 5
pt(G)= 2 th(G)= 5
g.delete_vertex(5) show(g) th(g) 
       
6
6
Z(g) 
       
2
2
# Example TH.2-4: th(G-v)=th(G) where v is an adjacent twin. 
       
g=graphs.FriendshipGraph(2) show(g) th(g) 
       
4
4
g.delete_vertex(0) show(g) th(g) 
       
4
4
# Example TH.2-5: th(G-e)=th(G)-1. 
       
# this is G g=graphs.CycleGraph(4) g.add_edge(0,2) show(g) th(g) 
       
4
4
# this is G-e g.delete_edge(0,2) show(g) th(g) 
       
3
3
# Example TH.2-6: th(G-e)= th(G)+1. 
       
g=graphs.FriendshipGraph(3) show(g) th(g) 
       
5
5
g.delete_edge(1,2) show(g) th(g) 
       
6
6
# Example TH.2-7: th(G/e)=th(G)+1. 
       
k3=graphs.CompleteGraph(3) p3=graphs.PathGraph(3) g=p3.cartesian_product(k3) g.relabel() show(g) th(g) 
       
5
5
g.contract_edge(4,5) show(g) th(g) 
       
6
6
#################################### # Section 11.3 # PSD throttling number is computed by the functions # thp(G)=th_+(G) and thpk(G); the latter retuns both th_+(G) and least k that achieves it. # There are examples for tight bounds on vertex and edge deletion. ################################### 
       
# Example TH.3-1: T_1 
       
T1=Graph({2:[1,3],3:[2,4,5],5:[6,7],7:[8,9],10:[9,11,12]}) show(T1) thp(T1) 
       
4
4
# Example TH.3-2: th_+(G-v) = th_+(G)-1 
       
nn=11 g=graphs.PathGraph(nn) gv=graphs.PathGraph(nn-1) print 'th_+(G) =',thp(g), 'th_+(G-v) =',thp(gv) 
       
th_+(G) = 5 th_+(G-v) = 4
th_+(G) = 5 th_+(G-v) = 4
# Example TH.3-3: th_+(G-e) = th_+(G)-1 
       
p4=graphs.PathGraph(4) k1=graphs.CompleteGraph(1) g=p4.join(k1) g.relabel() show(g) thp(g) 
       
4
4
g.delete_edge(1,4) show(g) thp(g) 
       
3
3
# Example TH.3-4: th_+(G-e) = th_+(G)+1 
       
g=graphs.CompleteBipartiteGraph(1,5) show(g) thp(g) 
       
2
2
g.delete_edge(0,1) show(g) thp(g) 
       
3
3
#################################### # Section 11.4 # Skew throttling number is computed by the functions # ths(G) and thsk(G); the latter retuns both th_-(G) and least k that achieves it. # There are examples for tight bounds on vertex and edge deletion. ################################### 
       
# Example TH.4-1: T_1 
       
T1=Graph({2:[1,3],3:[2,4,5],5:[6,7],7:[8,9],10:[9,11,12]}) show(T1) ths(T1) 
       
4
4
# Example TH.4-2: th_-(G-e)=th_-(G)-2 
       
p4=graphs.PathGraph(4) c4=graphs.CycleGraph(4) c4.relabel([4,5,6,7]) g=p4.union(c4) show(g) ths(g) 
       
4
4
g.delete_edge(4,7) show(g) ths(g) 
       
2
2
# Example TH.4-3: th_-(G/e) = th_-(G)+1 and th_-(G_e)=th_-(G)-1 
       
g=Graph({0:[1,2,3],1:[4]}) show(g) print Zskew(g), ths(g) 
       
1 3
1 3
g.subdivide_edge(0,2,1) show(g) print Zskew(g), ths(g) 
       
0 2
0 2
Observe that contracting edge {0,5} in the subdivided graph give the original graph 
       
#################################### # Section 11.5 # Power domination throttling is computed by the functions # thpd(G) and thpdk(G); the latter retuns both th_pd(G) and least k that achieves it. ################################### 
       
# Example TH.5-1: th_pd(T_1)=3 for T_1 in Example 11.5-1 
       
show(T1) thpd(T1) 
       
5
5
# Example TH.5-2: Connected graph with relatively high power domination throttling number 
       
g=Graph({1:[2,3,4,5],2:[3,6,7],3:[8,9],9:[10]}) show(g) thpd(g) 
       
5
5
# Example TH.5-3: th_pd(G-e)=th_pd(G)-1 
       
g=Graph({1:[2,4,6,9],10:[9,11],7:[6,8],3:[2,4],5:[2,4]}) show(g) thpd(g) 
       
5
5
g.delete_edge(2,3) show(g) thpd(g) 
       
4
4
# Example TH.5-4: th_pd(G-e)=th_pd(G)+1 
       
p4=graphs.PathGraph(4) k1=graphs.CompleteGraph(1) g=p4.join(k1) g.relabel() show(g) thpd(g) 
       
2
2
g.delete_edge(1,4) show(g) thpd(g) 
       
3
3
# Example TH.5-5: th_pd(G/e)=th_pd(G)+1 
       
g=Graph({1:[2,3,4,6],5:[4,10],6:[7,8,9],9:[4,10],11:[2],12:[7]}) show(g) thpd(g) 
       
4
4
g.contract_edge(4,9) show(g) thpd(g) 
       
5
5
#################################### # FUNCTIONS-HERE # These are the main function cells: evaluate each one #################################### 
       
# F-1 load main zero forcing functions # zero forcing, PSD forcing, and minimum rank code # developed by S. Butler, L. DeLoss, J. Grout, H.T. Hall, J. LaGrange, J.C.-H. Lin,T. McKay, J. Smith, G. Tims # updated and maintained by Jephian C.-H. Lin URL='https://raw.githubusercontent.com/jephianlin/mr_JG/py2/' files=['Zq_c.pyx','Zq.py','zero_forcing_64.pyx','zero_forcing_wavefront.pyx','minrank.py', 'inertia.py'] for f in files: print("Loading %s..."%f); load(URL+f) 
       
Loading Zq_c.pyx...
Compiling
/home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_ndIHkr.pyx..\
.
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling
/home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_sAhwEF.pyx..\
.
Loading zero_forcing_wavefront.pyx...
Compiling
/home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_tDPwHY.pyx..\
.
Loading minrank.py...
Loading inertia.py...
Loading Zq_c.pyx...
Compiling /home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_ndIHkr.pyx...
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling /home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_sAhwEF.pyx...
Loading zero_forcing_wavefront.pyx...
Compiling /home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_tDPwHY.pyx...
Loading minrank.py...
Loading inertia.py...
# F-2 other zero forcing functions # # Zero forcing number Z(G) def Z(g): z=len(zero_forcing_set_bruteforce(g)) return z # # All zero forcing sets of G of size k # input: a graph G and a positive integer k # output: a list of all zero forcing sets G of size k def ZFsets(G,k): ord=G.order() V = G.vertices() S = subsets(V,k) ZFS=[] for s in S: A=zerosgame(G,s) if len(A)==ord: ZFS.append(s) return ZFS 
       
# F-3 Function for propagation time of set S, pt(G,S) # adapted by Leslie Hogben from Steve Butler's skew propagation time code (F-10) # input: a graph G and a set S of vertices # output: pt(G,S) = prop time of S in G (if -1 is returned then S is not a zero forcing set) def ptz(G,S): V=set(G.vertices()) count = 0 done = False active = set(S) filled = set(S) for v in V: N=set(G.neighbors(v)) if v in active and N.issubset(filled): active.remove(v) if (v in filled) and (v not in active) and (len(N.intersection(filled)) == G.degree(v)-1): active.add(v) while not done: done = True new_active = copy(active) new_filled = copy(filled) for v in active: N=set(G.neighbors(v)) if len(N.intersection(filled)) == G.degree(v)-1: if done: done = False count += 1 N.symmetric_difference_update(N.intersection(filled)) u=N.pop() new_active.remove(v) new_active.add(u) new_filled.add(u) active = copy(new_active) filled = copy(new_filled) # print filled for v in V: N=set(G.neighbors(v)) if v in active and N.issubset(filled): active.remove(v) if (v in filled) and (v not in active) and (len(N.intersection(filled)) == G.degree(v)-1): active.add(v) if len(filled)==len(V): return count return -1 
       
# F-4 Functions for standard propagation time # ptk(G,k) = min prop time over zero forcing sets of G of size k # adapted by Leslie Hogben from Steve Butler's skew propagation time code (F-10) # this is the function used to compute propation time of G and also throttling # input: a graph G and a positive integer k >= Z(G) # output: pt(G,k) = min prop time over zero forcing sets G of size k (if order+1 is returned then k<Z(G)) def ptk(G,k): V=set(G.vertices()) sets = subsets(V,k) pt=len(V)+1 for S in sets: count = 0 done = False active = set(S) filled = set(S) for v in V: N=set(G.neighbors(v)) if v in active and N.issubset(filled): active.remove(v) if (v in filled) and (v not in active) and (len(N.intersection(filled)) == G.degree(v)-1): active.add(v) while not done: done = True new_active = copy(active) new_filled = copy(filled) for v in active: N=set(G.neighbors(v)) if len(N.intersection(filled)) == G.degree(v)-1: if done: done = False count += 1 N.symmetric_difference_update(N.intersection(filled)) u=N.pop() new_active.remove(v) new_active.add(u) new_filled.add(u) active = copy(new_active) filled = copy(new_filled) # print filled for v in V: N=set(G.neighbors(v)) if v in active and N.issubset(filled): active.remove(v) if (v in filled) and (v not in active) and (len(N.intersection(filled)) == G.degree(v)-1): active.add(v) if len(filled)==len(V): if count < pt: pt=count return pt # propagation time pt(G) # input: a graph G # output: propagation time pt(G)=pt(G,Z(G)) def pt(g): ptg=ptk(g,Z(g)) return ptg 
       
# F-5 Functions for standard throttling number # # These functions compute th(G) by minimizing th(G,k) # # input: a graph G # output: standard throttling number th(G) def th(g): z=Z(g) t=ptk(g,z) thz=z+t kmax=min(z+t-1,g.order()) for k in [z..kmax]: ptkg=ptk(g,k) if k+ptkg<thz: thz=k+ptkg return thz # input: a graph G # output: the set [th(G),ko] such that th(G)= th(G,ko) def thk(g): z=Z(g) ko=z t=ptk(g,z) thz=z+t kmax=min(z+t-1,g.order()) for k in [z..kmax]: ptkg=ptk(g,k) if k+ptkg<thz: thz=k+ptkg ko=k return [thz,ko] 
       
# F-6 Load PSD propagation time funtions, including pt_plus(G,S) = PSD propagation time of a set S # code by Nathan Warnberg # updated and maintained by Jephian C.-H. Lin load('https://raw.githubusercontent.com/jephianlin/zero_forcing/master/psd_prop_time_interval.py') 
       
xrange test passed
Loading Zq_c.pyx...
Compiling
/home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_p_lJwg.pyx..\
.
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling
/home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_lap2dl.pyx..\
.
Loading zero_forcing_wavefront.pyx...
Compiling
/home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_0hgkjK.pyx..\
.
Loading minrank.py...
Loading inertia.py...
xrange test passed
Loading Zq_c.pyx...
Compiling /home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_p_lJwg.pyx...
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling /home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_lap2dl.pyx...
Loading zero_forcing_wavefront.pyx...
Compiling /home/sageuser/8/.sage/temp/sage.math.iastate.edu/16895/tmp_0hgkjK.pyx...
Loading minrank.py...
Loading inertia.py...
# F-7 Functions for PSD propagation time # function to compute pt_+(G,k) = minimum PSD propagation time over sets of size k # input: a graph G and a positive integer k # output: pt_+(G,k) def ptpk(G,k): ord=G.order() V = G.vertices() S = subsets(V,k) ptp = -1 for s in S: ptps=pt_plus(G,s) if (ptp < 0): ptp=ptps if (ptps >= 0) and (ptps < ptp): ptp=ptps return ptp # function to compute pt_+(G) = pt_+(G,Z_+(G)) # input: a graph G # output: pt_+(G) def ptp(G): ptpg=ptpk(G,Zplus(G)) return ptpg 
       
# F-8 Functions for PSD throttling number # This function computes th_+(G)=thp(G) by minimizing thp(G,k) # input: a graph g # output: PSD throttling number thp(g) def thp(g): z=Zplus(g) t=ptpk(g,z) thz=z+t kmax=min(z+t-1,g.order()) for k in [z..kmax]: ptkg=ptpk(g,k) if k+ptkg<thz: thz=k+ptkg return thz # This function computes th_+(G)=thp(g) by minimizing thp(g,k) and # returns both th_+(G) and the least k such that th_+(G) = th_+(G,k) # input: a graph g # output: the set [thp(g),ko] such that thp(g)= thp(g,ko) def thpk(g): z=Zplus(g) ko=z t=ptpk(g,z) thz=z+t kmax=min(z+t-1,g.order()) for k in [z..kmax]: ptkg=ptpk(g,k) if k+ptkg<thz: thz=k+ptkg ko=k return [thz,ko] 
       
# F-9 main function cell zero forcing number loop graphs # skew forcing uses no loops # written by Jephian Lin def gzerosgame(g,F=[],B=[]): """ Return the derived set for a given graph g with set of banned edges B and a initial set of vertices. The derived set is given by doing generalized zero forcing process. That is, if y is the only white neighbor of x and xy is not banned, then x could force y into black. Input: g: a simple graph F: a list of vertices of g B: a list of tuples representing banned edges of g Output: A set of black vertices when zero forcing process stops. Examples: sage: gzerosgame(graphs.PathGraph(5),[0]) set([0, 1, 2, 3, 4]) sage: gzerosgame(graphs.PathGraph(5),[0],[(1,2)]) set([0, 1]) """ S=set(F) # suspicuous vertices Black_vertices=set(F) # current black vertices again=1 # iterate again or not while again==1: again=0 for x in S: N=set(g.neighbors(x)) D=N.difference(Black_vertices) # set of white neighbors if len(D)==1: for v in D: y=v # the only white neighbor if (((x,y) in B)==False) and (((y,x) in B)==False): again=1 S.remove(x) S.add(y) Black_vertices.add(y) break return(Black_vertices) def gZ_leq(graph, support=[], bannedset=[],i=None): """ For a given graph with support and banned set, if there is a zero forcing set of size i then return it; otherwise return False. Input: graph: a simple graph support: a list of vertices of g bannedset: a list of tuples representing banned edges of graph i: an integer, the function check gZ <= i or not Output: if F is a zero forcing set of size i and support is a subset of F, then return F False otherwise Examples: sage: gZ_leq(graphs.PathGraph(5),[],[],1) set([0]) sage: gZ_leq(graphs.PathGraph(5),[],[(0,1)],1) False """ if i < len(support): # print 'i cannot less than the cardinality of support' return False j=i-len(support) # additional number of black vertices VX=graph.vertices() order=graph.order() for y in support: VX.remove(y) # VX is the vertices outside support now for subset in Subsets(VX,j): test_set=set(support).union(subset) # the set is tested to be a zero forcing set outcome=gzerosgame(graph, test_set, bannedset) if len(outcome)==order: return test_set return False def find_gzfs(graph, support=[], bannedset=[], upper_bound=None, lower_bound=None): """ For a given graph with support and banned set, return the an optimal generalized zero forcing set. If upper_bound is less than the generalized zero forcing number then return ['wrong']. If lower_bound is greater than the generalized zero forcing number then the return value will not be correct Input: graph: a simple graph support: a list of vertices of g bannedset: a list of tuples representing banned edges of graph upper_bound: an integer supposed to be an upper bound of gZ. lower_bound: an integer supposed to be a lower bound of gZ. The two bounds may shorten the computation time. But one may leave it as default value if one is not sure. Output: if F is an optimal zero forcing set of size i then return F. If upper_bound is less than the general zero forcing number then return ['wrong']. Examples: sage: find_gzfs(graphs.PathGraph(5)) set([0]) sage: find_gzfs(graphs.PathGraph(5),[1],[(3,2)]) set([0, 1, 3]) """ VX=graph.vertices() order=graph.order() s=len(support) for y in support: VX.remove(y) # VX is the vertices outside support now if upper_bound==None: upper_bound=order # the default upper bound if lower_bound==None: lower_bound=len(VX) # temporary lower bound for v in VX: N=set(graph.neighbors(v)) D=N.difference(support) lower_bound=min([lower_bound,len(D)]) for v in support: N=set(graph.neighbors(v)) D=N.difference(support) lower_bound=min([lower_bound,len(D)-1]) lower_bound=lower_bound+s # the default lower bound i=upper_bound find=1 # does sage find a zero forcing set of size i outcome=['wrong'] # default outcome while i>=lower_bound and find==1: find=0 leq=gZ_leq(graph, support, bannedset,i) # check gZ <= i or not if leq!=False: outcome=leq find=1 i=i-1 return outcome def find_gZ(graph, support=[], bannedset=[], upper_bound=None, lower_bound=None): """ For a given graph with support and banned set, return the zero. upper_bound and lower_bound could be left as default value if one is not sure. Input: graph: a simple graph support: a list of vertices of g bannedset: a list of tuples representing banned edges of graph upper_bound: an integer supposed to be an upper bound of gZ. lower_bound: an integer supposed to be a lower bound of gZ. The two bounds may shorten the computation time. But one may leave it as default value if one is not sure. Output: the generalized zero forcing number Examples: sage: find_gZ(graphs.PathGraph(5)) 1 sage: find_gZ(graphs.PathGraph(5),[1],[(3,2)]) 3 """ return len(find_gzfs(graph, support, bannedset, upper_bound, lower_bound)) def X(g): """ For a given graph g, return the verices set X of a part of the bipartite used to compute the exhaustive zero forcing number. Input: g: a simple graph Output: a list of tuples ('a',i) for all vertices i of g Examples: sage: X(graphs.PathGraph(5)) [('a', 0), ('a', 1), ('a', 2), ('a', 3), ('a', 4)] """ return [('a',i) for i in g.vertices()] def Y(g): """ For a given graph g, return the verices set Y of the other part of the bipartite used to compute the exhaustive zero forcing number. Input: g: a simple graph Output: a list of tuples ('b',i) for all vertices i of g Examples: sage: Y(graphs.PathGraph(5)) [('b', 0), ('b', 1), ('b', 2), ('b', 3), ('b', 4)] """ return [('b',i) for i in g.vertices()] def tilde_bipartite(g,I=[]): """ For a given graph g and an index set I, return the bipartite graph \widetilde{G}_I used to compute the exhaustive zero forcing number. Input: g: a simple graph I: a list of vertices of g Output: the bipartite graph \widetilde{G}_I Examples: sage: h=tilde_bipartite(graphs.PathGraph(5),[1]) sage: h.vertices() [('a', 0), ('a', 1), ('a', 2), ('a', 3), ('a', 4), ('b', 0), ('b', 1), ('b', 2), ('b', 3), ('b', 4)] sage: h.edges() [(('a', 0), ('b', 1), None), (('a', 1), ('b', 0), None), (('a', 1), ('b', 1), None), (('a', 1), ('b', 2), None), (('a', 2), ('b', 1), None), (('a', 2), ('b', 3), None), (('a', 3), ('b', 2), None), (('a', 3), ('b', 4), None), (('a', 4), ('b', 3), None)] """ E0=[(('a',i), ('b',i)) for i in I] # edges given by I E1=[] # invariant edges for i in g.vertices(): for j in g.neighbors(i): E1.append((('a',i),('b',j))) h=Graph() h.add_vertices(X(g)) h.add_vertices(Y(g)) h.add_edges(E0) h.add_edges(E1) # h=(X union Y, E0 union E1) return h def find_EZ(g,bound=None): """ For a given graph g, return the exhaustive zero forcing number of g. A given bound may shorten the computation. Input: g: a simple graph bound: a integer as an upper bound. It could be left as default value if one is not sure. Output: the exhaustive zero forcing number (EZ) of g Examples: sage: find_EZ(graphs.PathGraph(5)) 1 sage: h=graphs.CycleGraph(5) sage: h.add_vertices([5,6,7,8,9]) sage: h.add_edges([(0,5),(1,6),(2,7),(3,8),(4,9)]) sage: find_EZ(h) # the EZ of a 5-sun 2 """ order=g.order() Z=find_gZ(g) # without support and banned set, the value is the original zero forcing number if bound==None: bound=Z # default upper bound gZ_bound=bound+order V=set(g.vertices()) e=-1 # temporary output for I in Subsets(V): leq=gZ_leq(tilde_bipartite(g,I),Y(g),[],e) # this avoid abundant computation if leq==False: e=find_gZ(tilde_bipartite(g,I),Y(g),[],gZ_bound,e+1) # in this case, we already know e+1-order<=gZ-order<=bound and so e+1<=gZ<=gZ_bound if e==gZ_bound: break return e-order # EZ=max-order def find_loopedZ(g,I): """ For a given graph g and the index of the vertices with loops, return the zero forcing number of this looped graph. Input: g: a simple graph, the underlying graph of the looped graph. I: the index of the vertices with loops. Output: the zero forcing number of this looped graph. Examples: sage: g = Graph({0:[1],1:[0]}); sage: I=[0,1]; sage: find_loopedZ(g,I) 1 sage: g = Graph({0:[1],1:[0]}); sage: I=[0]; sage: find_loopedZ(g,I) 0 """ return find_gZ(tilde_bipartite(g,I),Y(g),[])-g.order() def Zskew(g): # report Z_-(G) return find_loopedZ(g,[]) 
       
# F-10 Function for skew propagation time of a set # written by Steve Butler def prop_time_unlooped(G,S): """ Determine the propagation time of S in G assuming all the vertices are unlooped. """ V=set(G.vertices()) count = 0 done = False active = set(S) filled = set(S) for v in V: N=set(G.neighbors(v)) if v in active and N.issubset(filled): active.remove(v) if (v not in active) and (len(N.intersection(filled)) == G.degree(v)-1): active.add(v) while not done: done = True new_active = copy(active) new_filled = copy(filled) for v in active: N=set(G.neighbors(v)) if len(N.intersection(filled)) == G.degree(v)-1: if done: done = False count += 1 N.symmetric_difference_update(N.intersection(filled)) u=N.pop() new_active.remove(v) new_active.add(u) new_filled.add(u) active = copy(new_active) filled = copy(new_filled) for v in V: N=set(G.neighbors(v)) if v in active and N.issubset(filled): active.remove(v) if (v not in active) and (len(N.intersection(filled)) == G.degree(v)-1): active.add(v) if len(filled)==len(V): return count return -1 
       
# F-11 Functions for skew propagation time # function to compute pt_-(G,k) = minimum skew propagation time over sets of size k # input: a graph G and a positive integer k # output: pt_-(G,k) def ptsk(G,k): ord=G.order() V = G.vertices() S = subsets(V,k) ptm = -1 for s in S: ptms=prop_time_unlooped(G,s) if (ptm < 0): ptm=ptms if (ptms >= 0) and (ptms < ptm): ptm=ptms return ptm # function to compute pt_-(G) = pt_-(G,Zskew(G)) # input: a graph G and a positive integer k # output: pt_-(G,k) def pts(G): return ptsk(G,Zskew(G)) 
       
# F-12 Function for simultaneous leaf stripping algorithm # # input: a graph g # output: set of leaves and set of neigbors of leaves (if two leaves have same neighbor, neighbors set is smaller) def leafs(g): ord=g.order() k=0 leafs=[] nhbrs=[] for k in g.degree_iterator(labels=True): if k[1]==1: leafs.append(k[0]) nhbrs.append(g.neighbors(k[0])[0]) return leafs, nhbrs # # simultaneous leaf stripping to find skew propagaion time when Z_-(G)=0 # input: a graph h # output: pt_-(h) if Z_-(h)=0; -1 if Z_-(h)>0. def leafstripsimult(h): g=h.copy() hasleaf=True k=0 # k is the iteration counter while hasleaf: lfnbr=leafs(g) # lf is set of leaves and nehgbors k=k+1 # another iteration is running ord=g.order() ds=g.degree_sequence() # degree sequence is nonincreasing if ds[ord-1]==0: # isolated vertex return -1 if len(lfnbr[0])==0: # no leaf return -1 x=set(lfnbr[1]) # this is needed to remove repetitions in the list if len(lfnbr[0])!=len(x): # too many leaves for 1 neighbor return -1 if ds[0]==1: # successful completion but rK_2 contributes 1 not 2 to prop time return 2*k-1 g.delete_vertices(lfnbr[0]) # delete leaves g.delete_vertices(lfnbr[1]) # delete neighbors if g.order()==0: # successful completion return 2*k 
       
# F-13 Functions for skew throttling number # This function computes ths(g)=th_-(g) by minimizing ths(g,k) # input: a graph g # output: skew throttling number th(g) def ths(g): z=Zskew(g) t=ptsk(g,z) thz=z+t kmax=min(z+t-1,g.order()) for k in [z+1..kmax]: ptkg=ptsk(g,k) if k+ptkg<thz: thz=k+ptkg return thz # This function computes th_-(G)=ths(G) by minimizing ths(G,k) and # returns both th_-(G) and the least k such that th_-(G) = th_-(G,k) # input: a graph g # output: the set [ths(g),ko] such that ths(g)= ths(g,ko) def thsk(g): z=Zskew(g) ko=z t=ptsk(g,z) thz=z+t kmax=min(z+t-1,g.order()) for k in [z+1..kmax]: ptkg=ptsk(g,k) if k+ptkg<thz: thz=k+ptkg ko=k return [thz,ko] 
       
# WARNING: Power domination fuctions depend on zero forcing functions # Enter F-1, F-2, F-3, F-4, F-5 before entering the cells below 
       
# F-14 Power dominating set functions # by Brian Wissman # input: a graph G and a subset of its vertices V # output: true/false depending if V is a power dominating set of G def isPowerDominatingSet(G,V): N=[] for i in V: N+=G.neighbors(i) NVert=uniq(N+list(V)) A=zerosgame(G,NVert) if len(A)==G.order(): return True else: return False # input: a graph G # output: a power dominating set of G that achieves the power domination number def minPowerDominatingSet(G): V = G.vertices() for i in range(1,len(V)+1): S = subsets(V,i) for s in S: if isPowerDominatingSet(G,s): return s # input: a graph G # output: the power domination number of G, gamma_P(G) def PowerDom(G): V = G.vertices() for i in range(1,len(V)+1): S = subsets(V,i) for s in S: if isPowerDominatingSet(G,s): p=len(s) return p 
       
# F-15 All power dominating sets of size k # # input: a graph G and a positive integer k # output: a list of all power dominating sets G of size k def PDsets(G,k): ord=G.order() V = G.vertices() S = subsets(V,k) PDS=[] for s in S: if isPowerDominatingSet(G,s): PDS.append(s) return PDS 
       
# F-16 Functions for power propagation time # function to compute closed neighborhood N[S] # input: a graph G and a list S of vertices # output: set N[S] def nhbd(G,S): NS=set(S) for x in S: nx=G.neighbors(x) NS=NS.union(set(nx)) return NS # function to compute ppt(G,S) = power propagation time of a set S # input: a graph G and a subset S of vertices # output: ppt(G,S) def pptS(G,S): ord=G.order() V = G.vertices() NS=nhbd(G,S) ptNS=ptz(G,NS) if ptNS>-1: ptNS=ptNS+1 return ptNS # function to compute ppt(G,k) = minimum power propagation time over sets of size k # input: a graph G and a positive integer k # output: ppt(G,k) def pptk(G,k): ord=G.order() V = G.vertices() S = subsets(V,k) ptm = -1 for s in S: ptms=pptS(G,s) if (ptm < 0): ptm=ptms if (ptms >= 0) and (ptms < ptm): ptm=ptms return ptm # function to compute ppt(G) = ppt(G,gamma_P(G)) # input: a graph G and a positive integer k # output: ppt(G) def ppt(G): return pptk(G,PowerDom(G)) 
       
# F-17 Functions for power domination throttling number # # Compute domination number of graph g def dom(g): gds=g.dominating_set() ds=len(gds) return ds # # This function computes th_PD(G) by minimizing th_PD(G,k) for # gamma_P(G)<=k <= gamma(G) # # input: a graph G # output: power domination throttling number th_PD(G) def thpd(g): z=PowerDom(g) d=dom(g) t=pptk(g,z) thz=z+t ko=z for k in [z+1..d]: ptkg=pptk(g,k) if k+ptkg<thz: thz=k+ptkg ko=k return thz # This function computes th_PD(G) by minimizing th_PD(G,k) for # gamma_P(G)<=k <= gamma(G) # and reports both # th_PD^x(G) and ko such that th_PD(G)=th_PD^x(G,ko) # # input: a graph G and a bound b where it is known th_PD^x(G) <= b # (if know be known, use n) # output: returns power product throttling number th_PD^x(G) # and prints th_PD^x(G) and least k such that th_PD^x(G) = th_PD^x(G,k) def thpdk(g): z=PowerDom(g) d=dom(g) t=pptk(g,z) thz=z+t ko=z for k in [z+1..d]: ptkg=pptk(g,k) if k+ptkg<thz: thz=k+ptkg ko=k print 'th_PD(G)', thz, 'k',ko return thz