ACFHMTW-Product_Power_Throttling

1224 days ago by hogben

#################################### # This Sage worksheet is written by Leslie Hogben # based on code written by numerous other people, some of whom are listed. # It contains examples from the paper # Product throttling for power domination # by Sarah Anderson, Karen Collins, Daniela Ferrero, Leslie Hogben, Carolyn Mayer, Ann Trenk, Shanise Walker # Please report errors to hogben@aimath.org #################################### 
       
#################################### # To run the examples, you must got to the bottom of this sheet # and run all the function cells before running the examples. #################################### 
       
#################################### # Example 2.4 #################################### 
       
g=Graph({1:[0,2],3:[0,4],5:[0,6],7:[0,8],9:[0,10],0:[11],12:[11,13],14:[13,15],15:[14,16],17:[16]}) show(g) print 'dom set=', g.dominating_set(), 'gamma_p=',PowerDom(g), 'ppt(g)=',ppt(g) print thpdxb(g,18) 
       
dom set= [0, 1, 3, 5, 7, 9, 13, 16] gamma_p= 1 ppt(g)= 7
th_PD^x(G) 4 k 2
4
dom set= [0, 1, 3, 5, 7, 9, 13, 16] gamma_p= 1 ppt(g)= 7
th_PD^x(G) 4 k 2
4
#################################### # Example 6.2 #################################### 
       
p=Graph({18:[19]}) gp=p.cartesian_product(g) gp.relabel() show(gp) 
       
gds=gp.dominating_set() print 'dom set=', gds,'gamma(gp)=',len(gds) print 'gamma_p=',PowerDom(gp), 'ppt(gp)=',ppt(gp) 
       
dom set= [0, 11, 15, 20, 22, 24, 26, 28, 31, 35] gamma(gp)= 10
gamma_p= 2 ppt(gp)= 7
dom set= [0, 11, 15, 20, 22, 24, 26, 28, 31, 35] gamma(gp)= 10
gamma_p= 2 ppt(gp)= 7
pptS(gp,[0,18,15,33]) 
       
2
2
# Thus th_PD^x(gp) <= 8. # By Remark 3.8(2), need consider 2<= k <= 3 for th_PD^x(gp,k) = th_PD^x(gp) # Since ppt(gp,4) >= 2 and ppt(gp,{0,18,15,33})=2, ppt(gp,4)=2. # As above, ppt(gp,2)=7. 
       
'ppt(gp,3)=',pptk(gp,3) 
       
('ppt(gp,3)=', 4)
('ppt(gp,3)=', 4)
# Thus th_PD^x(gp) = 8. 
       
#################################### # Example 6.3 #################################### 
       
g=Graph({9:[1,17],2:[1,3,10],3:[4,11],4:[5,12],5:[6,13],6:[7,14],7:[8,15],8:[1,16]}) show(g) print 'dom set=', g.dominating_set(), 'gamma_p=',PowerDom(g), 'ppt(g)=',ppt(g) print thpdxb(g,17) 
       
dom set= [2, 3, 4, 5, 6, 7, 8, 9] gamma_p= 3 ppt(g)= 2
th_PD^x(G) 6 k 3
6
dom set= [2, 3, 4, 5, 6, 7, 8, 9] gamma_p= 3 ppt(g)= 2
th_PD^x(G) 6 k 3
6
# Thus th_PD^x(g) <= 6. 
       
p=Graph({18:[19]}) gp=p.cartesian_product(g) gp.relabel([1..34]) show(gp) 
       
gds=gp.dominating_set() print 'dom set=', gds,'gamma(gp)=',len(gds) print 'gamma_p=',PowerDom(gp), 'ppt(gp)=',ppt(gp) 
       
dom set= [1, 4, 11, 13, 15, 21, 23, 27, 31, 33, 34] gamma(gp)= 11
gamma_p= 3 ppt(gp)= 7
dom set= [1, 4, 11, 13, 15, 21, 23, 27, 31, 33, 34] gamma(gp)= 11
gamma_p= 3 ppt(gp)= 7
pptS(gp,[1,5,18,20,24]) 
       
