"""
Benchmark models
"""
import circadapt
import circadapt.plot
from circadapt.model import Model
from circadapt.adapt import ModelAdapt
import matplotlib.pyplot as plt
import numpy as np
class Patch(Model):
def __init__(self,
solver: str = None,
path_to_circadapt: str = None,
model_state: dict = None,
):
if solver is None:
solver = 'forward_euler'
ModelAdapt.__init__(self)
Model.__init__(self,
solver,
path_to_circadapt=path_to_circadapt,
model_state=model_state,
)
def build(self):
"""
Build pre-afterload model with Patch.
This model is designed for illustrating the Patch isolated from the
circulation and for benchmarking the Patch module.
"""
self.add_component('LoadExperiment', 'W')
self.add_component('Patch', 'P', 'W')
def set_reference(self):
self['Patch']['dt'] = 0.1
self['Patch']['time_act'] = 0.3
self['Patch']['Sf_act'] = 120e3
self['Patch']['v_max'] = 7
self['LoadExperiment']['Am_dead'] = 0
self['LoadExperiment']['Am_ref_afterload'] = 1.3 * self['Patch']['Am_ref']
self['LoadExperiment']['T_afterload'] = 50
self.set('Model.W.P.l_si', 1.0, float)
# update assert test settings
self.assert_test_settings['l_si_t0'] = self.get('Model.W.P.l_si',
float)
def assert_1_beat(self):
"""
Test status after 1 beat, used in unit tests.
This function can be tested after creating the model with the
implemented reference and after 1 beat has been calculated. In case
the initial model state deviates from the implemented reference, test
settings might be changed accordingly.
Test for:
l_si == l_si_init on t=t0
l_se_norm == 1 on t=tend
"""
l_si_t0 = self['Patch']['l_si'][0, 0]
assert (l_si_t0 == self.assert_test_settings['l_si_t0'])
l_si_tend = self['Patch']['l_si'][-1, 0]
l_s_tend = self['Patch']['l_s'][-1, 0]
l_se = self['Patch']['l_se0'][0]
l_se_norm = (l_s_tend-l_si_tend) / l_se
assert np.abs(l_se_norm-1) < 1e-9
def plot(self, fig=None):
if fig is None:
fig = plt.figure(1, clear=True, figsize=(8, 6))
m = 2
n = 2
t = self.get('Solver.t') * 1e3
ax1 = fig.add_subplot(m, n, 1)
ax1.plot(t, self['Patch']['l_s'])
ax1.plot(t, self['Patch']['l_si'])
ax1.legend(['l_s', 'l_si'])
ax1.set_ylabel('Sarcomere length [$\mu m$]')
ax2 = fig.add_subplot(m, n, 2)
ax2.plot(t, self['Patch']['C'], label='C')
ax2.set_ylabel('C-curve [-]')
ax3 = fig.add_subplot(m, n, 3)
ax3.plot(t, self['Patch']['Sf']*1e-3, label='Sf')
ax3.plot(t, self['Patch']['Sf_ecm']*1e-3, label='SfEcm')
ax3.set_ylabel('Stress [kPa]')
ax3.legend()
ax4 = fig.add_subplot(m, n, 4)
ax4.plot(t, self['Patch']['T'], label='Tension Wall')
ax4.set_ylabel('T [$N/m$]')
# Plot design
for ax in [ax1, ax2, ax3, ax4]:
ax.spines[['right', 'top']].set_visible(False)
ax3.set_xlabel('Time [ms]')
ax4.set_xlabel('Time [ms]')
plt.tight_layout()
[docs]
class Patch(Model):
def __init__(self,
solver: str = None,
path_to_circadapt: str = None,
model_state: dict = None,
):
if solver is None:
solver = 'forward_euler'
ModelAdapt.__init__(self)
Model.__init__(self,
solver,
path_to_circadapt=path_to_circadapt,
model_state=model_state,
)
[docs]
def build(self):
"""
Build pre-afterload model with Patch.
This model is designed for illustrating the Patch isolated from the
circulation and for benchmarking the Patch module.
"""
self.add_component('LoadExperiment', 'W')
self.add_component('Patch', 'P', 'W')
[docs]
def set_reference(self):
self['Patch']['dt'] = 0.1
self['Patch']['time_act'] = 0.3
self['Patch']['Sf_act'] = 120e3
self['Patch']['v_max'] = 7
self['LoadExperiment']['Am_dead'] = 0
self['LoadExperiment']['Am_ref_afterload'] = 1.3 * self['Patch']['Am_ref']
self['LoadExperiment']['T_afterload'] = 50
self.set('Model.W.P.l_si', 1.0, float)
# update assert test settings
self.assert_test_settings['l_si_t0'] = self.get('Model.W.P.l_si',
float)
[docs]
def assert_1_beat(self):
"""
Test status after 1 beat, used in unit tests.
This function can be tested after creating the model with the
implemented reference and after 1 beat has been calculated. In case
the initial model state deviates from the implemented reference, test
settings might be changed accordingly.
Test for:
l_si == l_si_init on t=t0
l_se_norm == 1 on t=tend
"""
l_si_t0 = self['Patch']['l_si'][0, 0]
assert (l_si_t0 == self.assert_test_settings['l_si_t0'])
l_si_tend = self['Patch']['l_si'][-1, 0]
l_s_tend = self['Patch']['l_s'][-1, 0]
l_se = self['Patch']['l_se0'][0]
l_se_norm = (l_s_tend-l_si_tend) / l_se
assert np.abs(l_se_norm-1) < 1e-9
[docs]
def plot(self, fig=None):
if fig is None:
fig = plt.figure(1, clear=True, figsize=(8, 6))
m = 2
n = 2
t = self.get('Solver.t') * 1e3
ax1 = fig.add_subplot(m, n, 1)
ax1.plot(t, self['Patch']['l_s'])
ax1.plot(t, self['Patch']['l_si'])
ax1.legend(['l_s', 'l_si'])
ax1.set_ylabel('Sarcomere length [$\mu m$]')
ax2 = fig.add_subplot(m, n, 2)
ax2.plot(t, self['Patch']['C'], label='C')
ax2.set_ylabel('C-curve [-]')
ax3 = fig.add_subplot(m, n, 3)
ax3.plot(t, self['Patch']['Sf']*1e-3, label='Sf')
ax3.plot(t, self['Patch']['Sf_ecm']*1e-3, label='SfEcm')
ax3.set_ylabel('Stress [kPa]')
ax3.legend()
ax4 = fig.add_subplot(m, n, 4)
ax4.plot(t, self['Patch']['T'], label='Tension Wall')
ax4.set_ylabel('T [$N/m$]')
# Plot design
for ax in [ax1, ax2, ax3, ax4]:
ax.spines[['right', 'top']].set_visible(False)
ax3.set_xlabel('Time [ms]')
ax4.set_xlabel('Time [ms]')
plt.tight_layout()
[docs]
class Chamber2022(Model):
def __init__(self,
solver: str = None,
path_to_circadapt: str = None,
model_state: dict = None,
):
if solver is None:
solver = 'forward_euler'
Model.__init__(self,
solver,
path_to_circadapt=path_to_circadapt,
model_state=model_state,
)
[docs]
def build(self):
"""
Build pre-afterload model with Patch2022.
This model is designed for illustrating the Patch2022 isolated from the
circulation and for benchmarking the Patch2022 module.
"""
self.add_component('Chamber2022', 'C')
self.add_component('Patch2022', 'P', 'C.wC')
[docs]
class Valve2022(Model):
def __init__(self,
solver: str = None,
path_to_circadapt: str = None,
model_state: dict = None,
model_options: dict = {},
):
if solver is None:
solver = 'forward_euler'
# set options
self.model_options = model_options
self.model_options['t_on'] = self.model_options.get('t_on', 0.1)
self.model_options['t_off'] = self.model_options.get('t_off', 0.35)
self.model_options['p_prox_low'] = self.model_options.get('p_prox_low', 0.)
self.model_options['p_prox_high'] = self.model_options.get('p_prox_high', 2000.)
self.model_options['p_dist'] = self.model_options.get('p_dist', 1000.)
self.model_options['soft_closure'] = self.model_options.get('soft_closure', False)
self.model_options['q0'] = self.model_options.get('q0', 1e-3)
self.model_options['A_leak'] = self.model_options.get('A_leak', 1e-10)
self.model_options['A_open'] = self.model_options.get('A_open', 1e-3)
self.model_options['A_prox'] = self.model_options.get('A_prox', 2*self.model_options['A_open'])
self.model_options['A_dist'] = self.model_options.get('A_dist', 2*self.model_options['A_open'])
Model.__init__(self,
solver,
path_to_circadapt=path_to_circadapt,
model_state=model_state,
)
[docs]
def build(self):
"""
Build benchmark for Valve2022.
This model has a Valve2022 module with a prescribed proximal and
distal pressure using the NodePressure module.
"""
self.add_component('NodePressure', 'NodeProx')
self.add_component('NodePressure', 'NodeDist')
self.add_component('Valve2022', 'V')
self.set_component('V.prox', 'NodeProx')
self.set_component('V.dist', 'NodeDist')
[docs]
def set_reference(self):
# get parameters
t_on = self.model_options.get('t_on')
t_off = self.model_options.get('t_off')
p_prox_low = self.model_options.get('p_prox_low')
p_prox_high = self.model_options.get('p_prox_high')
p_dist = self.model_options.get('p_dist')
soft_closure = self.model_options.get('soft_closure')
q0 = self.model_options.get('q0')
A_leak = self.model_options.get('A_leak')
A_open = self.model_options.get('A_open')
A_prox = self.model_options.get('A_prox')
A_dist = self.model_options.get('A_dist')
# Set parameters to model
self['General']['t_cycle'] = 0.4
self.set('Model.V.soft_closure', soft_closure)
self.set('Model.V.q', q0)
self.set('Model.NodeDist.Pressure.set_time', 0.)
self.set('Model.NodeDist.Pressure.set_value', p_dist)
self.set('Model.NodeProx.Pressure.set_time', 0.)
self.set('Model.NodeProx.Pressure.set_value', p_prox_low)
self.set('Model.NodeProx.Pressure.set_time', t_on-2e-9)
self.set('Model.NodeProx.Pressure.set_value', 0.)
self.set('Model.NodeProx.Pressure.set_time', t_on-1e-9)
self.set('Model.NodeProx.Pressure.set_value', p_prox_high)
self.set('Model.NodeProx.Pressure.set_time', t_off-2e-9)
self.set('Model.NodeProx.Pressure.set_value', p_prox_high)
self.set('Model.NodeProx.Pressure.set_time', t_off-1e-9)
self.set('Model.NodeProx.Pressure.set_value', p_prox_low)
self.set('Model.V.A_leak', A_leak)
self.set('Model.V.A_open', A_open)
self.set('Model.NodeProx.A', A_prox)
self.set('Model.NodeDist.A', A_dist)
[docs]
class Valve(Model):
def __init__(self,
solver: str = None,
path_to_circadapt: str = None,
model_state: dict = None,
model_options: dict = {},
):
if solver is None:
solver = 'forward_euler'
# set options
self.model_options = model_options
self.model_options['t_on'] = self.model_options.get('t_on', 0.1)
self.model_options['t_off'] = self.model_options.get('t_off', 0.35)
self.model_options['p_prox_low'] = self.model_options.get('p_prox_low', 0.)
self.model_options['p_prox_high'] = self.model_options.get('p_prox_high', 2000.)
self.model_options['p_dist'] = self.model_options.get('p_dist', 1000.)
self.model_options['soft_closure'] = self.model_options.get('soft_closure', False)
self.model_options['q0'] = self.model_options.get('q0', 1e-3)
self.model_options['A_leak'] = self.model_options.get('A_leak', 1e-10)
self.model_options['A_open'] = self.model_options.get('A_open', 1e-3)
self.model_options['A_prox'] = self.model_options.get('A_prox', 2*self.model_options['A_open'])
self.model_options['A_dist'] = self.model_options.get('A_dist', 2*self.model_options['A_open'])
Model.__init__(self,
solver,
path_to_circadapt=path_to_circadapt,
model_state=model_state,
)
[docs]
def build(self):
"""
Build benchmark for Valve.
This model has a Valve module with a prescribed proximal and
distal pressure using the NodePressure module.
"""
self.add_component('NodePressure', 'NodeProx')
self.add_component('NodePressure', 'NodeDist')
self.add_component('Valve', 'V')
self.set_component('V.prox', 'NodeProx')
self.set_component('V.dist', 'NodeDist')
[docs]
def set_reference(self):
# get parameters
t_on = self.model_options.get('t_on')
t_off = self.model_options.get('t_off')
p_prox_low = self.model_options.get('p_prox_low')
p_prox_high = self.model_options.get('p_prox_high')
p_dist = self.model_options.get('p_dist')
soft_closure = self.model_options.get('soft_closure')
q0 = self.model_options.get('q0')
A_leak = self.model_options.get('A_leak')
A_open = self.model_options.get('A_open')
A_prox = self.model_options.get('A_prox')
A_dist = self.model_options.get('A_dist')
# Set parameters to model
self['General']['t_cycle'] = 0.4
self.set('Model.V.soft_closure', soft_closure)
self.set('Model.V.q', q0)
self.set('Model.NodeDist.Pressure.set_time', 0.)
self.set('Model.NodeDist.Pressure.set_value', p_dist)
self.set('Model.NodeProx.Pressure.set_time', 0.)
self.set('Model.NodeProx.Pressure.set_value', p_prox_low)
self.set('Model.NodeProx.Pressure.set_time', t_on-2e-9)
self.set('Model.NodeProx.Pressure.set_value', 0.)
self.set('Model.NodeProx.Pressure.set_time', t_on-1e-9)
self.set('Model.NodeProx.Pressure.set_value', p_prox_high)
self.set('Model.NodeProx.Pressure.set_time', t_off-2e-9)
self.set('Model.NodeProx.Pressure.set_value', p_prox_high)
self.set('Model.NodeProx.Pressure.set_time', t_off-1e-9)
self.set('Model.NodeProx.Pressure.set_value', p_prox_low)
self.set('Model.V.A_leak', A_leak)
self.set('Model.V.A_open', A_open)
self.set('Model.NodeProx.A', A_prox)
self.set('Model.NodeDist.A', A_dist)