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:
model1, protocol, x = myokit.load('dn-1985-if-gna.mmt')
sim = myokit.Simulation(model1, protocol)
exporter= myokit.formats.python.PythonExporter()
exporter.model('dn-1985-if-gna_py.mmt, model1)
And i get the NotImplemented Error:
NotImplementedError Traceback (most recent call last)
c:\More_Program_Files\Coding\Jupyter_files\2023_cellml_myokit\2023_07_04_myokit_tut.ipynb Cell 48 in ()
1 exporter= myokit.formats.python.PythonExporter()
----> 2 a = exporter.model('dn-1985-if-gna_copy_py.mmt', sim)
File c:\More_Program_Files\Anaconda3\lib\site-packages\myokit\formats\__init__.py:61, in Exporter.model(self, path, model)
54 def model(self, path, model):
55 """
56 Exports a :class:`myokit.Model`.
57
58 The output will be stored in the **file** ``path``. A
59 :class:`myokit.ExportError` will be raised if any errors occur.
60 """
---> 61 raise NotImplementedError
NotImplementedError:
I suppose that i am possibly not creating the correct class object, but I have tried also this:
exporter= myokit.formats.python.PythonExporter
exporter.model('dn-1985-if-gna_py.mmt', model1)
The error is:
TypeError Traceback (most recent call last)
c:\More_Program_Files\Coding\Jupyter_files\2023_cellml_myokit\2023_07_04_myokit_tut.ipynb Cell 48 in ()
1 exporter= myokit.formats.python.PythonExporter
----> 2 exporter.model('dn-1985-if-gna_copy_py.mmt', model)
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
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.
I downloaded a sample br-1977.mmt model from the documentation and I used this following script to export the python code.
import myokit
model = myokit.load("./br-1977.mmt")
print(model)
# (<Model(Beeler-Reuter-1977)>, <myokit._protocol.Protocol object at 0x7fd77bb52bd0>, "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['engine.time'], d['membrane.V'])\nplt.show()\n\n")
exporter = myokit.formats.python.PythonExporter()
exporter.runnable(model=model[0], path="./br-1977")
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.
#!/usr/bin/env python
#
# Generated on 2023-05-13 02:47:17
#
import math
#
# Components and variables
#
class CEngine(object):
def __init__(self):
self.time = None
self.pace = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
def init(self):
"""
Resets the state variables to their initial values
"""
def update(self):
"""
Re-calculates all values for the current time and state
"""
c_engine.pace = engine.pace
c_engine.time = engine.time
class CIk1(object):
def __init__(self):
self.IK1 = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
def init(self):
"""
Resets the state variables to their initial values
"""
def update(self):
"""
Re-calculates all values for the current time and state
"""
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))))
class CIna(object):
def __init__(self):
self.gNaBar = None
self.gNaC = None
self.ENa = None
self.INa = None
self.m = None
self.d_m = None
self.ina_m_alpha = None
self.ina_m_beta = None
self.h = None
self.d_h = None
self.ina_h_alpha = None
self.ina_h_beta = None
self.j = None
self.d_j = None
self.ina_j_alpha = None
self.ina_j_beta = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
self.ENa = 50.0
self.gNaBar = 4.0
self.gNaC = 0.003
def init(self):
"""
Resets the state variables to their initial values
"""
self.m = 0.01
self.h = 0.99
self.j = 0.98
def update(self):
"""
Re-calculates all values for the current time and state
"""
self.ina_h_alpha = 0.126 * math.exp((-0.25) * (c_membrane.V + 77.0))
self.ina_h_beta = 1.7 / (1.0 + math.exp((-0.082) * (c_membrane.V + 22.5)))
self.d_h = self.ina_h_alpha * (1.0 - self.h) - self.ina_h_beta * self.h
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)))
self.ina_j_beta = 0.3 / (1.0 + math.exp((-0.1) * (c_membrane.V + 32.0)))
self.d_j = self.ina_j_alpha * (1.0 - self.j) - self.ina_j_beta * self.j
self.ina_m_alpha = (c_membrane.V + 47.0) / (1.0 - math.exp((-0.1) * (c_membrane.V + 47.0)))
self.ina_m_beta = 40.0 * math.exp((-0.056) * (c_membrane.V + 72.0))
self.d_m = self.ina_m_alpha * (1.0 - self.m) - self.ina_m_beta * self.m
self.INa = (self.gNaBar * self.m ** 3.0 * self.h * self.j + self.gNaC) * (c_membrane.V - self.ENa)
class CIsi(object):
def __init__(self):
self.gsBar = None
self.Es = None
self.Isi = None
self.d = None
self.d_d = None
self.isi_d_alpha = None
self.isi_d_beta = None
self.f = None
self.d_f = None
self.isi_f_alpha = None
self.isi_f_beta = None
self.Cai = None
self.d_cai = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
self.gsBar = 0.09
def init(self):
"""
Resets the state variables to their initial values
"""
self.d = 0.003
self.f = 0.99
self.Cai = 2e-07
def update(self):
"""
Re-calculates all values for the current time and state
"""
self.Es = (-82.3) - 13.0287 * math.log(self.Cai)
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)
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)
self.d_d = self.isi_d_alpha * (1.0 - self.d) - self.isi_d_beta * self.d
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)
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)
self.d_f = self.isi_f_alpha * (1.0 - self.f) - self.isi_f_beta * self.f
self.Isi = self.gsBar * self.d * self.f * (c_membrane.V - self.Es)
self.d_cai = (-1e-07) * self.Isi + 0.07 * (1e-07 - self.Cai)
class CIx1(object):
def __init__(self):
self.Ix1 = None
self.x1 = None
self.d_x1 = None
self.ix1_x1_alpha = None
self.ix1_x1_beta = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
def init(self):
"""
Resets the state variables to their initial values
"""
self.x1 = 0.0004
def update(self):
"""
Re-calculates all values for the current time and state
"""
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))
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)
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)
self.d_x1 = self.ix1_x1_alpha * (1.0 - self.x1) - self.ix1_x1_beta * self.x1
class CStimulus(object):
def __init__(self):
self.amplitude = None
self.IStim = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
self.amplitude = 25.0
def init(self):
"""
Resets the state variables to their initial values
"""
def update(self):
"""
Re-calculates all values for the current time and state
"""
self.IStim = c_engine.pace * self.amplitude
class CMembrane(object):
def __init__(self):
self.C = None
self.V = None
self.d_v = None
self._constants()
self.init()
def _constants(self):
"""
Sets the constant values
"""
self.C = 1.0
def init(self):
"""
Resets the state variables to their initial values
"""
self.V = -84.622
def update(self):
"""
Re-calculates all values for the current time and state
"""
self.d_v = (-(1.0 / self.C)) * (c_ik1.IK1 + c_ix1.Ix1 + c_ina.INa + c_isi.Isi - c_stimulus.IStim)
#
# Engine component
#
class Engine(object):
"""
Calculates the derivatives in the current state
"""
def __init__(self):
self.pace = 0.0
self.time = 0.0
def update(self):
c_engine.update()
c_ik1.update()
c_ina.update()
c_isi.update()
c_ix1.update()
c_stimulus.update()
c_membrane.update()
# Remaining equations
#
# Create objects, set initial values
#
def init():
""" (Re-)Initializes the model """
global c_engine
global c_ik1
global c_ina
global c_isi
global c_ix1
global c_stimulus
global c_membrane
global engine
c_engine = CEngine()
c_ik1 = CIk1()
c_ina = CIna()
c_isi = CIsi()
c_ix1 = CIx1()
c_stimulus = CStimulus()
c_membrane = CMembrane()
engine = Engine()
#
# Update function (rhs function, takes a single step)
#
def update(stepSize):
""" Calculates all derivatives, update state, advances time """
engine.update()
c_membrane.V += stepSize * c_membrane.d_v
c_ina.m += stepSize * c_ina.d_m
c_ina.h += stepSize * c_ina.d_h
c_ina.j += stepSize * c_ina.d_j
c_isi.d += stepSize * c_isi.d_d
c_isi.f += stepSize * c_isi.d_f
c_ix1.x1 += stepSize * c_ix1.d_x1
c_isi.Cai += stepSize * c_isi.d_cai
#
# State vector returning function
#
def state():
""" Returns the state vector """
return [c_membrane.V,
c_ina.m,
c_ina.h,
c_ina.j,
c_isi.d,
c_isi.f,
c_ix1.x1,
c_isi.Cai]
#
# State vector printing function
#
def print_state():
""" Prints the current state to the screen """
f = "{:<10} {:<20} {:<20}"
print("-"*80)
print(f.format("Name", "State value", "Derivative"))
f = "{: <10} {:< 20.13e} {:< 20.13e}"
print(f.format("membrane.V", c_membrane.V, c_membrane.d_v))
print(f.format("ina.m", c_ina.m, c_ina.d_m))
print(f.format("ina.h", c_ina.h, c_ina.d_h))
print(f.format("ina.j", c_ina.j, c_ina.d_j))
print(f.format("isi.d", c_isi.d, c_isi.d_d))
print(f.format("isi.f", c_isi.f, c_isi.d_f))
print(f.format("ix1.x1", c_ix1.x1, c_ix1.d_x1))
print(f.format("isi.Cai", c_isi.Cai, c_isi.d_cai))
#
# Test step function
#
def test_step():
""" Calculates and prints the initial derivatives """
init()
engine.update()
print_state()
#
# Pacing
#
class Protocol(object):
""" Holds an ordered set of ProtocolEvent objects """
def __init__(self):
super(Protocol, self).__init__()
self.head = None
def add(self, e):
""" Schedules an event """
if self.head is None:
self.head = e
return
if e.start < self.head.start:
e.next = self.head
self.head = e
return
f = self.head
while (f.next is not None) and (e.start > f.next.start):
f = f.next
e.next = f.next
f.next = e
def pop(self):
""" Returns the next event """
e = self.head
if self.head is not None:
self.head = self.head.next
return e
class ProtocolEvent(object):
def __init__(self, level, start, duration, period=0, multiplier=0):
super(ProtocolEvent, self).__init__()
self.level = float(level)
self.start = float(start)
self.duration = float(duration)
if self.duration <= 0:
raise Exception('Duration must be greater than zero')
self.period = float(period)
if self.period < 0:
raise Exception('multiplier must be zero or greater')
self.multiplier = int(multiplier)
if self.multiplier < 0:
raise Exception('Multiplier must be zero or greater')
if self.period == 0 and self.multiplier > 0:
raise Exception('Non-periodic event cannot occur more than once')
self.next = None
def pacing_protocol():
pacing = Protocol()
pacing.add(ProtocolEvent(1.0, 100.0, 0.5, 1000.0, 0))
return pacing
#
# Solver
#
def beat(stepSmall = 0.005, stepLarge = 0.01):
"""
Simulates a single beat
"""
tmin = 0
tmax = 1000
# Feedback
outInt = int(math.ceil(tmax / 10.0))
outPos = engine.time
outVal = 0
# Logging
logInt = 1
logPos = engine.time
log = []
# Stepsize
stepSize = stepSmall
hadPulse = False
vInit = c_membrane.V
# Pacing
pacing = pacing_protocol()
next = pacing.head
while next and next.start < tmin:
next = next.next
if next.start < tmin:
next = None
fire = None
fireDown = 0
stopTime = min(next.start, tmin + stepSize)
print('Starting integration with step sizes ' + str(stepSmall) + ' and ' + str(stepLarge) + '.')
while engine.time < tmax:
update(stopTime - engine.time)
engine.time = stopTime
# Event over
if (fire and engine.time >= fireDown):
engine.pace = 0
fire = None
# New event
if (next and engine.time >= next.start):
fire = next
next = next.next
engine.pace = fire.level
fireDown = fire.start + fire.duration
if fire.period > 0:
if fire.multiplier == 1:
fire.period = 0
else:
if fire.multiplier > 1:
fire.multiplier -=1
fire.start += fire.period
pacing.add(fire)
next = pacing.head
# User feedback
if engine.time >= outPos and outVal < 100:
print(str(outVal) + "%")
outVal += 10
outPos += outInt
# Logging
if engine.time >= logPos:
log.append((engine.time, state()))
logPos += logInt
# Step size update
if fire: # or c_membrane.V > -70:
if stepSize != stepSmall:
print("Small steps")
stepSize = stepSmall
else:
if stepSize != stepLarge:
print("Big steps")
stepSize = stepLarge
# Set next time
stopTime = engine.time + stepSize
if fire and fireDown < stopTime:
stopTime = fireDown
if next and next.start < stopTime:
stopTime = next.start
if logPos < stopTime:
stopTime = logPos
print("100% done")
print("t = " + str(engine.time))
print_state()
return log
#
# Run if loaded as main script
#
if __name__ == '__main__':
small = 0.005
large = 0.01
go = True
done = False
while go:
go = False
try:
init()
data = beat(small, large)
done = True
except ArithmeticError as e:
print('Arithmetic error occurred')
y = 'Continue with smaller stepsize? (y/n): '
try:
y = raw_input(y) # Python 2
except NameError:
y = input(y) # Python 3
if y.lower()[0:1] == 'y':
small /= 2
large /= 2
go = True
if done:
print('Showing result...')
x = []
y = []
for time, state in data:
x.append(time)
y.append(state[0])
import matplotlib.pyplot as py
plot = py.plot(x, y)
py.show()