2
2
# Thus th_PD^x(gp) <= 10. # By Remark 3.8(2), need consider 3<= k <= 5 for th_PD^x(gp,k) = th_PD^x(gp) # and k=3, k=5 have been done. 
       
'ppt(gp,4)=',pptk(gp,4) 
       
('ppt(gp,4)=', 3)
('ppt(gp,4)=', 3)
# Thus th_PD^x(gp) = 10. 
       
 
       
#################################### # Theorem 6.11 #################################### 
       
#################################### # Cartesian products of Pn[]Pm for n = 4,5,6 # # For n = 4 case, check m = 5,6,9 # For n = 5 case, check m = 8,12 # For n = 6 case, check m = 6,10 #################################### 
       
p4 = graphs.PathGraph(4) p5 = graphs.PathGraph(5) p6 = graphs.PathGraph(6) c4 = graphs.CycleGraph(4) c5 = graphs.CycleGraph(5) c6 = graphs.CycleGraph(6) 
       
# n = 4: 
       
print 'domination number p4[]p5 =', dom(p4.cartesian_product(p5)) print 'throttle product p4[]p5=', thpdx(p4.cartesian_product(p5)) print 'domination number p4[]p6=', dom(p4.cartesian_product(p6)) print 'throttle product p4[]p6=', thpdx(p4.cartesian_product(p6)) print 'domination number p4[]p9=', dom(p4.cartesian_product(graphs.PathGraph(9))) print 'throttle product p4[]p9=', thpdx(p4.cartesian_product(graphs.PathGraph(9))) 
       
domination number p4[]p5 = 6
throttle product p4[]p5= 6
domination number p4[]p6= 7
throttle product p4[]p6= 7
domination number p4[]p9= 10
throttle product p4[]p9= 10
domination number p4[]p5 = 6
throttle product p4[]p5= 6
domination number p4[]p6= 7
throttle product p4[]p6= 7
domination number p4[]p9= 10
throttle product p4[]p9= 10
print 'domination number c4[]c5 =', dom(c4.cartesian_product(graphs.CycleGraph(5))) print 'domination number c4[]c6 =', dom(c4.cartesian_product(graphs.CycleGraph(6))) print 'domination number c4[]c9 =', dom(c4.cartesian_product(graphs.CycleGraph(9))) print 'domination number p4[]c6 =', dom(p4.cartesian_product(graphs.CycleGraph(6))) 
       
domination number c4[]c5 = 5
domination number c4[]c6 = 6
domination number c4[]c9 = 9
domination number p4[]c6 = 6
domination number c4[]c5 = 5
domination number c4[]c6 = 6
domination number c4[]c9 = 9
domination number p4[]c6 = 6
print 'domination number p4[]c5 =', dom(p4.cartesian_product(graphs.CycleGraph(5))) print 'throttle product p4[]c5=', thpdx(p4.cartesian_product(graphs.CycleGraph(5))) print 'domination number p4[]c9=', dom(p4.cartesian_product(graphs.CycleGraph(9))) print 'throttle product p4[]c9=', thpdx(p4.cartesian_product(graphs.CycleGraph(9))) 
       
domination number p4[]c5 = 6
throttle product p4[]c5= 6
domination number p4[]c9= 10
throttle product p4[]c9= 10
domination number p4[]c5 = 6
throttle product p4[]c5= 6
domination number p4[]c9= 10
throttle product p4[]c9= 10
# n = 5: 
       
print 'domination number p5[]c8 =', dom(p5.cartesian_product(graphs.CycleGraph(8))), ', ceil(5m/4)=', 5*8/4 print 'domination number p5[]c12 =', dom(p5.cartesian_product(graphs.CycleGraph(12))), ', ceil(5m/4)=', 5*12/4 print 'domination number c5[]p8 =', dom(c5.cartesian_product(graphs.PathGraph(8))), ', ceil(5m/4)=', 5*8/4 print 'domination number c5[]p12 =', dom(c5.cartesian_product(graphs.PathGraph(12))), ', ceil(5m/4)=', 5*12/4 print 'domination number c5[]c8 =', dom(c5.cartesian_product(graphs.CycleGraph(8))), ', ceil(5m/4)=', 5*8/4 print 'domination number c5[]c12 =', dom(c5.cartesian_product(graphs.CycleGraph(12))), ', ceil(5m/4)=', 5*12/4 
       
