Product_Throttling

1137 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 illustrating results in the paper # Survey of product throttling # 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. # Look for "FUNCTIONS-HERE" # # Exampkes are given for both kinds of product throttling for # power domination, standard zero forcing, and PSD forcing. # There is no material for Cops and Robbers. #################################### 
       
#################################### # Section 3 # Illustrate th_pd^*(G) = gamma(G) for several families of graphs #################################### 
       
# cycle g=graphs.CycleGraph(8) print thpda(g), dom(g) 
       
(3, 3) 3
(3, 3) 3
# gridgraph p3=graphs.PathGraph(3) p5=graphs.PathGraph(5) g=p3.cartesian_product(p5) print thpda(g), dom(g) 
       
(4, 4) 4
(4, 4) 4
#################################### # Section 5 # Illustrate th^x(G) = order n of G and th^*(G) = k(G,1) = least k such that pt(G,k)=1 #################################### 
       
# th^x(G) # Since it is known that th^x(G) = order n of G and no set with n/2 < k <n has thZx(G,k) <= n, and thZx(G,n)=n, # thZx(G) is written as a demo and prints lists of k, pt(G,k), and thZx(G,k) for k=Z(G),...,n/2 for comarison # and returns (thZx(G,ko),ko) such that ko is least k having thZx(G)=thZx(G,k)=n 
       
p=graphs.PathGraph(6) thZx(p) 
       
1 5 6
2 2 6
3 1 6
(6, 1)
1 5 6
2 2 6
3 1 6
(6, 1)
g=graphs.CycleGraph(7) thZx(g) 
       
2 3 8
3 2 9
(7, 7)
2 3 8
3 2 9
(7, 7)
g=graphs.CycleGraph(10) thZx(g) 
       
2 4 10
3 4 15
4 2 12
5 2 15
(10, 2)
2 4 10
3 4 15
4 2 12
5 2 15
(10, 2)
g=graphs.CompleteBipartiteGraph(4,5) thZx(g) 
       
(9, 9)
(9, 9)
# th^a(G) # Since it is known that th^*(G) = k(G,1) = least k such that pt(G,k)=1, # thZx(G) is written as a demo and p prints lists of k, pt(G,k), and thZa(G,k) for k=Z(G),..,k(G,1) # and returns (thx(G,ko),ko) such that thZa(G)=thZa(G,ko)=k(G,1) 
       
# Ex 5.8 th^a(P_n)=n/2 for n even p=graphs.PathGraph(10) thZa(p) 
       
1 9 9
2 4 8
3 3 9
4 2 8
5 1 5
(5, 5)
1 9 9
2 4 8
3 3 9
4 2 8
5 1 5
(5, 5)
# th^a(P_n) < th^a(C_n) and th^a(P_n) = th^a(C_n) are both possible 
       
g=graphs.CycleGraph(10) thZa(g) 
       
2 4 8
3 4 12
4 2 8
5 2 10
6 1 6
(6, 6)
2 4 8
3 4 12
4 2 8
5 2 10
6 1 6
(6, 6)
p=graphs.PathGraph(8) thZa(p) 
       
1 7 7
2 3 6
3 2 6
4 1 4
(4, 4)
1 7 7
2 3 6
3 2 6
4 1 4
(4, 4)
g=graphs.CycleGraph(8) thZa(g) 
       
2 3 6
3 3 9
4 1 4
(4, 4)
2 3 6
3 3 9
4 1 4
(4, 4)
# complete graph 
       
g=graphs.CompleteGraph(10) thZa(g) 
       
9 1 9
(9, 9)
9 1 9
(9, 9)
# complete bipartite graph 
       
g=graphs.CompleteBipartiteGraph(5,8) thZa(g) 
       
11 1 11
(11, 11)
11 1 11
(11, 11)
# circulant 
       
g=graphs.CirculantGraph(9,[1,2]) show(g) thZa(g) 
       
4 3 12
5 2 10
6 2 12
7 1 7
(7, 7)
4 3 12
5 2 10
6 2 12
7 1 7
(7, 7)
#################################### # Section 7 # Illustrate th_pd^*(G) and th_pd^x(G) #################################### 
       
# 7.1 th_pd^*(G) for additional families 
       
# hypercube for d in [1..4]: g=graphs.CubeGraph(d) print d, thpda(g), dom(g) 
       
