Receptor Tools#
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'
def print_graph(G):
r"""Prints the vertices and edges of the graph `G`.
"""
print('vertices:',G.vertices(sort=True))
print('edges:',G.edges(sort=True))
mydoc(print_graph)
print_graph(G)
Prints the vertices and edges of the graph G.
def generator(A,rowsum=True):
r"""Creates the generator matrix 'Q' from the adjacency matrix 'A'. By default the sum along each row of 'Q' is zero ('rowsum=True'). For the sum along each column of 'Q' to be zero, set 'rowsum=False'.
"""
if rowsum:
return A-diagonal_matrix(sum(A.T))
else:
return A-diagonal_matrix(sum(A))
mydoc(generator)
generator(A, rowsum=True)
Creates the generator matrix ‘Q’ from the adjacency matrix ‘A’. By default the sum along each row of ‘Q’ is zero (‘rowsum=True’). For the sum along each column of ‘Q’ to be zero, set ‘rowsum=False’.
def hill_diagramatic_method(Q):
r"""Returns the spanning tree polynomials rooted in each vertex."""
n = Q.nrows()
if Q.ncols() != n:
raise ValueError
z = [0]*n
for i in range(n):
a = [ j for j in range(n) ]
a.remove(i)
z[i] = (-1)^(n-1)*Q[a,a].determinant().simplify_full()
return z
mydoc(hill_diagramatic_method)
hill_diagramatic_method(Q)
Returns the spanning tree polynomials rooted in each vertex.
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
def add_edge_monomials(G0,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):
r"""
Add monomials to edges of a graph.
The add_edge_monomials function takes a graph G, as well as 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.
INPUT:
- ``G`` -- graph object (default: `Graph()`);
- ``method`` -- integer (default: ``integer``);
OUTPUT:
- The graph with monomials as edges
"""
G = G0.copy()
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
mydoc(add_edge_monomials)
add_edge_monomials(G0, 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'], ring=False, short_name=False)
Add monomials to edges of a graph.
The add_edge_monomials function takes a graph G, as well as 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.
INPUT:
G– graph object (default:Graph());method– integer (default:integer);
OUTPUT:
The graph with monomials as edges
def add_edge_monomials_ver2(G,edge_labels='default',prefix=''):
r"""
Add monomials to edges of a graph.
[NEED TO WRITE DESCRIPTION]
"""
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
mydoc(add_edge_monomials_ver2)
add_edge_monomials_ver2(G, edge_labels='default', prefix='')
Add monomials to edges of a graph.
[NEED TO WRITE DESCRIPTION]
def enumerate_allosteric_parameters(G=graphs.HouseGraph(),**kwargs):
r"""
Enumerate allosteric parameters of a receptor model.
[NEED TO WRITE DESCRIPTION]
INPUT:
- ``G`` -- graph object (default: `Graph()`);
- ``method`` -- integer (default: ``integer``);
OUTPUT:
- The graph with monomials as vertices
"""
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(graph_border=True,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(graph_border=True,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)
mydoc(enumerate_allosteric_parameters)
enumerate_allosteric_parameters(G=House Graph: Graph on 5 vertices, **kwargs)
Enumerate allosteric parameters of a receptor model.
[NEED TO WRITE DESCRIPTION]
INPUT:
G– graph object (default:Graph());method– integer (default:integer);
OUTPUT:
The graph with monomials as vertices
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)
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,min_ones=1):
# min_ones is the minimum number of ones required in output
import numpy as np
import random
if n_states<=min_ones:
raise ValueError("Variable n_states must be greater than variable min_ones.")
while True:
x=np.random.randint(2, size=n_states)
if sum(x)>=min_ones:
break
return x
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(G, k=2, edge_labels='cannonical'):
r"""
Construct Cartesian power of a graph.
[NEED TO WRITE DESCRIPTION]
"""
# 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_monomials_ver2(Gk,edge_labels='cannonical')
#Gk = add_edge_monomials(Gk)
return Gk
mydoc(cartesian_power)
cartesian_power(G, k=2, edge_labels='cannonical')
Construct Cartesian power of a graph.
[NEED TO WRITE DESCRIPTION]
def reduced_cartesian_power(G, k=2, edge_labels='cannonical',prefix='',independent=false):
r"""
Construct reduced Cartesian power of a graph.
[NEED TO WRITE DESCRIPTION]
"""
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
mydoc(reduced_cartesian_power)
reduced_cartesian_power(G, k=2, edge_labels='cannonical', prefix='', independent=False)
Construct reduced Cartesian power of a graph.
[NEED TO WRITE DESCRIPTION]
def reduced_cartesian_power_ver2(G, k, contexts=True, edge_labels=True):
V=G.vertices(sort=True)
if V[0]!=0: # construct reduced Cartesian power using monomials
f=1
for kk in range(k):
f=f*sum(G.vertices(sort=True))
fmon = f.monomials()
fmon.reverse()
if G.is_directed():
Gk = DiGraph()
else:
Gk = Graph()
for fi in fmon:
for fj in fmon:
df=fj/fi
aj=df.numerator()
ai=df.denominator()
if G.has_edge(ai,aj):
if G.is_directed():
coeff=fi.degree(ai)
else:
coeff=1
if G.edge_label(ai,aj):
# tuple
# Gk.add_edge(fi,fj,(coeff,G.edge_label(ai,aj),fj.gcd(fi)))
# product
if contexts:
Gk.add_edge(fi,fj,coeff*G.edge_label(ai,aj)*fj.gcd(fi))
else:
Gk.add_edge(fi,fj,coeff*G.edge_label(ai,aj))
else:
Gk.add_edge(fi,fj,(coeff,(ai,aj),fj.gcd(fi)))
#print(fi,fj,(ai,aj),fj.gcd(fi))
else: # construct reduced Cartesian power using tuples of integers
Gk = cartesian_power(G,k)
for v in Gk.vertices(sort=True): # do reduction by merging vertices
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])
vcount=[0]*Gk.order()
for i in range(Gk.order()):
vc=[0]*G.order() # yes, G, not Gk, the number of vertices in monomial model
for x in Gk.vertices(sort=True)[i]:
vc[x]=vc[x]+1
vcount[i]=tuple(vc)
Gk.relabel(vcount)
if edge_labels:
for e in Gk.edges(sort=True):
#print(e[0],e[1],coefficient_transition_context_from_tuple(e[0],e[1]))
Gk.set_edge_label(e[0],e[1],coefficient_transition_context_from_tuple(e[0],e[1]))
return Gk
mydoc(reduced_cartesian_power_ver2)
reduced_cartesian_power_ver2(G, k, contexts=True, edge_labels=True)
<IPython.core.display.Markdown object>
def combinatorial_laplacian(G,combinatorial_coefficients=false):
r"""Construct the combinatorial Laplacian of a graph.
[NEED TO WRITE DESCRIPTION]
"""
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
mydoc(combinatorial_laplacian)
combinatorial_laplacian(G, combinatorial_coefficients=False)
Construct the combinatorial Laplacian of a graph.
[NEED TO WRITE DESCRIPTION]
def tree_polynomial(G,combinatorial_coefficients=false):
r"""
Construct the tree polynomial of a graph.
[NEED TO WRITE DESCRIPTION]
"""
L=combinatorial_laplacian(G,combinatorial_coefficients)
T=det(L[1:,1:]).expand()
return T
mydoc(tree_polynomial)
tree_polynomial(G, combinatorial_coefficients=False)
Construct the tree polynomial of a graph.
[NEED TO WRITE DESCRIPTION]
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