# 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
|
|
|
|