/
stoichBondGraph.py
159 lines (126 loc) · 5.1 KB
/
stoichBondGraph.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
"""Convert a stoichiometric representation to BondGraphTools"""
import BondGraphTools as bgt
import numpy as np
import datetime
def model(s,filename=None):
""" Converts a stoichiometric matrix to the BondGraphTools bond graph representation
"""
## Extract info
N = s["N"] # Stoichiometric matrix
## Copy or create forward and reverse stoichiometric matrices
if "Nf" in s:
#print("Using Nf")
Nf = s["Nf"]
else:
Nf = -N*(N<0)
if "Nr" in s:
Nr = s["Nr"]
else:
Nr = N*(N>0)
shape = N.shape
n_X = shape[0] # Number of species (states)
n_V = shape[1] # Number of reactions (flows)
species = s["species"] # List of species names
reaction = s["reaction"] # List of reaction names
if "name" in s.keys():
name = s["name"]
else:
name = "none"
## Sanity check
if not ((n_X == len(species))):
print("Number of species (",len(species),") not equal to number of rows (", n_X,").")
if not ((n_V == len(reaction))):
print("Number of reactions (",len(reaction),") not equal to number of columns (", n_V,").")
## Create model header
indent = " "
delim = '"""'
headStr = ("import BondGraphTools as bgt\n"
"import sympy as sp\n\n"
"def model():\n"
"{0}{3} Acausal bond graph {1}.py\n"
"{0}Created by stoichBondGraph at {2}\n\n"
"{0}Usage:\n"
"{0}import {1}; model = {1}.model()\n"
"{0}{3}\n\n"
'{0}model = bgt.new(name="{1}")\n'
)
head = headStr.format(indent,name, datetime.datetime.now().ctime(),delim)
## print(head)
## Formatting
CeStr = ("\n{0}### Ce:{1}\n"
"\n{0}## Component Ce:{2}{1}\n"
"{0}K_{1} = sp.symbols('K_{1}')\n"
"{0}RT = sp.symbols('RT')\n"
"{0}{2}{1} = bgt.new('Ce',name='{1}',value={{'k':K_{1},'R':RT,'T':1}},library='BioChem')\n"
"{0}model.add({2}{1})\n"
)
junStr = ("\n{0}## Junction {1}:{2}\n"
"{0}{2} = bgt.new('{1}')\n"
"{0}model.add({2})\n"
)
bondStr = ("\n{0}## Bond from {1} to {2}\n"
"{0}bgt.connect({1},{2})\n"
)
conCR = ("\n{0}### Ce:{1} to Re:{2}\n")
conRC = ("\n{0}### Re:{1} to Ce:{2}\n")
prepend = "comp_"
junName = "{1}{0}_jun" # Junction name
junF = "{1}{0}_junF" # Forward junction name
junR = "{1}{0}_junR" # Reverse junction name
specStr = indent+"### Species\n"
## Create species
for spec in species:
specStr += CeStr.format(indent,spec,prepend)
specStr += junStr.format(indent,"0",junName.format(spec,prepend))
specStr += bondStr.format(indent,junName.format(spec,prepend),prepend+spec)
## print(specStr)
## Create reactions
ReStr = ("\n{0}### Re:{1}\n"
"\n{0}## Component Re:{2}{1}\n"
"{0}kappa_{1} = sp.symbols('kappa_{1}')\n"
"{0}RT = sp.symbols('RT')\n"
"{0}{2}{1} = bgt.new('Re',name='{1}',value={{'r':kappa_{1},'R':RT,'T':1}},library='BioChem')\n"
"{0}model.add({2}{1})\n"
)
reacStr = indent+"### Reactions\n"
for reac in reaction:
reacStr += ReStr.format(indent,reac,prepend)
reacStr += junStr.format(indent,"1",junF.format(reac,prepend))
portF = "("+prepend+reac+",0)"
reacStr += bondStr.format(indent,junF.format(reac,prepend),portF)
reacStr += junStr.format(indent,"1",junR.format(reac,prepend))
portR = "("+prepend+reac+",1)"
reacStr += bondStr.format(indent,portR,junR.format(reac,prepend))
## print(reacStr)
## Create the connections
connectStr = indent+"### Connections\n"
for i,row in enumerate(Nf):
for j,coeff in enumerate(row):
if (coeff>0):
for k in range(coeff):
#print(species[i],reaction[j])
connectStr += conCR.format(indent,species[i],reaction[j])
connectStr += bondStr.format(indent,junName.format(species[i],prepend),junF.format(reaction[j],prepend))
for i,row in enumerate(Nr):
for j,coeff in enumerate(row):
if (coeff>0):
for k in range(coeff):
connectStr += conRC.format(indent,reaction[j],species[i])
connectStr += bondStr.format(indent,junR.format(reaction[j],prepend),junName.format(species[i],prepend))
# if (coeff>0):
# for k in range(abs(coeff)):
# connectStr += conRC.format(indent,reaction[j],species[i])
# connectStr += bondStr.format(indent,junR.format(reaction[j],prepend),junName.format(species[i],prepend))
## print(connectStr)
## Print to file
if filename is None:
Filename = name+".py"
else:
Filename = filename+".py"
f = open(Filename,'w')
f.write(head)
f.write(specStr)
f.write(reacStr)
f.write(connectStr)
f.write("\n"+indent+"return model\n")
f.close