My Tools#
mydoc(myfun,method='both')
def mydoc(myfun,method='both'):
r"""
Display the **signature** of `myfun` followed by the **docstring** of `myfun`.
Set `method ='signature'` or `'method = 'docstring'` or `'method = 'both'`
"""
import inspect
from IPython.display import display, Markdown, Latex, display_markdown
if method == 'signature' or method == 'both':
from inspect import signature
display(Markdown("**`"+str(myfun.__name__)+str(signature(myfun))+"`**"))
if method == 'docstring' or method == 'both':
display(Markdown(inspect.getdoc(myfun)))
mydoc(myfun, method='both')
Display the signature of myfun followed by the docstring of myfun.
Set method ='signature' or 'method = 'docstring' or 'method = 'both'
my_graph_show(G)#
def my_graph_show(G):
r"""Show a graph Greg's way.
Prints the vertices and edges of the graph `G`.
"""
print('vertices:',G.vertices(sort=True))
print('edges:',G.edges(sort=True))
my_graph_show(G)
Show a graph Greg’s way.
Prints the vertices and edges of the graph G.
add_vertex_monomials(G=Graph(),method='integer',ring=False)
add_vertex_monomials(G=Graph on 0 vertices, method='integer', ring=False)
Add monomials to vertices of a graph.
The add_vertex_monomials function takes a graph G, as well as optional parameters method and ring. The function creates a new graph H with vertices labeled by monomials. The monomials are chosen based on the number of vertices in G. If the method parameter is set to ‘alpha’ and the number of vertices in G is less than or equal to 10, the monomials are chosen as alphabetical letters (‘a’ to ‘k’). Otherwise, the monomials are chosen as strings of the form ‘a0’, ‘a1’, …, ‘an-1’, where n is the number of vertices in G. The function then adds the vertices from G to H using the monomials as labels, and adds the edges from G to H using the monomials as endpoints. If the ring parameter is set to True, the function also creates a polynomial ring V with the chosen monomials and ‘invlex’ order, and returns both H and V. Otherwise, it returns only H.
INPUT:
G– graph object (default:Graph());method– integer (default:integer);
OUTPUT: The graph with monomials as vertices
EXAMPLES:
This example illustrates … ::
sage: A = ModuliSpace()
sage: A.point(2,3)
xxx
def replace_powers(s,nchar=1):
result=''
i=0
while i<len(s):
if s[i:i+2]=='^2':
result+=s[i-nchar:i]
i+=2
else:
result+=s[i]
i+=1
return result
def table_multiply(Alpha,Beta):
AlphaBeta=[]
for A,B in zip(Alpha,Beta):
AlphaBeta.append([a*b for a,b in zip(A,B)])
return AlphaBeta
def repmat(value,n,*args):
# n rows and m columns
if args:
m=args[0]
else:
m=n
return [[value] * m for _ in range(n)]
def zeros(n,*args):
if args:
m=args[0]
else:
m=n
return repmat(0,n,m)
def ones(n,*args):
if args:
m=args[0]
else:
m=n
return repmat(1,n,m)
The function add_edge_monomials takes in a graph object G and optional parameters method, edge_vars, ring, and short_name.
If method is set to ‘integer’, the function creates a polynomial ring using the given edge variables and assigns variables to the edges of the graph. The edge variables can be represented either as ‘e’ followed by the first vertex label or the first and second vertex labels concatenated. If the vertex labels are integers and the short_name parameter is set to True, the edge variables are created using only the first vertex label.
If method is set to ‘alpha’, the function creates a polynomial ring using the given edge variables and assigns variables to the edges of the graph in reverse order. The number of edge variables used is determined by the size of the graph.
The ring parameter, if set to True, injects the polynomial variables into the global namespace and returns the graph and the polynomial ring. Otherwise, it simply returns the graph.
def add_edge_monomials(G,method='integer',
edge_vars=['b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z'], # no "a" here!
ring=False,short_name=False):
if method=='integer': # method parameter is set to 'integer'
if G.vertices(sort=True)[0]==0: # checking if the vertex labels are integers
if short_name:
edge_vars = ['e%s' %e[0] for e in G.edges(sort=True)] # create a list of edge variables using the first vertex label
else:
edge_vars = ['e%s%s' %(e[0],e[1]) for e in G.edges(sort=True)] # create a list of edge variables using both vertex labels
E=PolynomialRing(ZZ,names=edge_vars,order='invlex') # create a polynomial ring using the edge variables
else:
if short_name:
edge_vars = ['e%s' %e[0] for e in G.edges(sort=True)] # create a list of edge variables using the first vertex label
edge_vars = [ ev.replace('a','') for ev in edge_vars ] # remove 'a' from the edge variables
edge_subscripts = ['%s' %(e[0],e[1]) for e in G.edges(sort=True)] # create a list of edge subscripts using both vertex labels
edge_subscripts = [ ss.replace('a','') for ss in edge_subscripts ] # remove 'a' from edge subscripts
else:
edge_vars = ['e_%s%s' %(e[0],e[1]) for e in G.edges(sort=True)] # create a list of edge variables using both vertex labels
edge_subscripts = ['%s%s' %(e[0],e[1]) for e in G.edges(sort=True)] # create a list of edge subscripts using both vertex labels
edge_names = []
for sub in edge_subscripts:
edge_names.append('e_{\\mathit{' + sub.replace('e','e_') + '}}') # create a list of latex names for the edge variables
E=PolynomialRing(ZZ,names=edge_vars,order='invlex') # create a polynomial ring using the edge variables
E._latex_names = edge_names # assign latex names to the polynomial variables
gen=0;
for e in G.edges(sort=True):
G.add_edge(e[0],e[1],E.gen(gen)) # add an edge with the corresponding variable to the graph
gen=gen+1
else: # method parameter is set to 'alpha'
edge_vars = edge_vars[:G.size()] # get the first G.size() edge variables
E=PolynomialRing(ZZ,names=edge_vars,order='invlex') # create a polynomial ring using the edge variables
evars = list(E.gens()[:G.size()]) # create a list of the polynomial variables
evars.reverse()
for e in G.edges(sort=True):
G.add_edge(e[0],e[1],evars.pop()) # add an edge with the corresponding variable to the graph
if ring:
# E.inject_variables() # inject the polynomial variables into the global namespace
return (G, E) # return the graph and the polynomial ring
else:
return G # return the graph
def enumerate_allosteric_parameters(G=graphs.HouseGraph(),**kwargs):
defaultKwargs = { 'method': 'integer', 'verbose': False, 'show': False}
kwargs = { **defaultKwargs, **kwargs }
#print(kwargs)
my_method = kwargs['method']
if my_method not in ['integer','alpha']:
raise ValueError("enumerate_allosteric_parameters says: method must be 'integer' or 'alpha'.")
verbose = kwargs['verbose']
myshow = kwargs['show']
if my_method=='alpha':
(G, V) = add_vertex_monomials(G,method='alpha',ring=True)
Groot=V.gen(0)
else:
Groot=0
if verbose or myshow:
G.show(edge_labels=False)
(BFSVertexList,BFSTree) = G.lex_BFS(tree=True,initial_vertex=Groot)
(T,E)=add_edge_monomials(BFSTree,method=my_method,ring=True,short_name=True)
T.set_pos(G.get_pos())
if verbose or myshow:
T.show(edge_labels=True)
if verbose:
show(E)
F = [E(0)]
for v in BFSVertexList[1:]: # every element except the first (the root 0, which has been set to E(0))
P = T.all_paths(v, BFSVertexList[0], use_multiedges=False, report_edges=True, labels=True)
P = P[0] # there is only one path, so take the first
f=0
for p in P:
f = f+p[2]
F.append(f)
nf = len(F) # this should be equal to G.order()
if verbose:
show(F)
KappaVars = []; KappaNames = []
for e in E.gens():
estr = str(e)
if len(estr)!=1:
estr=estr.replace('e','')
KappaVars.append('kappa_'+estr)
KappaNames.append('\\kappa_{\\mathit{' + estr + '}}')
if verbose:
print(KappaVars)
print(KappaNames)
EtaProb = zeros(nf); EtaMons = zeros(nf); EtaCoef = zeros(nf)
for i0 in range(len(F)):
for i1 in range(len(F)):
fprod=F[i0]*F[i1]
EtaProb[i0][i1]=fprod
EtaMons[i0][i1]=fprod.monomials()
EtaCoef[i0][i1]=fprod.coefficients()
if verbose:
show(table(EtaProb))
show(table(EtaMons))
show(table(EtaCoef))
UniqueMons = []
for mon in sorted(set(flatten(EtaMons))):
UniqueMons.append(mon)
if verbose:
show(UniqueMons)
EtaVars = []; EtaNames = []
for mon in UniqueMons:
vvar = str(mon)
vvar = vvar.replace('*', '')
if my_method=='alpha':
nchar=1
else:
nchar=2
vvar = replace_powers(vvar,nchar)
if my_method=='integer':
vvar = vvar.replace('e','')
EtaVars.append('eta_'+vvar)
EtaNames.append('\\eta_{\\mathit{' + vvar + '}}')
if verbose:
print(EtaVars)
print(EtaNames)
A=PolynomialRing(ZZ,names=KappaVars+EtaVars,order='invlex')
# A.inject_variables()
A._latex_names = KappaNames+EtaNames
# https://ask.sagemath.org/question/8202/how-to-give-latex-names-to-generators-of-polynomial-rings/
d_vars=dict(zip(list(E.gens())+UniqueMons,list(A.gens())))
d_vars[E(1)]=1
if verbose:
print(d_vars)
KappaProb = zeros(nf); KappaMons = zeros(nf); KappaCoef = zeros(nf);
Kappa = repmat(A(1),nf)
for i0 in range(len(F)):
for i1 in range(0,i0):
Kappa[i0][i1]=0
for i1 in range(i0,len(F)):
fsum = F[i0]+F[i1]
KappaProb[i0][i1]=fsum
KappaMons[i0][i1]=fsum.monomials()
KappaCoef[i0][i1]=fsum.coefficients()
for m in range(len(KappaMons[i0][i1])):
mon = KappaMons[i0][i1][m]
Kappa[i0][i1]*= d_vars[mon]^KappaCoef[i0][i1][m]
if i0!=i1:
Kappa[i0][i1]*=2
if verbose:
print(KappaProb)
print(KappaMons)
print(KappaCoef)
print(Kappa)
Eta = repmat(A(1),nf)
for i0 in range(len(F)):
for i1 in range(0, i0):
Eta[i0][i1] = 0
for i1 in range(i0, len(F)):
for m in range(len(EtaMons[i0][i1])):
mon = EtaMons[i0][i1][m]
Eta[i0][i1] *= d_vars[mon] ** EtaCoef[i0][i1][m]
if verbose:
print(Eta)
KappaEta = table_multiply(Kappa,Eta)
if verbose or myshow:
show(table(KappaEta))
return (G, T, KappaEta, A)
def normalize(phi):
import numpy as np
phi_sum = np.sum(phi)
return np.array([np.divide(p,phi_sum) for p in phi])
def random_binary_array(n_states):
import numpy as np
import random
return np.random.randint(2, size=n_states)
def make_random_q(n_states):
import numpy as np
import random
q = np.random.uniform(low=0,high=10,size=n_states)
q[0]=1
return q
def pade(phi_list,q_list,x_list):
import numpy as np
num = np.zeros(x_list.size)
den = np.zeros(x_list.size)
for i, (phi, q) in enumerate(zip(phi_list,q_list)):
num = np.add(num,[phi*q*x**i for x in x_list])
den = np.add(den,[q*x**i for x in x_list])
return np.divide(num,den)
def make_simulated_data(phi,q,xlogmin=-3,xlogmax=3,xnum=20,noise=0):
import numpy as np
import random
if len(q)!=len(phi):
raise ValueError("make_simulated_data says: phi and q are not the same length.")
if q[0]!=1:
print('Warning: Setting q0 to 1.')
q[0]=1
x = np.logspace(xlogmin,xlogmax,xnum)
y = pade(phi,q,x) + np.random.normal(size=x.size, scale=noise)
return (x,y)
def dimerize(phi):
n = len(phi)
Phi = zeros(n)
for i0 in range(n):
for i1 in range(i0,n):
Phi[i0][i1] = phi[i0]*phi[i1]
if i0==i1:
Phi[i0][i1]*=2
return Phi
def make_symbolic_dimer_binding_curve(Phi,KappaEta):
n=len(Phi)
num=0; den = 0;
for i0 in range(n):
for i1 in range(i0,n):
num += Phi[i0][i1]*KappaEta[i0][i1]
den += KappaEta[i0][i1]
return num/den
def transition_and_context(f,t):
# f=from tuple, t=to tuple
if len(f)!=len(t):
raise ValueError('Tuples must be the same size.')
f=sorted(list(f));
t=sorted(list(t));
def remove_common_item():
for tt in t:
for ff in f:
if ff==tt:
f.remove(ff)
t.remove(ff)
return ff
return []
context=[];
while 1:
c=remove_common_item()
if c==[]:
return f[0],t[0],context
else:
context.append(c)
end
def combinatorial_coefficient(f,t):
# f=from tuple, t=to tuple
ff,tt,context=transition_and_context(f,t)
coeff=1
for c in context:
if c==ff:
coeff=coeff+1;
return coeff
def Cartesian_power_original(G, k, edge_labels='cannonical'):
# Make Cartesian power G^k (unreduced)
Gk=G
for i in range(k-1):
Gk = Gk.cartesian_product(G)
# Make each vertex a tuple
vflat=list(range(Gk.order()));
for i in range(Gk.order()):
v=Gk.vertices(sort=True)[i]
#print(i,v,flatten(v))
vflat[i]=tuple(flatten(v))
Gk.relabel(vflat)
if edge_labels=='cannonical':
Gk = add_edge_labels_ver2(Gk,edge_labels='cannonical')
#Gk = add_edge_monomials(Gk)
return Gk
def reduced_Cartesian_power_original(G, k, edge_labels='cannonical',prefix='',independent=false):
Gk = Cartesian_power(G, k, edge_labels)
for v in Gk.vertices(sort=True):
for u in Gk.vertices(sort=True):
sv=tuple(sorted(v))
su=tuple(sorted(u))
if v!=u and sv==su and v<u and Gk.has_vertex(v) and Gk.has_vertex(u):
Gk.merge_vertices([v,u])
if edge_labels=='induced':
for e in Gk.edges(sort=True):
fr,to,context=transition_and_context(e[0],e[1])
if prefix=='':
Gk.set_edge_label(e[0],e[1],'(%s,%s)%s' %(fr,to,context))
else:
if independent:
Gk.set_edge_label(e[0],e[1],prefix+'%s%s' %(fr,to))
else:
ctxt=''
for k in context:
if G.order() <= 9:
ctxt=ctxt+str(k)
else:
ctxt=ctxt+'_'+str(k)
Gk.set_edge_label(e[0],e[1],prefix+'%s%s_%s' %(fr,to,ctxt))
return Gk
def add_edge_labels_ver2(G,edge_labels='default',prefix=''):
if edge_labels=='default':
if G.to_simple().size() <= 26:
edge_labels=['a','b','c','d','e','f','g','h','i','j','k','l','m','n','o','p','q','r','s','t','u','v','w','x','y','z']
else:
edge_labels=['e%s' %(i) for i in range(G.to_simple().size())]
remaining_edge_labels=edge_labels
for e in G.edges(sort=True):
if edge_labels=='cannonical':
if G.order() <= 9:
if prefix=='':
G.set_edge_label(e[0],e[1],'(%s,%s)' %(e[0],e[1]))
else:
G.set_edge_label(e[0],e[1],'%s%s%s' %(prefix,e[0],e[1]))
else:
if prefix=='':
G.set_edge_label(e[0],e[1],'(%s,%s)' %(e[0],e[1]))
else:
G.set_edge_label(e[0],e[1],'%s_%s_%s' %(prefix,e[0],e[1]))
else:
if G.is_directed():
if e[0]<e[1]:
edge_label=remaining_edge_labels.pop(0)
G.set_edge_label(e[0],e[1],edge_label)
G.set_edge_label(e[1],e[0],edge_label+'bar')
else:
edge_label=remaining_edge_labels.pop(0)
G.set_edge_label(e[0],e[1],edge_label)
return G
def combinatorial_laplacian(G,combinatorial_coefficients=false):
v=G.vertices(sort=True)
A = matrix(SR, G.order(), lambda i,j: G.edge_label(v[i],v[j]) if G.has_edge(v[i],v[j]) else 0)
if combinatorial_coefficients:
for i in range(G.order()):
for j in range(G.order()):
if G.has_edge(v[i],v[j]):
A[i,j]=combinatorial_coefficient(v[i],v[j])*A[i,j]
D = diagonal_matrix(sum(A.transpose()))
L=D-A
return L
def tree_polynomial(G,combinatorial_coefficients=false):
L=combinatorial_laplacian(G,combinatorial_coefficients)
T=det(L[1:,1:]).expand()
return T
def symmetric_directed(G):
D = G.to_directed()
for e in D.edges(sort=True):
if e[2]!=None and e[0]>e[1]:
D.set_edge_label(e[0],e[1],e[2]+'bar')
return D