Cycle fluxes in receptor models#
UNDER CONSTRUCTION#
%%capture
%run receptor_tools.ipynb
As discussed earlier in Markov chain tree theorem, the Markov chain tree theorem (also called Hill’s diagrammatic method) provides algebraic expressions for the steady-state probabilities of a receptor model. The function hill_diagrammatic_method() below takes a symbolic generator matrix for such a Markov chain as input and, using the Markov chain tree theorem, calculates the steady-state probability of each state.
The receptor models considered here have state-transition diagrams with the topology of a simple symmetric directed graph (connected, no loops). The Markov chain tree theorem always applies to receptor models with these properties. However, there are Markov chains that cannot be analyzed in this way (e.g., Markov chains with an infinite number of states).
Receptor model with topology of the House Graph#
Consider receptor model that has more than one cycle, namely, the house graph.
var('p0 p1 p2 p3 p4 a01 a10 a02 a20 a04 a40 a12 a21 a23 a32 a34 a43')
d = {0: {1:a01, 2:a02, 4:a04}, 1: {0:a10, 2:a12}, 2: {1:a21, 0:a20, 3:a23}, 3: {2:a32, 4:a34}, 4: {3:a43, 0:a40}};
G = DiGraph(d,weighted=True)
vertex_positions = {0: (0, 0), 1: (1, 1.41), 2: (2, 0), 3: (2,-2), 4: (0,-2)}
G.plot(figsize=8,edge_labels=True,pos=vertex_positions,graph_border=True)
A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.2.6 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.
If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.
Traceback (most recent call last): File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
return _run_code(code, main_globals, None,
File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
exec(code, run_globals)
File "/usr/lib/python3/dist-packages/sage/repl/ipython_kernel/__main__.py", line 3, in <module>
IPKernelApp.launch_instance(kernel_class=SageKernel)
File "/usr/lib/python3/dist-packages/traitlets/config/application.py", line 846, in launch_instance
app.start()
File "/usr/lib/python3/dist-packages/ipykernel/kernelapp.py", line 677, in start
self.io_loop.start()
File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 199, in start
self.asyncio_loop.run_forever()
File "/usr/lib/python3.10/asyncio/base_events.py", line 603, in run_forever
self._run_once()
File "/usr/lib/python3.10/asyncio/base_events.py", line 1909, in _run_once
handle._run()
File "/usr/lib/python3.10/asyncio/events.py", line 80, in _run
self._context.run(self._callback, *self._args)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 461, in dispatch_queue
await self.process_one()
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 450, in process_one
await dispatch(*args)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 357, in dispatch_shell
await result
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 652, in execute_request
reply_content = await reply_content
File "/usr/lib/python3/dist-packages/ipykernel/ipkernel.py", line 353, in do_execute
res = shell.run_cell(code, store_history=store_history, silent=silent)
File "/usr/lib/python3/dist-packages/ipykernel/zmqshell.py", line 532, in run_cell
return super().run_cell(*args, **kwargs)
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2914, in run_cell
result = self._run_cell(
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2960, in _run_cell
return runner(coro)
File "/usr/lib/python3/dist-packages/IPython/core/async_helpers.py", line 78, in _pseudo_sync_runner
coro.send(None)
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3185, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3377, in run_ast_nodes
if (await self.run_code(code, result, async_=asy)):
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3457, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "/tmp/ipykernel_11042/3441472480.py", line 5, in <module>
G.plot(figsize=Integer(8),edge_labels=True,pos=vertex_positions,graph_border=True)
File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 262, in __call__
format_dict, md_dict = self.compute_format_data(result)
File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 151, in compute_format_data
return self.shell.display_formatter.format(result)
File "/usr/lib/python3/dist-packages/sage/repl/display/formatter.py", line 181, in format
sage_format, sage_metadata = self.dm.displayhook(obj)
File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 825, in displayhook
plain_text, rich_output = self._rich_output_formatter(obj, dict())
File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 643, in _rich_output_formatter
rich_output = self._call_rich_repr(obj, rich_repr_kwds)
File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 603, in _call_rich_repr
return obj._rich_repr_(self)
File "/usr/lib/python3/dist-packages/sage/plot/graphics.py", line 1000, in _rich_repr_
return display_manager.graphics_from_save(
File "/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py", line 731, in graphics_from_save
save_function(filename, **kwds)
File "/usr/lib/python3/dist-packages/sage/misc/decorators.py", line 410, in wrapper
return func(*args, **kwds)
File "/usr/lib/python3/dist-packages/sage/plot/graphics.py", line 3296, in save
from matplotlib import rcParams
File "/usr/lib/python3/dist-packages/matplotlib/__init__.py", line 109, in <module>
from . import _api, _version, cbook, docstring, rcsetup
File "/usr/lib/python3/dist-packages/matplotlib/rcsetup.py", line 27, in <module>
from matplotlib.colors import Colormap, is_color_like
File "/usr/lib/python3/dist-packages/matplotlib/colors.py", line 56, in <module>
from matplotlib import _api, cbook, scale
File "/usr/lib/python3/dist-packages/matplotlib/scale.py", line 23, in <module>
from matplotlib.ticker import (
File "/usr/lib/python3/dist-packages/matplotlib/ticker.py", line 136, in <module>
from matplotlib import transforms as mtransforms
File "/usr/lib/python3/dist-packages/matplotlib/transforms.py", line 46, in <module>
from matplotlib._path import (
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
AttributeError: _ARRAY_API not found
/usr/lib/python3/dist-packages/sage/repl/rich_output/display_manager.py:608: RichReprWarning: Exception in _rich_repr_ while displaying object: numpy.core.multiarray failed to import
warnings.warn(
A module that was compiled using NumPy 1.x cannot be run in
NumPy 2.2.6 as it may crash. To support both 1.x and 2.x
versions of NumPy, modules must be compiled with NumPy 2.0.
Some module may need to rebuild instead e.g. with 'pybind11>=2.12'.
If you are a user of the module, the easiest solution will be to
downgrade to 'numpy<2' or try to upgrade the affected module.
We expect that some modules will need time to support NumPy 2.
Traceback (most recent call last): File "/usr/lib/python3.10/runpy.py", line 196, in _run_module_as_main
return _run_code(code, main_globals, None,
File "/usr/lib/python3.10/runpy.py", line 86, in _run_code
exec(code, run_globals)
File "/usr/lib/python3/dist-packages/sage/repl/ipython_kernel/__main__.py", line 3, in <module>
IPKernelApp.launch_instance(kernel_class=SageKernel)
File "/usr/lib/python3/dist-packages/traitlets/config/application.py", line 846, in launch_instance
app.start()
File "/usr/lib/python3/dist-packages/ipykernel/kernelapp.py", line 677, in start
self.io_loop.start()
File "/usr/lib/python3/dist-packages/tornado/platform/asyncio.py", line 199, in start
self.asyncio_loop.run_forever()
File "/usr/lib/python3.10/asyncio/base_events.py", line 603, in run_forever
self._run_once()
File "/usr/lib/python3.10/asyncio/base_events.py", line 1909, in _run_once
handle._run()
File "/usr/lib/python3.10/asyncio/events.py", line 80, in _run
self._context.run(self._callback, *self._args)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 461, in dispatch_queue
await self.process_one()
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 450, in process_one
await dispatch(*args)
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 357, in dispatch_shell
await result
File "/usr/lib/python3/dist-packages/ipykernel/kernelbase.py", line 652, in execute_request
reply_content = await reply_content
File "/usr/lib/python3/dist-packages/ipykernel/ipkernel.py", line 353, in do_execute
res = shell.run_cell(code, store_history=store_history, silent=silent)
File "/usr/lib/python3/dist-packages/ipykernel/zmqshell.py", line 532, in run_cell
return super().run_cell(*args, **kwargs)
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2914, in run_cell
result = self._run_cell(
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 2960, in _run_cell
return runner(coro)
File "/usr/lib/python3/dist-packages/IPython/core/async_helpers.py", line 78, in _pseudo_sync_runner
coro.send(None)
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3185, in run_cell_async
has_raised = await self.run_ast_nodes(code_ast.body, cell_name,
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3377, in run_ast_nodes
if (await self.run_code(code, result, async_=asy)):
File "/usr/lib/python3/dist-packages/IPython/core/interactiveshell.py", line 3457, in run_code
exec(code_obj, self.user_global_ns, self.user_ns)
File "/tmp/ipykernel_11042/3441472480.py", line 5, in <module>
G.plot(figsize=Integer(8),edge_labels=True,pos=vertex_positions,graph_border=True)
File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 262, in __call__
format_dict, md_dict = self.compute_format_data(result)
File "/usr/lib/python3/dist-packages/IPython/core/displayhook.py", line 151, in compute_format_data
return self.shell.display_formatter.format(result)
File "/usr/lib/python3/dist-packages/sage/repl/display/formatter.py", line 186, in format
if (not isinstance(obj, (IPYTHON_NATIVE_TYPES, Figure)) and
File "/usr/lib/python3/dist-packages/matplotlib/__init__.py", line 109, in <module>
from . import _api, _version, cbook, docstring, rcsetup
File "/usr/lib/python3/dist-packages/matplotlib/rcsetup.py", line 27, in <module>
from matplotlib.colors import Colormap, is_color_like
File "/usr/lib/python3/dist-packages/matplotlib/colors.py", line 56, in <module>
from matplotlib import _api, cbook, scale
File "/usr/lib/python3/dist-packages/matplotlib/scale.py", line 23, in <module>
from matplotlib.ticker import (
File "/usr/lib/python3/dist-packages/matplotlib/ticker.py", line 136, in <module>
from matplotlib import transforms as mtransforms
File "/usr/lib/python3/dist-packages/matplotlib/transforms.py", line 46, in <module>
from matplotlib._path import (
---------------------------------------------------------------------------
AttributeError Traceback (most recent call last)
AttributeError: _ARRAY_API not found
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
/tmp/ipykernel_11042/3441472480.py in <module>
3 G = DiGraph(d,weighted=True)
4 vertex_positions = {Integer(0): (Integer(0), Integer(0)), Integer(1): (Integer(1), RealNumber('1.41')), Integer(2): (Integer(2), Integer(0)), Integer(3): (Integer(2),-Integer(2)), Integer(4): (Integer(0),-Integer(2))}
----> 5 G.plot(figsize=Integer(8),edge_labels=True,pos=vertex_positions,graph_border=True)
/usr/lib/python3/dist-packages/IPython/core/displayhook.py in __call__(self, result)
260 self.start_displayhook()
261 self.write_output_prompt()
--> 262 format_dict, md_dict = self.compute_format_data(result)
263 self.update_user_ns(result)
264 self.fill_exec_result(result)
/usr/lib/python3/dist-packages/IPython/core/displayhook.py in compute_format_data(self, result)
149
150 """
--> 151 return self.shell.display_formatter.format(result)
152
153 # This can be set to True by the write_output_prompt method in a subclass
/usr/lib/python3/dist-packages/sage/repl/display/formatter.py in format(self, obj, include, exclude)
184 # use Sage rich output for any except those native to IPython, but only
185 # if it is not plain and dull
--> 186 if (not isinstance(obj, (IPYTHON_NATIVE_TYPES, Figure)) and
187 not set(sage_format.keys()).issubset([PLAIN_TEXT])):
188 return sage_format, sage_metadata
/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport.__instancecheck__ (build/cythonized/sage/misc/lazy_import.c:7695)()
914 True
915 """
--> 916 return isinstance(x, self.get_object())
917
918 def __subclasscheck__(self, x):
/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport.get_object (build/cythonized/sage/misc/lazy_import.c:2612)()
215 if likely(self._object is not None):
216 return self._object
--> 217 return self._get_object()
218
219 cpdef _get_object(self):
/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport._get_object (build/cythonized/sage/misc/lazy_import.c:3073)()
255 if self._feature:
256 raise FeatureNotPresentError(self._feature, reason=f'Importing {self._name} failed: {e}')
--> 257 raise
258
259 name = self._as_name
/usr/lib/python3/dist-packages/sage/misc/lazy_import.pyx in sage.misc.lazy_import.LazyImport._get_object (build/cythonized/sage/misc/lazy_import.c:2935)()
251
252 try:
--> 253 self._object = getattr(__import__(self._module, {}, {}, [self._name]), self._name)
254 except ImportError as e:
255 if self._feature:
/usr/lib/python3/dist-packages/matplotlib/__init__.py in <module>
107 # cbook must import matplotlib only within function
108 # definitions, so it is safe to import from it here.
--> 109 from . import _api, _version, cbook, docstring, rcsetup
110 from matplotlib.cbook import MatplotlibDeprecationWarning, sanitize_sequence
111 from matplotlib.cbook import mplDeprecation # deprecated
/usr/lib/python3/dist-packages/matplotlib/rcsetup.py in <module>
25 from matplotlib import _api, cbook
26 from matplotlib.cbook import ls_mapper
---> 27 from matplotlib.colors import Colormap, is_color_like
28 from matplotlib.fontconfig_pattern import parse_fontconfig_pattern
29 from matplotlib._enums import JoinStyle, CapStyle
/usr/lib/python3/dist-packages/matplotlib/colors.py in <module>
54 import matplotlib as mpl
55 import numpy as np
---> 56 from matplotlib import _api, cbook, scale
57 from ._color_data import BASE_COLORS, TABLEAU_COLORS, CSS4_COLORS, XKCD_COLORS
58
/usr/lib/python3/dist-packages/matplotlib/scale.py in <module>
21 import matplotlib as mpl
22 from matplotlib import _api, docstring
---> 23 from matplotlib.ticker import (
24 NullFormatter, ScalarFormatter, LogFormatterSciNotation, LogitFormatter,
25 NullLocator, LogLocator, AutoLocator, AutoMinorLocator,
/usr/lib/python3/dist-packages/matplotlib/ticker.py in <module>
134 import matplotlib as mpl
135 from matplotlib import _api, cbook
--> 136 from matplotlib import transforms as mtransforms
137
138 _log = logging.getLogger(__name__)
/usr/lib/python3/dist-packages/matplotlib/transforms.py in <module>
44
45 from matplotlib import _api
---> 46 from matplotlib._path import (
47 affine_transform, count_bboxes_overlapping_bbox, update_path_extents)
48 from .path import Path
ImportError: numpy.core.multiarray failed to import
Q = generator(G.weighted_adjacency_matrix(sparse=False))
z = hill_diagrammatic_method(Q)
for i in range(5):
print('z[%s] ='%i, f'{expand(z[i])}')
z[0] = a10*a20*a32*a40 + a12*a20*a32*a40 + a10*a21*a32*a40 + a10*a20*a34*a40 + a12*a20*a34*a40 + a10*a21*a34*a40 + a10*a23*a34*a40 + a12*a23*a34*a40 + a10*a20*a32*a43 + a12*a20*a32*a43 + a10*a21*a32*a43
z[1] = a01*a20*a32*a40 + a01*a21*a32*a40 + a02*a21*a32*a40 + a01*a20*a34*a40 + a01*a21*a34*a40 + a02*a21*a34*a40 + a01*a23*a34*a40 + a01*a20*a32*a43 + a01*a21*a32*a43 + a02*a21*a32*a43 + a04*a21*a32*a43
z[2] = a02*a10*a32*a40 + a01*a12*a32*a40 + a02*a12*a32*a40 + a02*a10*a34*a40 + a01*a12*a34*a40 + a02*a12*a34*a40 + a02*a10*a32*a43 + a04*a10*a32*a43 + a01*a12*a32*a43 + a02*a12*a32*a43 + a04*a12*a32*a43
z[3] = a02*a10*a23*a40 + a01*a12*a23*a40 + a02*a12*a23*a40 + a04*a10*a20*a43 + a04*a12*a20*a43 + a04*a10*a21*a43 + a02*a10*a23*a43 + a04*a10*a23*a43 + a01*a12*a23*a43 + a02*a12*a23*a43 + a04*a12*a23*a43
z[4] = a04*a10*a20*a32 + a04*a12*a20*a32 + a04*a10*a21*a32 + a04*a10*a20*a34 + a04*a12*a20*a34 + a04*a10*a21*a34 + a02*a10*a23*a34 + a04*a10*a23*a34 + a01*a12*a23*a34 + a02*a12*a23*a34 + a04*a12*a23*a34
In this case, each of the five relative probabilities is a multinomial with 10 terms because there are 10 rooted spanning trees associated to each state (vertex). Each term has 4 factors because any spanning tree of the house graph has 4 edges.
def myHouseGraph():
H = DiGraph()
vertex_vars = ['a0','a1','a2','a3','a4']
edge_vars = ['a','b','c','d','e','f','A','B','C','D','E','F']
V=PolynomialRing(ZZ,names=vertex_vars,order='invlex')
V.inject_variables()
E=PolynomialRing(ZZ,names=edge_vars,order='lex')
E.inject_variables()
H.add_edge(a0,a1,E.gen(edge_vars.index('a')))
H.add_edge(a1,a0,E.gen(edge_vars.index('A')))
H.add_edge(a0,a2,E.gen(edge_vars.index('c')))
H.add_edge(a2,a0,E.gen(edge_vars.index('C')))
H.add_edge(a0,a4,E.gen(edge_vars.index('f')))
H.add_edge(a4,a0,E.gen(edge_vars.index('F')))
H.add_edge(a1,a2,E.gen(edge_vars.index('b')))
H.add_edge(a2,a1,E.gen(edge_vars.index('B')))
H.add_edge(a2,a3,E.gen(edge_vars.index('d')))
H.add_edge(a3,a2,E.gen(edge_vars.index('D')))
H.add_edge(a3,a4,E.gen(edge_vars.index('e')))
H.add_edge(a4,a3,E.gen(edge_vars.index('E')))
return H, E
H, P = myHouseGraph()
H.show(edge_labels=True,layout='spring')
Defining a0, a1, a2, a3, a4
Defining a, b, c, d, e, f, A, B, C, D, E, F
def tree_polynomial(G,P=None,root=None):
if P==None:
edge_vars = []
f = 0
#edge_coeff = []
for e in G.edges(sort=True):
f += e[2]
#print(e,e[2].monomials()[0],list(e[2].monomials()[0]))
print(f.variables())
edge_vars = list(f.variables())
# edge_coeff.append(e[2]/e[2].monomials()[0])
P = PolynomialRing(ZZ,names=edge_vars)
Trees = []
if G.is_directed():
if root != None:
T = P(0)
for b in G.in_branchings(source=root,spanning=True):
t = P(1)
for e in b.edges(sort=True):
t *= P(e[2])
T += t
return T
else:
for r in G.vertices(sort=True):
#print('r=',r)
T = P(0)
for b in G.in_branchings(source=r,spanning=True):
#b.show(figsize=10,edge_labels=true)
t = P(1)
for e in b.edges(sort=True):
#print(e)
#t *= P(e[2].monomials()[0])
t *= P(e[2])
T += t
Trees.append(T)
#return Trees
return Trees
else:
T = 0
for b in G.spanning_trees():
t = 1
for e in b.edges(sort=True):
t *= G.edge_label(e[0],e[1])
T += t
return P(T)
TT = tree_polynomial(H,P)
T = tree_polynomial(H,P,a0)
show(T)
\(\displaystyle b d e F + b e C F + b C D E + b C D F + d e A F + e A B F + e A C F + A B D E + A B D F + A C D E + A C D F\)
for v in H.vertices(sort=True):
T = tree_polynomial(H,P,v)
print(v,':',T)
a0 : b*d*e*F + b*e*C*F + b*C*D*E + b*C*D*F + d*e*A*F + e*A*B*F + e*A*C*F + A*B*D*E + A*B*D*F + A*C*D*E + A*C*D*F
a1 : a*d*e*F + a*e*B*F + a*e*C*F + a*B*D*E + a*B*D*F + a*C*D*E + a*C*D*F + c*e*B*F + c*B*D*E + c*B*D*F + f*B*D*E
a2 : a*b*e*F + a*b*D*E + a*b*D*F + b*c*e*F + b*c*D*E + b*c*D*F + b*f*D*E + c*e*A*F + c*A*D*E + c*A*D*F + f*A*D*E
a3 : a*b*d*E + a*b*d*F + b*c*d*E + b*c*d*F + b*d*f*E + b*f*C*E + c*d*A*E + c*d*A*F + d*f*A*E + f*A*B*E + f*A*C*E
a4 : a*b*d*e + b*c*d*e + b*d*e*f + b*e*f*C + b*f*C*D + c*d*e*A + d*e*f*A + e*f*A*B + e*f*A*C + f*A*B*D + f*A*C*D
def cycle_fluxes(G,P,verbose=True):
Cvert3 = []
Cvert2 = []
for C in G.all_simple_cycles():
if len(C) > 3:
Cvert3.append(C)
else:
Cvert2.append(C)
Cvert2.sort()
Cvert3Sort = Cvert3.copy()
Cvert3Sort.sort()
Cvert3SortP = Cvert3Sort.copy()
for i in range(len(Cvert3Sort)/2):
Cvert3SortP.remove(Cvert3SortP[i][::-1])
Cvert3SortP2 = []
for i in range(3,100): # arbitrary max cycle length
Cvert3SortP2 += [ C for C in Cvert3SortP if len(C)==i ]
Cvert3SortP = Cvert3SortP2.copy()
Cvert3SortM = []
for i in range(len(Cvert3SortP)):
Cvert3SortM.append(Cvert3SortP[i][::-1])
CvertP = Cvert3SortP + Cvert2
CvertM = Cvert3SortM + Cvert2
ncycles = len(CvertP)
CedgeP = []
CedgeM = []
ThetaP = []
ThetaM = []
for i in range(ncycles):
CPe = [ G.edge_label(CvertP[i][j],CvertP[i][j+1]) for j in range(len(CvertP[i])-1)]
CMe = [ G.edge_label(CvertM[i][j],CvertM[i][j+1]) for j in range(len(CvertM[i])-1)]
TPe = prod(CPe)
TMe = prod(CMe)
CedgeP.append(CPe)
ThetaP.append(TPe)
CedgeM.append(CMe)
ThetaM.append(TMe)
Tkappa = []
ThetaPMTkappa = []
for i in range(ncycles):
CP=CvertP[i]
#print('CP=',CP)
H = G.copy()
#show(H,edge_labels=True) ######
H.allow_multiple_edges(True)
#print([None]+CP)
H.merge_vertices([None]+CP) # None becomes vertex 0 (as opposed to a0, a1, ...)
if H.order():
#show(H,edge_labels=True) ######
T = tree_polynomial(H,P,root=0)
else:
T = P(1)
Tkappa.append(T)
ThetaPMTkappa.append((ThetaP[i]-ThetaM[i])*Tkappa[i])
if verbose:
print('\ntotal of',len(Cvert2),'+ 2 *',len(Cvert3)/2,'cycles.\n')
for i in range(ncycles):
H = G.copy()
for e in H.edges(sort=True):
H.set_edge_label(e[0],e[1],0)
for j in range(len(CvertP[i])-1):
H.set_edge_label(CvertP[i][j],CvertP[i][j+1],1)
H.set_edge_label(CvertM[i][j],CvertM[i][j+1],-1)
if ThetaP[i] == ThetaM[i]:
print(i,'=>',CvertP[i],'-',CvertM[i],'=>',CedgeP[i],'-',CedgeM[i],'=>',ThetaP[i]-ThetaM[i],'@ (',Tkappa[i],')')
else:
print(i,'=>',CvertP[i],'-',CvertM[i],'=>',CedgeP[i],'-',CedgeM[i],'=>',(ThetaP[i]-ThetaM[i]).factor(),'@ (',Tkappa[i],')')
H.show(color_by_label=True,edge_labels=True)
return ncycles, CedgeP, CedgeM, ThetaP, ThetaM, Tkappa, ThetaPMTkappa
Trees = tree_polynomial(H)
T = [0]*len(Trees)
for i in range(len(Trees)):
T[i] = Trees[i]
print('T[%s] = ' %(i),T[i])
TT = sum(T)
(a, b, c, d, e, f, A, B, C, D, E, F)
T[0] = A*B*D*E + b*C*D*E + A*C*D*E + b*d*e*F + d*e*A*F + e*A*B*F + b*e*C*F + e*A*C*F + A*B*D*F + b*C*D*F + A*C*D*F
T[1] = a*B*D*E + c*B*D*E + f*B*D*E + a*C*D*E + a*d*e*F + a*e*B*F + c*e*B*F + a*e*C*F + a*B*D*F + c*B*D*F + a*C*D*F
T[2] = a*b*D*E + b*c*D*E + b*f*D*E + c*A*D*E + f*A*D*E + a*b*e*F + b*c*e*F + c*e*A*F + a*b*D*F + b*c*D*F + c*A*D*F
T[3] = a*b*d*E + b*c*d*E + b*d*f*E + c*d*A*E + d*f*A*E + f*A*B*E + b*f*C*E + f*A*C*E + a*b*d*F + b*c*d*F + c*d*A*F
T[4] = a*b*d*e + b*c*d*e + b*d*e*f + c*d*e*A + d*e*f*A + e*f*A*B + b*e*f*C + e*f*A*C + f*A*B*D + b*f*C*D + f*A*C*D
ncycles, CedgeP, CedgeM, ThetaP, ThetaM, Tkappa, ThetaPMTkappa = cycle_fluxes(H,P,verbose=True)
for i in range(ncycles):
print('(',ThetaP[i],'-',ThetaM[i],') @ (',Tkappa[i],')')
total of 6 + 2 * 3 cycles.
0 => [a0, a1, a2, a0] - [a0, a2, a1, a0] => [a, b, C] - [c, B, A] => (-1) * (-a*b*C + c*A*B) @ ( e*F + D*E + D*F )
1 => [a0, a2, a3, a4, a0] - [a0, a4, a3, a2, a0] => [c, d, e, F] - [f, E, D, C] => (-1) * (-c*d*e*F + f*C*D*E) @ ( b + A )
2 => [a0, a1, a2, a3, a4, a0] - [a0, a4, a3, a2, a1, a0] => [a, b, d, e, F] - [f, E, D, B, A] => (-1) * (-a*b*d*e*F + f*A*B*D*E) @ ( 1 )
3 => [a0, a1, a0] - [a0, a1, a0] => [a, A] - [a, A] => 0 @ ( d*e*F + e*B*F + e*C*F + B*D*E + B*D*F + C*D*E + C*D*F )
4 => [a0, a2, a0] - [a0, a2, a0] => [c, C] - [c, C] => 0 @ ( b*e*F + b*D*E + b*D*F + e*A*F + A*D*E + A*D*F )
5 => [a0, a4, a0] - [a0, a4, a0] => [f, F] - [f, F] => 0 @ ( b*d*e + b*e*C + b*C*D + d*e*A + e*A*B + e*A*C + A*B*D + A*C*D )
6 => [a1, a2, a1] - [a1, a2, a1] => [b, B] - [b, B] => 0 @ ( a*e*F + a*D*E + a*D*F + c*e*F + c*D*E + c*D*F + f*D*E )
7 => [a2, a3, a2] - [a2, a3, a2] => [d, D] - [d, D] => 0 @ ( a*b*E + a*b*F + b*c*E + b*c*F + b*f*E + c*A*E + c*A*F + f*A*E )
8 => [a3, a4, a3] - [a3, a4, a3] => [e, E] - [e, E] => 0 @ ( a*b*d + b*c*d + b*d*f + b*f*C + c*d*A + d*f*A + f*A*B + f*A*C )
( a*b*C - c*A*B ) @ ( e*F + D*E + D*F )
( c*d*e*F - f*C*D*E ) @ ( b + A )
( a*b*d*e*F - f*A*B*D*E ) @ ( 1 )
( a*A - a*A ) @ ( d*e*F + e*B*F + e*C*F + B*D*E + B*D*F + C*D*E + C*D*F )
( c*C - c*C ) @ ( b*e*F + b*D*E + b*D*F + e*A*F + A*D*E + A*D*F )
( f*F - f*F ) @ ( b*d*e + b*e*C + b*C*D + d*e*A + e*A*B + e*A*C + A*B*D + A*C*D )
( b*B - b*B ) @ ( a*e*F + a*D*E + a*D*F + c*e*F + c*D*E + c*D*F + f*D*E )
( d*D - d*D ) @ ( a*b*E + a*b*F + b*c*E + b*c*F + b*f*E + c*A*E + c*A*F + f*A*E )
( e*E - e*E ) @ ( a*b*d + b*c*d + b*d*f + b*f*C + c*d*A + d*f*A + f*A*B + f*A*C )