1 (1, 1) 1
2 (2, 2) 2
3 (2, 2) 2
4 (4, 4) 4
1 (1, 1) 1
2 (2, 2) 2
3 (2, 2) 2
4 (4, 4) 4
# necklace k5a=graphs.CompleteGraph(5) k5a.delete_edge([1,4]) k5b=graphs.CompleteGraph(5) k5b.relabel([5,6,7,8,9]) k5b.delete_edge([6,9]) k5c=graphs.CompleteGraph(5) k5c.relabel([10,11,12,13,14]) k5c.delete_edge([11,14]) h=k5a.union(k5b) g=h.union(k5c) g.add_edges([(4,6),(9,11),(14,1)]) show(g) print thpda(g), dom(g) 
       
(3, 3) 3
(3, 3) 3
# G(3,3,4) p6a=graphs.PathGraph(6) p6a.relabel([0,3,4,5,6,1]) p6b=graphs.PathGraph(6) p6b.relabel([0,7,8,9,10,1]) h1=p6a.union(p6b) p6c=graphs.PathGraph(6) p6c.relabel([0,11,12,13,14,1]) h2=h1.union(p6c) p6d=graphs.PathGraph(6) p6d.relabel([1,15,16,17,18,2]) h3=h2.union(p6d) p6e=graphs.PathGraph(6) p6e.relabel([1,19,20,21,22,2]) h4=h3.union(p6e) p6f=graphs.PathGraph(6) p6f.relabel([1,23,24,25,26,2]) h5=h4.union(p6f) p6g=graphs.PathGraph(6) p6g.relabel([2,27,28,29,30,0]) h6=h5.union(p6g) p6h=graphs.PathGraph(6) p6h.relabel([2,31,32,33,34,0]) h7=h6.union(p6h) p6i=graphs.PathGraph(6) p6i.relabel([2,35,36,37,38,0]) g=h7.union(p6i) show(g) 
       
dom(g) 
       
12
12
pdg=PowerDom(g) print pdg pptg=ppt(g) print pptg pdg*pptg 
       
2
5
10
2
5
10
thpda(g) 
       
(6, 3)
(6, 3)
# 7.2 th_pd^x(G) 
       
# cycle g=graphs.CycleGraph(10) thpdx(g) 
       
(6, 1)
(6, 1)
# complete graph g=graphs.CompleteGraph(15) thpdx(g) 
       
(2, 1)
(2, 1)
# complete bipartite graphs 
       
g=graphs.CompleteBipartiteGraph(1,12) thpdx(g) 
       
(2, 1)
(2, 1)
g=graphs.CompleteBipartiteGraph(2,12) thpdx(g) 
       
(3, 1)
(3, 1)
g=graphs.CompleteBipartiteGraph(3,12) thpdx(g) 
       
(4, 2)
(4, 2)
g1=Graph({1:[2,3,4],2:[5,6],6:[7]}) show(g1) thpdx(g1) 
       
(6, 2)
(6, 2)
g2=g1.copy() g2.relabel([8..14]) g3=g1.copy() g3.relabel([15..21]) h=g1.union(g2) g=h.union(g3) g.add_edges([(1,9),(2,8),(8,15)]) show(g) 
       
thpdx(g) 
       
(18, 6)
(18, 6)
# grid graphs 
       
p2=graphs.PathGraph(2) p4=graphs.PathGraph(4) g=p2.cartesian_product(p4) thpdx(g) 
       
(5, 1)
(5, 1)
p2=graphs.PathGraph(2) p6=graphs.PathGraph(6) g=p2.cartesian_product(p6) thpdx(g) 
       
(6, 2)
(6, 2)
#################################### # Section 8 # Illustrate th_+^x(G) and th_+^*(G) #################################### 
       
# 8.1 th_+^x(G) 
       
# th_+^x(T) = 1 + rad T using k = 1 if T is a tree 
       
g=graphs.PathGraph(7) print thpx(g), g.radius() 
       
(4, 1) 3
(4, 1) 3
g=Graph({1:[0,2],3:[0,4],5:[0,6],7:[0,8],9:[0,10],0:[11],12:[11,13],14:[13]}) show(g) print thpx(g), g.radius() 
       
(4, 1) 3
(4, 1) 3
# cycles do interesting things 
       
g=graphs.CycleGraph(6) thpx(g) 
       
(4, 2)
(4, 2)
g=graphs.CycleGraph(11) thpx(g) 
       
