Error with exporting a myokit model to python 导出myokit模型到Python时出现错误

huangapple go评论94阅读模式
英文:

Error with exporting a myokit model to python

问题

I am trying to export a myokit model to python file using a myokit.formats.python.PythonExporter from here: https://myokit.readthedocs.io/en/stable/api_formats/python.html?highlight=export python# The code is like this:

  1. model1, protocol, x = myokit.load('dn-1985-if-gna.mmt')
  2. sim = myokit.Simulation(model1, protocol)
  3. exporter= myokit.formats.python.PythonExporter()
  4. exporter.model('dn-1985-if-gna_py.mmt', model1)

And I get the NotImplemented Error:

  1. NotImplementedError:

I suppose that I am possibly not creating the correct class object, but I have also tried this:

  1. exporter = myokit.formats.python.PythonExporter
  2. exporter.model('dn-1985-if-gna_py.mmt', model1)

The error is:

  1. TypeError: model() missing 1 required positional argument: 'model'

P.S. I have tried a couple of models, for example, this .mmt file https://dropmefiles.com/e3xjm

英文:

I am trying to export a myokit model to python file using a myokit.formats.python.PythonExporter from here:
https://myokit.readthedocs.io/en/stable/api_formats/python.html?highlight=export python#

The code is like this:

  1. model1, protocol, x = myokit.load('dn-1985-if-gna.mmt')
  2. sim = myokit.Simulation(model1, protocol)
  3. exporter= myokit.formats.python.PythonExporter()
  4. exporter.model('dn-1985-if-gna_py.mmt, model1)

And i get the NotImplemented Error:

  1. NotImplementedError Traceback (most recent call last)
  2. c:\More_Program_Files\Coding\Jupyter_files23_cellml_myokit23_07_04_myokit_tut.ipynb Cell 48 in ()
  3. 1 exporter= myokit.formats.python.PythonExporter()
  4. ----> 2 a = exporter.model('dn-1985-if-gna_copy_py.mmt', sim)
  5. File c:\More_Program_Files\Anaconda3\lib\site-packages\myokit\formats\__init__.py:61, in Exporter.model(self, path, model)
  6. 54 def model(self, path, model):
  7. 55 """
  8. 56 Exports a :class:`myokit.Model`.
  9. 57
  10. 58 The output will be stored in the **file** ``path``. A
  11. 59 :class:`myokit.ExportError` will be raised if any errors occur.
  12. 60 """
  13. ---> 61 raise NotImplementedError
  14. NotImplementedError:

I suppose that i am possibly not creating the correct class object, but I have tried also this:

  1. exporter= myokit.formats.python.PythonExporter
  2. exporter.model('dn-1985-if-gna_py.mmt', model1)

The error is:

  1. TypeError Traceback (most recent call last)
  2. c:\More_Program_Files\Coding\Jupyter_files23_cellml_myokit23_07_04_myokit_tut.ipynb Cell 48 in ()
  3. 1 exporter= myokit.formats.python.PythonExporter
  4. ----> 2 exporter.model('dn-1985-if-gna_copy_py.mmt', model)
  5. TypeError: model() missing 1 required positional argument: 'model'

P.S. I have tried a couple of models, for exapmle, this .mmt file https://dropmefiles.com/e3xjm

答案1

得分: 1

这是您提供的Python代码的翻译:

  1. #!/usr/bin/env python
  2. #
  3. # 生成日期:2023年5月13日02:47:17
  4. #
  5. import math
  6. #
  7. # 组件和变量
  8. #
  9. class CEngine(object):
  10. def __init__(self):
  11. self.time = None
  12. self.pace = None
  13. self._constants()
  14. self.init()
  15. def _constants(self):
  16. """
  17. 设置常量值
  18. """
  19. def init(self):
  20. """
  21. 将状态变量重置为其初始值
  22. """
  23. def update(self):
  24. """
  25. 重新计算当前时间和状态的所有值
  26. """
  27. c_engine.pace = engine.pace
  28. c_engine.time = engine.time
  29. class CIk1(object):
  30. def __init__(self):
  31. self.IK1 = None
  32. self._constants()
  33. self.init()
  34. def _constants(self):
  35. """
  36. 设置常量值
  37. """
  38. def init(self):
  39. """
  40. 将状态变量重置为其初始值
  41. """
  42. def update(self):
  43. """
  44. 重新计算当前时间和状态的所有值
  45. """
  46. self.IK1 = 0.35 * (4.0 * (math.exp(0.04 * (c_membrane.V + 85.0)) - 1.0) / (math.exp(0.08 * (c_membrane.V + 53.0)) + math.exp(0.04 * (c_membrane.V + 53.0))) + 0.2 * (c_membrane.V + 23.0) / (1.0 - math.exp((-0.04) * (c_membrane.V + 23.0))))
  47. class CIna(object):
  48. def __init__(self):
  49. self.gNaBar = None
  50. self.gNaC = None
  51. self.ENa = None
  52. self.INa = None
  53. self.m = None
  54. self.d_m = None
  55. self.ina_m_alpha = None
  56. self.ina_m_beta = None
  57. self.h = None
  58. self.d_h = None
  59. self.ina_h_alpha = None
  60. self.ina_h_beta = None
  61. self.j = None
  62. self.d_j = None
  63. self.ina_j_alpha = None
  64. self.ina_j_beta = None
  65. self._constants()
  66. self.init()
  67. def _constants(self):
  68. """
  69. 设置常量值
  70. """
  71. self.ENa = 50.0
  72. self.gNaBar = 4.0
  73. self.gNaC = 0.003
  74. def init(self):
  75. """
  76. 将状态变量重置为其初始值
  77. """
  78. self.m = 0.01
  79. self.h = 0.99
  80. self.j = 0.98
  81. def update(self):
  82. """
  83. 重新计算当前时间和状态的所有值
  84. """
  85. self.ina_h_alpha = 0.126 * math.exp((-0.25) * (c_membrane.V + 77.0))
  86. self.ina_h_beta = 1.7 / (1.0 + math.exp((-0.082) * (c_membrane.V + 22.5)))
  87. self.d_h = self.ina_h_alpha * (1.0 - self.h) - self.ina_h_beta * self.h
  88. self.ina_j_alpha = 0.055 * math.exp((-0.25) * (c_membrane.V + 78.0)) / (1.0 + math.exp((-0.2) * (c_membrane.V + 78.0)))
  89. self.ina_j_beta = 0.3 / (1.0 + math.exp((-0.1) * (c_membrane.V + 32.0)))
  90. self.d_j = self.ina_j_alpha * (1.0 - self.j) - self.ina_j_beta * self.j
  91. self.ina_m_alpha = (c_membrane.V + 47.0) / (1.0 - math.exp((-0.1) * (c_membrane.V + 47.0)))
  92. self.ina_m_beta = 40.0 * math.exp((-0.056) * (c_membrane.V + 72.0))
  93. self.d_m = self.ina_m_alpha * (1.0 - self.m) - self.ina_m_beta * self.m
  94. self.INa = (self.gNaBar * self.m ** 3.0 * self.h * self.j + self.gNaC) * (c_membrane.V - self.ENa)
  95. class CIsi(object):
  96. def __init__(self):
  97. self.gsBar = None
  98. self.Es = None
  99. self.Isi = None
  100. self.d = None
  101. self.d_d = None
  102. self.isi_d_alpha = None
  103. self.isi_d_beta = None
  104. self.f = None
  105. self.d_f = None
  106. self.isi_f_alpha = None
  107. self.isi_f_beta = None
  108. self.Cai = None
  109. self.d_cai = None
  110. self._constants()
  111. self.init()
  112. def _constants(self):
  113. """
  114. 设置常量值
  115. """
  116. self.gsBar = 0.09
  117. def init(self):
  118. """
  119. 将状态变量重置为其初始值
  120. """
  121. self.d = 0.003
  122. self.f = 0.99
  123. self.Cai = 2e-07
  124. def update(self):
  125. """
  126. 重新计算当前时间和状态的所有值
  127. """
  128. self.Es = (-82.3) - 13.0287 * math.log(self.Cai)
  129. self.isi_d_alpha = 0.095 * math.exp((-0.01) * (c_membrane.V + (-5.0))) / (math.exp((-0.072) * (c_membrane.V + (-5.0))) + 1.0)
  130. self.isi_d_beta = 0.07 * math.exp((-0.017) * (c_membrane.V + 44.0)) / (math.exp(0.05 * (c_membrane.V + 44.0)) + 1.0)
  131. self.d_d = self.isi_d_alpha * (1.0 - self.d) - self.isi_d_beta * self.d
  132. self.isi_f_alpha = 0.012 * math.exp((-0.008) * (c_membrane.V + 28.0)) / (math.exp(0.15 * (c_membrane.V + 28.0)) + 1.0)
  133. self.isi_f_beta = 0.0065 * math.exp((-0.02) * (c_membrane.V + 30.0)) / (math.exp((-0.2) * (c_membrane.V + 30.0)) + 1.0)
  134. self.d_f = self.isi_f_alpha * (1.0 - self.f) - self.isi_f_beta
  135. <details>
  136. <summary>英文:</summary>
  137. I have never used this library before so there could be a lot of mistakes in my comment but hopefully it might show you the correct path.
  138. I downloaded a sample `br-1977.mmt` model from the documentation and I used this following script to export the python code.
  139. ```python
  140. import myokit
  141. model = myokit.load(&quot;./br-1977.mmt&quot;)
  142. print(model)
  143. # (&lt;Model(Beeler-Reuter-1977)&gt;, &lt;myokit._protocol.Protocol object at 0x7fd77bb52bd0&gt;, &quot;import matplotlib.pyplot as plt\nimport myokit\n\n#\n# This example file plots a single AP using the model develope by Beeler\n# and Reuter.\n#\n\n# Get model and protocol, create simulation\nm = get_model()\np = get_protocol()\ns = myokit.Simulation(m, p)\n\n# Run simulation\nd = s.run(1000)\n\n# Display the result\nplt.figure()\nplt.plot(d[&#39;engine.time&#39;], d[&#39;membrane.V&#39;])\nplt.show()\n\n&quot;)
  144. exporter = myokit.formats.python.PythonExporter()
  145. exporter.runnable(model=model[0], path=&quot;./br-1977&quot;)

Here the loaded model is a tuple, so I passed the first item from the tuple as it seemed to be the model object required by the function.

And I got this as the generated code.

  1. #!/usr/bin/env python
  2. #
  3. # Generated on 2023-05-13 02:47:17
  4. #
  5. import math
  6. #
  7. # Components and variables
  8. #
  9. class CEngine(object):
  10. def __init__(self):
  11. self.time = None
  12. self.pace = None
  13. self._constants()
  14. self.init()
  15. def _constants(self):
  16. &quot;&quot;&quot;
  17. Sets the constant values
  18. &quot;&quot;&quot;
  19. def init(self):
  20. &quot;&quot;&quot;
  21. Resets the state variables to their initial values
  22. &quot;&quot;&quot;
  23. def update(self):
  24. &quot;&quot;&quot;
  25. Re-calculates all values for the current time and state
  26. &quot;&quot;&quot;
  27. c_engine.pace = engine.pace
  28. c_engine.time = engine.time
  29. class CIk1(object):
  30. def __init__(self):
  31. self.IK1 = None
  32. self._constants()
  33. self.init()
  34. def _constants(self):
  35. &quot;&quot;&quot;
  36. Sets the constant values
  37. &quot;&quot;&quot;
  38. def init(self):
  39. &quot;&quot;&quot;
  40. Resets the state variables to their initial values
  41. &quot;&quot;&quot;
  42. def update(self):
  43. &quot;&quot;&quot;
  44. Re-calculates all values for the current time and state
  45. &quot;&quot;&quot;
  46. self.IK1 = 0.35 * (4.0 * (math.exp(0.04 * (c_membrane.V + 85.0)) - 1.0) / (math.exp(0.08 * (c_membrane.V + 53.0)) + math.exp(0.04 * (c_membrane.V + 53.0))) + 0.2 * (c_membrane.V + 23.0) / (1.0 - math.exp((-0.04) * (c_membrane.V + 23.0))))
  47. class CIna(object):
  48. def __init__(self):
  49. self.gNaBar = None
  50. self.gNaC = None
  51. self.ENa = None
  52. self.INa = None
  53. self.m = None
  54. self.d_m = None
  55. self.ina_m_alpha = None
  56. self.ina_m_beta = None
  57. self.h = None
  58. self.d_h = None
  59. self.ina_h_alpha = None
  60. self.ina_h_beta = None
  61. self.j = None
  62. self.d_j = None
  63. self.ina_j_alpha = None
  64. self.ina_j_beta = None
  65. self._constants()
  66. self.init()
  67. def _constants(self):
  68. &quot;&quot;&quot;
  69. Sets the constant values
  70. &quot;&quot;&quot;
  71. self.ENa = 50.0
  72. self.gNaBar = 4.0
  73. self.gNaC = 0.003
  74. def init(self):
  75. &quot;&quot;&quot;
  76. Resets the state variables to their initial values
  77. &quot;&quot;&quot;
  78. self.m = 0.01
  79. self.h = 0.99
  80. self.j = 0.98
  81. def update(self):
  82. &quot;&quot;&quot;
  83. Re-calculates all values for the current time and state
  84. &quot;&quot;&quot;
  85. self.ina_h_alpha = 0.126 * math.exp((-0.25) * (c_membrane.V + 77.0))
  86. self.ina_h_beta = 1.7 / (1.0 + math.exp((-0.082) * (c_membrane.V + 22.5)))
  87. self.d_h = self.ina_h_alpha * (1.0 - self.h) - self.ina_h_beta * self.h
  88. self.ina_j_alpha = 0.055 * math.exp((-0.25) * (c_membrane.V + 78.0)) / (1.0 + math.exp((-0.2) * (c_membrane.V + 78.0)))
  89. self.ina_j_beta = 0.3 / (1.0 + math.exp((-0.1) * (c_membrane.V + 32.0)))
  90. self.d_j = self.ina_j_alpha * (1.0 - self.j) - self.ina_j_beta * self.j
  91. self.ina_m_alpha = (c_membrane.V + 47.0) / (1.0 - math.exp((-0.1) * (c_membrane.V + 47.0)))
  92. self.ina_m_beta = 40.0 * math.exp((-0.056) * (c_membrane.V + 72.0))
  93. self.d_m = self.ina_m_alpha * (1.0 - self.m) - self.ina_m_beta * self.m
  94. self.INa = (self.gNaBar * self.m ** 3.0 * self.h * self.j + self.gNaC) * (c_membrane.V - self.ENa)
  95. class CIsi(object):
  96. def __init__(self):
  97. self.gsBar = None
  98. self.Es = None
  99. self.Isi = None
  100. self.d = None
  101. self.d_d = None
  102. self.isi_d_alpha = None
  103. self.isi_d_beta = None
  104. self.f = None
  105. self.d_f = None
  106. self.isi_f_alpha = None
  107. self.isi_f_beta = None
  108. self.Cai = None
  109. self.d_cai = None
  110. self._constants()
  111. self.init()
  112. def _constants(self):
  113. &quot;&quot;&quot;
  114. Sets the constant values
  115. &quot;&quot;&quot;
  116. self.gsBar = 0.09
  117. def init(self):
  118. &quot;&quot;&quot;
  119. Resets the state variables to their initial values
  120. &quot;&quot;&quot;
  121. self.d = 0.003
  122. self.f = 0.99
  123. self.Cai = 2e-07
  124. def update(self):
  125. &quot;&quot;&quot;
  126. Re-calculates all values for the current time and state
  127. &quot;&quot;&quot;
  128. self.Es = (-82.3) - 13.0287 * math.log(self.Cai)
  129. self.isi_d_alpha = 0.095 * math.exp((-0.01) * (c_membrane.V + (-5.0))) / (math.exp((-0.072) * (c_membrane.V + (-5.0))) + 1.0)
  130. self.isi_d_beta = 0.07 * math.exp((-0.017) * (c_membrane.V + 44.0)) / (math.exp(0.05 * (c_membrane.V + 44.0)) + 1.0)
  131. self.d_d = self.isi_d_alpha * (1.0 - self.d) - self.isi_d_beta * self.d
  132. self.isi_f_alpha = 0.012 * math.exp((-0.008) * (c_membrane.V + 28.0)) / (math.exp(0.15 * (c_membrane.V + 28.0)) + 1.0)
  133. self.isi_f_beta = 0.0065 * math.exp((-0.02) * (c_membrane.V + 30.0)) / (math.exp((-0.2) * (c_membrane.V + 30.0)) + 1.0)
  134. self.d_f = self.isi_f_alpha * (1.0 - self.f) - self.isi_f_beta * self.f
  135. self.Isi = self.gsBar * self.d * self.f * (c_membrane.V - self.Es)
  136. self.d_cai = (-1e-07) * self.Isi + 0.07 * (1e-07 - self.Cai)
  137. class CIx1(object):
  138. def __init__(self):
  139. self.Ix1 = None
  140. self.x1 = None
  141. self.d_x1 = None
  142. self.ix1_x1_alpha = None
  143. self.ix1_x1_beta = None
  144. self._constants()
  145. self.init()
  146. def _constants(self):
  147. &quot;&quot;&quot;
  148. Sets the constant values
  149. &quot;&quot;&quot;
  150. def init(self):
  151. &quot;&quot;&quot;
  152. Resets the state variables to their initial values
  153. &quot;&quot;&quot;
  154. self.x1 = 0.0004
  155. def update(self):
  156. &quot;&quot;&quot;
  157. Re-calculates all values for the current time and state
  158. &quot;&quot;&quot;
  159. self.Ix1 = self.x1 * 0.8 * (math.exp(0.04 * (c_membrane.V + 77.0)) - 1.0) / math.exp(0.04 * (c_membrane.V + 35.0))
  160. self.ix1_x1_alpha = 0.0005 * math.exp(0.083 * (c_membrane.V + 50.0)) / (math.exp(0.057 * (c_membrane.V + 50.0)) + 1.0)
  161. self.ix1_x1_beta = 0.0013 * math.exp((-0.06) * (c_membrane.V + 20.0)) / (math.exp((-0.04) * (c_membrane.V + 333.0)) + 1.0)
  162. self.d_x1 = self.ix1_x1_alpha * (1.0 - self.x1) - self.ix1_x1_beta * self.x1
  163. class CStimulus(object):
  164. def __init__(self):
  165. self.amplitude = None
  166. self.IStim = None
  167. self._constants()
  168. self.init()
  169. def _constants(self):
  170. &quot;&quot;&quot;
  171. Sets the constant values
  172. &quot;&quot;&quot;
  173. self.amplitude = 25.0
  174. def init(self):
  175. &quot;&quot;&quot;
  176. Resets the state variables to their initial values
  177. &quot;&quot;&quot;
  178. def update(self):
  179. &quot;&quot;&quot;
  180. Re-calculates all values for the current time and state
  181. &quot;&quot;&quot;
  182. self.IStim = c_engine.pace * self.amplitude
  183. class CMembrane(object):
  184. def __init__(self):
  185. self.C = None
  186. self.V = None
  187. self.d_v = None
  188. self._constants()
  189. self.init()
  190. def _constants(self):
  191. &quot;&quot;&quot;
  192. Sets the constant values
  193. &quot;&quot;&quot;
  194. self.C = 1.0
  195. def init(self):
  196. &quot;&quot;&quot;
  197. Resets the state variables to their initial values
  198. &quot;&quot;&quot;
  199. self.V = -84.622
  200. def update(self):
  201. &quot;&quot;&quot;
  202. Re-calculates all values for the current time and state
  203. &quot;&quot;&quot;
  204. self.d_v = (-(1.0 / self.C)) * (c_ik1.IK1 + c_ix1.Ix1 + c_ina.INa + c_isi.Isi - c_stimulus.IStim)
  205. #
  206. # Engine component
  207. #
  208. class Engine(object):
  209. &quot;&quot;&quot;
  210. Calculates the derivatives in the current state
  211. &quot;&quot;&quot;
  212. def __init__(self):
  213. self.pace = 0.0
  214. self.time = 0.0
  215. def update(self):
  216. c_engine.update()
  217. c_ik1.update()
  218. c_ina.update()
  219. c_isi.update()
  220. c_ix1.update()
  221. c_stimulus.update()
  222. c_membrane.update()
  223. # Remaining equations
  224. #
  225. # Create objects, set initial values
  226. #
  227. def init():
  228. &quot;&quot;&quot; (Re-)Initializes the model &quot;&quot;&quot;
  229. global c_engine
  230. global c_ik1
  231. global c_ina
  232. global c_isi
  233. global c_ix1
  234. global c_stimulus
  235. global c_membrane
  236. global engine
  237. c_engine = CEngine()
  238. c_ik1 = CIk1()
  239. c_ina = CIna()
  240. c_isi = CIsi()
  241. c_ix1 = CIx1()
  242. c_stimulus = CStimulus()
  243. c_membrane = CMembrane()
  244. engine = Engine()
  245. #
  246. # Update function (rhs function, takes a single step)
  247. #
  248. def update(stepSize):
  249. &quot;&quot;&quot; Calculates all derivatives, update state, advances time &quot;&quot;&quot;
  250. engine.update()
  251. c_membrane.V += stepSize * c_membrane.d_v
  252. c_ina.m += stepSize * c_ina.d_m
  253. c_ina.h += stepSize * c_ina.d_h
  254. c_ina.j += stepSize * c_ina.d_j
  255. c_isi.d += stepSize * c_isi.d_d
  256. c_isi.f += stepSize * c_isi.d_f
  257. c_ix1.x1 += stepSize * c_ix1.d_x1
  258. c_isi.Cai += stepSize * c_isi.d_cai
  259. #
  260. # State vector returning function
  261. #
  262. def state():
  263. &quot;&quot;&quot; Returns the state vector &quot;&quot;&quot;
  264. return [c_membrane.V,
  265. c_ina.m,
  266. c_ina.h,
  267. c_ina.j,
  268. c_isi.d,
  269. c_isi.f,
  270. c_ix1.x1,
  271. c_isi.Cai]
  272. #
  273. # State vector printing function
  274. #
  275. def print_state():
  276. &quot;&quot;&quot; Prints the current state to the screen &quot;&quot;&quot;
  277. f = &quot;{:&lt;10} {:&lt;20} {:&lt;20}&quot;
  278. print(&quot;-&quot;*80)
  279. print(f.format(&quot;Name&quot;, &quot;State value&quot;, &quot;Derivative&quot;))
  280. f = &quot;{: &lt;10} {:&lt; 20.13e} {:&lt; 20.13e}&quot;
  281. print(f.format(&quot;membrane.V&quot;, c_membrane.V, c_membrane.d_v))
  282. print(f.format(&quot;ina.m&quot;, c_ina.m, c_ina.d_m))
  283. print(f.format(&quot;ina.h&quot;, c_ina.h, c_ina.d_h))
  284. print(f.format(&quot;ina.j&quot;, c_ina.j, c_ina.d_j))
  285. print(f.format(&quot;isi.d&quot;, c_isi.d, c_isi.d_d))
  286. print(f.format(&quot;isi.f&quot;, c_isi.f, c_isi.d_f))
  287. print(f.format(&quot;ix1.x1&quot;, c_ix1.x1, c_ix1.d_x1))
  288. print(f.format(&quot;isi.Cai&quot;, c_isi.Cai, c_isi.d_cai))
  289. #
  290. # Test step function
  291. #
  292. def test_step():
  293. &quot;&quot;&quot; Calculates and prints the initial derivatives &quot;&quot;&quot;
  294. init()
  295. engine.update()
  296. print_state()
  297. #
  298. # Pacing
  299. #
  300. class Protocol(object):
  301. &quot;&quot;&quot; Holds an ordered set of ProtocolEvent objects &quot;&quot;&quot;
  302. def __init__(self):
  303. super(Protocol, self).__init__()
  304. self.head = None
  305. def add(self, e):
  306. &quot;&quot;&quot; Schedules an event &quot;&quot;&quot;
  307. if self.head is None:
  308. self.head = e
  309. return
  310. if e.start &lt; self.head.start:
  311. e.next = self.head
  312. self.head = e
  313. return
  314. f = self.head
  315. while (f.next is not None) and (e.start &gt; f.next.start):
  316. f = f.next
  317. e.next = f.next
  318. f.next = e
  319. def pop(self):
  320. &quot;&quot;&quot; Returns the next event &quot;&quot;&quot;
  321. e = self.head
  322. if self.head is not None:
  323. self.head = self.head.next
  324. return e
  325. class ProtocolEvent(object):
  326. def __init__(self, level, start, duration, period=0, multiplier=0):
  327. super(ProtocolEvent, self).__init__()
  328. self.level = float(level)
  329. self.start = float(start)
  330. self.duration = float(duration)
  331. if self.duration &lt;= 0:
  332. raise Exception(&#39;Duration must be greater than zero&#39;)
  333. self.period = float(period)
  334. if self.period &lt; 0:
  335. raise Exception(&#39;multiplier must be zero or greater&#39;)
  336. self.multiplier = int(multiplier)
  337. if self.multiplier &lt; 0:
  338. raise Exception(&#39;Multiplier must be zero or greater&#39;)
  339. if self.period == 0 and self.multiplier &gt; 0:
  340. raise Exception(&#39;Non-periodic event cannot occur more than once&#39;)
  341. self.next = None
  342. def pacing_protocol():
  343. pacing = Protocol()
  344. pacing.add(ProtocolEvent(1.0, 100.0, 0.5, 1000.0, 0))
  345. return pacing
  346. #
  347. # Solver
  348. #
  349. def beat(stepSmall = 0.005, stepLarge = 0.01):
  350. &quot;&quot;&quot;
  351. Simulates a single beat
  352. &quot;&quot;&quot;
  353. tmin = 0
  354. tmax = 1000
  355. # Feedback
  356. outInt = int(math.ceil(tmax / 10.0))
  357. outPos = engine.time
  358. outVal = 0
  359. # Logging
  360. logInt = 1
  361. logPos = engine.time
  362. log = []
  363. # Stepsize
  364. stepSize = stepSmall
  365. hadPulse = False
  366. vInit = c_membrane.V
  367. # Pacing
  368. pacing = pacing_protocol()
  369. next = pacing.head
  370. while next and next.start &lt; tmin:
  371. next = next.next
  372. if next.start &lt; tmin:
  373. next = None
  374. fire = None
  375. fireDown = 0
  376. stopTime = min(next.start, tmin + stepSize)
  377. print(&#39;Starting integration with step sizes &#39; + str(stepSmall) + &#39; and &#39; + str(stepLarge) + &#39;.&#39;)
  378. while engine.time &lt; tmax:
  379. update(stopTime - engine.time)
  380. engine.time = stopTime
  381. # Event over
  382. if (fire and engine.time &gt;= fireDown):
  383. engine.pace = 0
  384. fire = None
  385. # New event
  386. if (next and engine.time &gt;= next.start):
  387. fire = next
  388. next = next.next
  389. engine.pace = fire.level
  390. fireDown = fire.start + fire.duration
  391. if fire.period &gt; 0:
  392. if fire.multiplier == 1:
  393. fire.period = 0
  394. else:
  395. if fire.multiplier &gt; 1:
  396. fire.multiplier -=1
  397. fire.start += fire.period
  398. pacing.add(fire)
  399. next = pacing.head
  400. # User feedback
  401. if engine.time &gt;= outPos and outVal &lt; 100:
  402. print(str(outVal) + &quot;%&quot;)
  403. outVal += 10
  404. outPos += outInt
  405. # Logging
  406. if engine.time &gt;= logPos:
  407. log.append((engine.time, state()))
  408. logPos += logInt
  409. # Step size update
  410. if fire: # or c_membrane.V &gt; -70:
  411. if stepSize != stepSmall:
  412. print(&quot;Small steps&quot;)
  413. stepSize = stepSmall
  414. else:
  415. if stepSize != stepLarge:
  416. print(&quot;Big steps&quot;)
  417. stepSize = stepLarge
  418. # Set next time
  419. stopTime = engine.time + stepSize
  420. if fire and fireDown &lt; stopTime:
  421. stopTime = fireDown
  422. if next and next.start &lt; stopTime:
  423. stopTime = next.start
  424. if logPos &lt; stopTime:
  425. stopTime = logPos
  426. print(&quot;100% done&quot;)
  427. print(&quot;t = &quot; + str(engine.time))
  428. print_state()
  429. return log
  430. #
  431. # Run if loaded as main script
  432. #
  433. if __name__ == &#39;__main__&#39;:
  434. small = 0.005
  435. large = 0.01
  436. go = True
  437. done = False
  438. while go:
  439. go = False
  440. try:
  441. init()
  442. data = beat(small, large)
  443. done = True
  444. except ArithmeticError as e:
  445. print(&#39;Arithmetic error occurred&#39;)
  446. y = &#39;Continue with smaller stepsize? (y/n): &#39;
  447. try:
  448. y = raw_input(y) # Python 2
  449. except NameError:
  450. y = input(y) # Python 3
  451. if y.lower()[0:1] == &#39;y&#39;:
  452. small /= 2
  453. large /= 2
  454. go = True
  455. if done:
  456. print(&#39;Showing result...&#39;)
  457. x = []
  458. y = []
  459. for time, state in data:
  460. x.append(time)
  461. y.append(state[0])
  462. import matplotlib.pyplot as py
  463. plot = py.plot(x, y)
  464. py.show()

huangapple
  • 本文由 发表于 2023年5月13日 08:28:18
  • 转载请务必保留本文链接:https://go.coder-hub.com/76240593.html
匿名

发表评论

匿名网友

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

确定