domination number p5[]c8 = 10 , ceil(5m/4)= 10
domination number p5[]c12 = 15 , ceil(5m/4)= 15
domination number c5[]p8 = 10 , ceil(5m/4)= 10
domination number c5[]p12 = 14 , ceil(5m/4)= 15
domination number c5[]c8 = 10 , ceil(5m/4)= 10
domination number c5[]c12 = 13 , ceil(5m/4)= 15
domination number p5[]c8 = 10 , ceil(5m/4)= 10
domination number p5[]c12 = 15 , ceil(5m/4)= 15
domination number c5[]p8 = 10 , ceil(5m/4)= 10
domination number c5[]p12 = 14 , ceil(5m/4)= 15
domination number c5[]c8 = 10 , ceil(5m/4)= 10
domination number c5[]c12 = 13 , ceil(5m/4)= 15
# p5[]p8: 
       
print 'domination number =', dom(p5.cartesian_product(graphs.PathGraph(8))) print 'throttle product =', thpdx(p5.cartesian_product(graphs.PathGraph(8))) 
       
domination number = 11
throttle product = 11
domination number = 11
throttle product = 11
# n = 6: 
       
print 'domination number p6[]c6 =', dom(p6.cartesian_product(graphs.CycleGraph(6))), ', ceil(6m/4)=', 6*6/4 print 'domination number p6[]c10 =', dom(p6.cartesian_product(graphs.CycleGraph(10))), ', ceil(6m/4)=', 6*10/4 print 'domination number c6[]p6 =', dom(c6.cartesian_product(graphs.PathGraph(6))), ', ceil(6m/4)=', 6*6/4 print 'domination number c6[]p10 =', dom(c6.cartesian_product(graphs.PathGraph(10))), ', ceil(6m/4)=', 6*10/4 print 'domination number c6[]c6 =', dom(c6.cartesian_product(graphs.CycleGraph(6))), ', ceil(6m/4)=', 6*6/4 print 'domination number c6[]c10 =', dom(c6.cartesian_product(graphs.CycleGraph(10))), ', ceil(6m/4)=', 6*10/4 
       
domination number p6[]c6 = 9 , ceil(6m/4)= 9
domination number p6[]c10 = 15 , ceil(6m/4)= 15
domination number c6[]p6 = 9 , ceil(6m/4)= 9
domination number c6[]p10 = 14 , ceil(6m/4)= 15
domination number c6[]c6 = 8 , ceil(6m/4)= 9
domination number c6[]c10 = 14 , ceil(6m/4)= 15
domination number p6[]c6 = 9 , ceil(6m/4)= 9
domination number p6[]c10 = 15 , ceil(6m/4)= 15
domination number c6[]p6 = 9 , ceil(6m/4)= 9
domination number c6[]p10 = 14 , ceil(6m/4)= 15
domination number c6[]c6 = 8 , ceil(6m/4)= 9
domination number c6[]c10 = 14 , ceil(6m/4)= 15
# p6[]p6 
       
print 'domination number =', dom(p6.cartesian_product(graphs.PathGraph(6))) print 'throttle product =', thpdx(p6.cartesian_product(graphs.PathGraph(6))) 
       
domination number = 10
throttle product = 10
domination number = 10
throttle product = 10
# Special cases p5[]p12 and p6[]p10: # # For p5[]p12 and p6[]p10, the domination # is 16 and the size of a minimum power dominating set is 2. 
       
print 'domination number p5[]p12 =', dom(p5.cartesian_product(graphs.PathGraph(12))) print 'domination number p6[]p10 =', dom(p6.cartesian_product(graphs.PathGraph(10))) print 'gamma_P(p5[]p12)=', len(minPowerDominatingSet(p5.cartesian_product(graphs.PathGraph(12)))) print 'gamma_P(p6[]p10)=', len(minPowerDominatingSet(p6.cartesian_product(graphs.PathGraph(10)))) 
       