(8, 2)
(8, 2)
g=graphs.CycleGraph(15) thpx(g) 
       
(9, 3)
(9, 3)
# Remark 8.4 # analysis of order 4 graphs (computations below) # need to list: 4K_1, K_2 U 2K_1, K_3 U K_1, K_4 # thpx=4 but covered but covered by other case: P_3 U K_1, 2K_2, paw, C_4, diamond # thpx<4:K_1,3, P_3 
       
want=[] for g in graphs.nauty_geng("%s %s:%s "%(4,0,6)): want.append(g.canonical_label().graph6_string()); print len(want) def thpx4(i): g=Graph(want[i]) print thpx(g) show(g) 
       
11
11
thpx4(0) 
       
(4, 4)
(4, 4)
thpx4(1) 
       
(4, 4)
(4, 4)
thpx4(2) 
       
(4, 2)
(4, 2)
thpx4(3) 
       
(2, 1)
(2, 1)
thpx4(4) 
       
(4, 2)
(4, 2)
thpx4(5) 
       
(3, 1)
(3, 1)
thpx4(6) 
       
(4, 4)
(4, 4)
thpx4(7) 
       
(4, 2)
(4, 2)
thpx4(8) 
       
(4, 2)
(4, 2)
thpx4(9) 
       
(4, 2)
(4, 2)
thpx4(10) 
       
(4, 4)
(4, 4)
# 8.2 th_+^a(G) 
       
# for paths and cycles, th_+^a(G) = th_c^a(G) = th_pd^a(G) = gamma(G) = ceiling(n/3) with propagation time 1 
       
g=graphs.CycleGraph(6) thpa(g) 
       
(2, 2)
(2, 2)
g=graphs.CycleGraph(11) thpa(g) 
       
(4, 4)
(4, 4)
g=graphs.CycleGraph(15) thpa(g) 
       
(5, 5)
(5, 5)
# for trees T, th_+^a(T) = th_+^c(T) but th_+^pd(T) may differ 
       
g=Graph({1:[0,2],3:[0,4],5:[0,6],7:[0,8],9:[0,10],0:[11],12:[11,13],14:[13]}) show(g) 
       
thpa(g) 
       
(3, 1)
(3, 1)
thpda(g) 
       
(4, 1)
(4, 1)
#################################### # FUNCTIONS-HERE # These are the main function cells: evaluate each one # before trying to run the examples above. #################################### 
       
#################################### # Standard zero forcing, propagation, and throttling #################################### 
       
# 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/7/.sage/temp/sage.math.iastate.edu/27269/tmp_RN993J.pyx..\
.
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling
/home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_wqq3I9.pyx..\
.
Loading zero_forcing_wavefront.pyx...
Compiling
/home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_rfuftF.pyx..\
.
Loading minrank.py...
Loading inertia.py...
Loading Zq_c.pyx...
Compiling /home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_RN993J.pyx...
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling /home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_wqq3I9.pyx...
Loading zero_forcing_wavefront.pyx...
Compiling /home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_rfuftF.pyx...
Loading minrank.py...
Loading inertia.py...
# Additional function for zero forcing number using brute_force def Z(g): z=len(zero_forcing_set_bruteforce(g)) return z 
       
# Function for standard propagation time of set S, pt(G,S) # adapted by Leslie Hogben from Steve Butler's skew propagation time code # 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 
       
# More functions for standard propagation time # pt(G,k) = min prop time over zero forcing sets of G of size k # 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) savept=len(V)+1 for S in sets: pt=ptz(G,S) if 0 <= pt < savept: savept=pt return savept # propagation time pt(G) # input: a graph G # output: propagation time pt(G)=pt(G,Z(G)) def ptg(g): return ptk(g,Z(g)) 
       
# Function for product-x Z throttling number # # This function computes th^x(G) by minimizing th^x(G,k) # It is known theoretically that this is always the order of the graph. # It contains a print statement to illustrate this # # input: a graph g # output: product throttling number th^x(g) and least k such that th^x(g) = th^x(g,k) def thZx(g): nn=g.order() saveko=Z(g) ptko=ptk(g,saveko) savept=ptko for ko in [saveko..nn/2]: ptko=ptk(g,ko) print ko, ptko, ko*(1+ptko) if ko*(1+ptko)<saveko*(1+savept): saveko=ko savept=ptko if saveko*(1+savept)>nn: saveko=nn savept=0 return saveko*(1+savept), saveko 
       
# Functions for product-* Z throttling number # # This function computes th^*(G) by minimizing th^*(G,k) # It is known theoretically that this is the least ko such that pt(G,k0)=1. # It contains a print statement to illustrate this # # input: a graph G # output: product throttling number th^*(G) and least k such that th^*(G) = th^*(G,k) def thZa(g): ko=Z(g) nn=g.order() ptko=ptk(g,ko) savept=ptko saveko=ko while ptko >1: print ko, ptko, ko*ptko ko=ko+1 ptko=ptk(g,ko) if ko*ptko<saveko*savept: saveko=ko savept=ptko print ko, ptko, ko*ptko return saveko*savept, saveko 
       
#################################### # Power domination, propagation, and throttling #################################### 
       
# 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 
       
# 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)) 
       
# compute domination number of graph g def dom(g): gds=g.dominating_set() ds=len(gds) return ds 
       
# Function for product-x power throttling number # # This function computes th_pd^x(G) by minimizing th_pd^x(G,k) over k in {1,2,...,2/3gamma(G),gamma(G)} # # input: a graph g # output: product throttling number th_pd^x(g) and least k such that th_pd^x(g) = th_pd^x(g,k) def thpdx(g): nn=g.order() saveko=PowerDom(g) savept=pptk(g,saveko) gam=dom(g) for ko in [saveko+1..2*gam/3]: ptko=pptk(g,ko) if ko*(1+ptko)<saveko*(1+savept): saveko=ko savept=ptko if 2*gam<saveko*(1+savept): saveko=gam savept=1 if nn<saveko*(1+savept): saveko=nn savept=0 return saveko*(1+savept), saveko 
       
# Function for product-* power throttling number # # This function computes th_+^*(G) by minimizing th_+^*(G,k) # # input: a graph G # output: product throttling number th_pd^*(G) and k such that th_pd^*(G) = th_pd^*(G,k) # if th_pd^*(G)=gamma(G) then k=gamma(G) is returned; othersie least k is returned def thpda(g): nn=g.order() gam=dom(g) saveko=gam savept=1 kmin=PowerDom(g) for ko in [kmin..gam/2]: ptko=pptk(g,ko) if ko*ptko<saveko*savept: saveko=ko savept=ptko if saveko*savept <= gam/2: return saveko*savept, saveko return saveko*savept, saveko 
       
#################################### # PSD zero forcing, propagation, and throttling #################################### 
       
# 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/7/.sage/temp/sage.math.iastate.edu/27269/tmp_ol3fED.pyx..\
.
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling
/home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_flmr5w.pyx..\
.
Loading zero_forcing_wavefront.pyx...
Compiling
/home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_HAh60W.pyx..\
.
Loading minrank.py...
Loading inertia.py...
xrange test passed
Loading Zq_c.pyx...
Compiling /home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_ol3fED.pyx...
Loading Zq.py...
Loading zero_forcing_64.pyx...
Compiling /home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_flmr5w.pyx...
Loading zero_forcing_wavefront.pyx...
Compiling /home/sageuser/7/.sage/temp/sage.math.iastate.edu/27269/tmp_HAh60W.pyx...
Loading minrank.py...
Loading inertia.py...
# 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 
       
# Function for product-x PSD throttling number # # This function computes th_+^x(G) by minimizing th_+^x(G,k) # # input: a graph g # output: product throttling number th^x(g) and least k such that th^x(g) = th^x(g,k) def thpx(g): nn=g.order() saveko=Zplus(g) savept=ptpk(g,saveko) for ko in [saveko+1..nn/2]: ptko=ptpk(g,ko) if ko*(1+ptko)<saveko*(1+savept): saveko=ko savept=ptko if nn<saveko*(1+savept): saveko=nn savept=0 return saveko*(1+savept), saveko 
       
# Function for product-* PSD throttling number # # This function computes th_+^*(G) by minimizing th_+^*(G,k) # # input: a graph G # output: product throttling number th_+^*(G) and least k such that th_+^*(G) = th_+^*(G,k) def thpa(g): ko=Zplus(g) nn=g.order() ptko=ptpk(g,ko) savept=ptko saveko=ko while ptko >1: ko=ko+1 ptko=ptpk(g,ko) if ko*ptko<saveko*savept: saveko=ko savept=ptko return saveko*savept, saveko