domination number p5[]p12 = 16
domination number p6[]p10 = 16
gamma_P(p5[]p12)= 2
gamma_P(p6[]p10)= 2
domination number p5[]p12 = 16
domination number p6[]p10 = 16
gamma_P(p5[]p12)= 2
gamma_P(p6[]p10)= 2
# For p5[]p12 and p6[]p10, we need only check that ppt(G,|S|) > t for the following (|S|,t) pairs: (3,5), (5,3) 
       
print 'ppt(p5[]p12,3)=', pptk((graphs.PathGraph(5)).cartesian_product(graphs.PathGraph(12)),3) print 'ppt(p6[]p10,3)=', pptk((graphs.PathGraph(6)).cartesian_product(graphs.PathGraph(10)),3) 
       
ppt(p5[]p12,3)= 6
ppt(p6[]p10,3)= 9
ppt(p5[]p12,3)= 6
ppt(p6[]p10,3)= 9
print 'ppt(p5[]p12,5)=', pptk((graphs.PathGraph(5)).cartesian_product(graphs.PathGraph(12)),5) 
       
ppt(p5[]p12,5)= 4
ppt(p5[]p12,5)= 4
print 'ppt(p6[]p10,5)=', pptk((graphs.PathGraph(6)).cartesian_product(graphs.PathGraph(10)),5) 
       
ppt(p6[]p10,5)= 5
ppt(p6[]p10,5)= 5
#################################### # The function cells are below. # Enter each function cell in order before # trying to run the examples above. #################################### 
       
# F1 # 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/1/.sage/temp/sage.math.iastate.edu/1054/tmp_rHjEmY.pyx...
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling
/home/sageuser/1/.sage/temp/sage.math.iastate.edu/1054/tmp_89lBLg.pyx...
Loading zero_forcing_wavefront.pyx...
Compiling
/home/sageuser/1/.sage/temp/sage.math.iastate.edu/1054/tmp_BvPjca.pyx...
Loading minrank.py...
Loading inertia.py...
Loading Zq_c.pyx...
Compiling /home/sageuser/1/.sage/temp/sage.math.iastate.edu/1054/tmp_rHjEmY.pyx...
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling /home/sageuser/1/.sage/temp/sage.math.iastate.edu/1054/tmp_89lBLg.pyx...
Loading zero_forcing_wavefront.pyx...
Compiling /home/sageuser/1/.sage/temp/sage.math.iastate.edu/1054/tmp_BvPjca.pyx...
Loading minrank.py...
Loading inertia.py...
# F2 # Function for standard 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 
       
# F3 # More functions for standard propagation time # pt(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 
       
# F4 # 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 
       
# F5 # Find power dominating sets of order k # input: a graph G and a and an integer k >= gamma_p(G) # output: printout of all power dominating sets of G of cardinality k def kpds(G,k): V = G.vertices() S = subsets(V,k) for s in S: if isPowerDominatingSet(G,s): print s 
       
# F6 # 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)) 
       
# F7 # compute domination number of graph g def dom(g): gds=g.dominating_set() ds=len(gds) return ds 
       
# F8 # Functions for product power throttling number # # This function computes th_PD^x(G) by minimizing th_PD^x(G,k) # # input: a graph G # output: power product throttling number th_PD^x(G)nd least k such that th_PD^x(G) = th_PD^x(G,k) def thpdx(g): z=PowerDom(g) t=pptk(g,z) thz=z*t kmax=floor(dom(g)/2) for k in [z..kmax]: ptkg=pptk(g,k) if k*ptkg<thz: thz=k*ptkg ko=k thz = min(thz, dom(g)) return thz # This function computes th_PD^x(G) by minimizing th_PD^x(G,k) for # gamma_P(G)<=k <= min(gamma(G),b)/2 # and reports both # th_PD^x(G) and ko such that th_PD^x(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 thpdxb(g,b): d=dom(g) if d<b: b=d z=PowerDom(g) t=pptk(g,z) ko=z thz=z*t kmax=floor(b/2) for k in [z..kmax]: ptkg=pptk(g,k) if k*ptkg<thz: thz=k*ptkg ko=k thz = min(thz, dom(g)) if thz==dom(g): ko=1 print 'th_PD^x(G)', thz, 'k',ko return thz