1"""
2Tutorial CircAdapt VanOsta2022.
3
4November 2022, by Nick van Osta
5
6The goal of this tutorial is to understand the CircAdapt framework and to use
7the VanOsta2022 model. This tutorial assumes little to no knowledge about
8python. Therefore, basic python conventions and syntax will be discussed.
9This tutorial assumes the installation is followed as described on the wiki
10(https://wiki.circadapt.org/index.php?title=Circadapt_in_Python). This uses
11Python >3.9 installed with anaconda and editted in Spyder. Other ways are
12possible, but might not be in line with this tutorial.
13
14Content
15-------
16 1. Basics of python
17 2. Load the model
18 3. Plot global hemodynamics
19 4. Change parameters
20 5. Multipatch and local dynamics
21 6. Save and Load
22"""
23
24# Uncomment next lines if not installed
25# import sys
26# sys.path.append('../../../src/')
27
28import circadapt
29
30# Uncomment next lines if not installed
31# circadapt.DEFAULT_PATH_TO_CIRCADAPT = "../../../CircAdapt_Library/out/build/x64-Release/CircAdaptLib.dll"
32
33# %% 1. Basics of python
34print('1. Basics of python')
35
36# Always start the document with importing modules.
37# Numpy is used for mathematics, Matplotlib for plots. These are conventionally
38# imported as np and plt.
39import numpy as np
40import matplotlib.pyplot as plt
41
42# Only import the class we need in this tutorial.
43from circadapt.model import VanOsta2023
44# model = VanOsta2022()
45
46# alternatively, you can import the whole package, but it changes the way you
47# make the object.
48# import circadapt
49# model = circadapt.model.VanOsta2022()
50
51# Parameters types are automaticaly set or changed by the interpreter, but it
52# is good to create an integer when you need an integer and float when needed.
53i = 1 # integer
54f = 1. # float
55b = True # bool
56l = [1, 2, 3] # list
57d = {'a': 1, 'b': 2} # dictionary
58
59# get data from the list and dictionary
60first_item_of_list = l[0]
61item_from_dictionary = d['a']
62
63# the use of numpy is advised for more complex use and for calculation
64numpy_array = np.array(l)
65print('Find if array is 2: ', (numpy_array == 2))
66print('Multiply array with 2: ', (numpy_array * 2))
67
68# In spyder, you can place bullits. While debugging, ipython will stop at these
69# bullits. You can also press f9 to run a single line or selection and press
70# crtl+<enter> to run a block seperated by # %%
71
72# %% 2. load model
73print('\n 2. Load model. ')
74# Load predefined model with predefined parameterization
75# More information on this model can be found here:
76# https://wiki.circadapt.org/index.php?title=VanOsta2022
77model = VanOsta2023()
78
79# CircAdapt tries to follow the syntax of python and numpy as much as possible.
80# The object can be handled as a dictionary. Content can be printed in the
81# console, and printed on request in the ipython console.
82print('The result of printing the object gives information about the '
83 'components: ')
84print(model)
85
86# components can also be retrieved as a list of strings
87components = model.components
88print('Components of this model: ', components, '\n')
89
90# Similar to the object itself, components can be printed
91print('Patches in this model: ')
92print(model['Patch'])
93
94# Each component points to multiple c++ objects of that component type. The
95# objects can also be obtained using
96objects = model['Patch'].objects
97parameters = model['Patch'].parameters
98signals = model['Patch'].signals
99
100# Signals are not stored, so they are only available after running a beat.
101# Therefore the model should run. You can either run a number of beats, or run
102# until the model is hemodynamically stable.
103model.run(5)
104model.run(stable=True)
105
106# %% 3. Plot global hemodynamics
107# Here is an example code to plot the PV loop
108# First we open a figure. Assigning this figure to a variable is optional, but
109# is useful for design purposes.
110fig = plt.figure(1)
111
112# get volume and pressure of LV
113Vlv = model['Cavity']['V'][:, 'cLv']*1e6
114plv = model['Cavity']['p'][:, 'cLv']*7.5e-3
115
116# get volume and pressure of RV
117Vrv = model['Cavity']['V'][:, 'cRv']*1e6
118prv = model['Cavity']['p'][:, 'cRv']*7.5e-3
119
120# You can also use location names to get/set signals and parameters
121# For this, use only the last part of the full object name, e.g. cLv for
122# Model.Peri.TriSeg.cLv. You can get one signal or multiple signals
123Vlv = model['Cavity']['V'][:, 'cLv']*1e6
124Vrv = model['Cavity']['V'][:, 'cRv']*1e6
125pressure = model['Cavity']['p'][:, ['cLv', 'cRv']]*7.5e-3
126
127# you can split the two pressure signals into two parameters using the
128# following line. First transpose the pressure such that the first axis sets
129# the signals
130plv, prv = pressure.T
131
132# Now we plot the two lines.
133line1 = plt.plot(Vlv, plv, c='k', label='Lv')
134line2 = plt.plot(Vrv, prv, c='r', label='Rv')
135plt.ylabel('Pressure [mmHg]')
136plt.xlabel('Volume [mL]')
137plt.legend()
138
139# %% 4. Change parameters
140# Now reduce the contractility of all 3 ventricular walls
141model['Patch']['Sf_act'][2:] = 60e3
142
143# We can also explicitly change the ventricular wall patches
144model['Patch']['Sf_act'][['pLv0', 'pSv0', 'pRv0']] = 60e3
145
146# or do it one by one
147model['Patch']['Sf_act']['pLv0'] = 60e3
148model['Patch']['Sf_act']['pSv0'] = 60e3
149model['Patch']['Sf_act']['pRv0'] = 60e3
150
151# run the simulation and plot
152model.run(stable=True)
153plt.plot(model['Cavity']['V'][:, 'cLv']*1e6, model['Cavity']['p'][:, 'cLv']*7.5e-3, 'k--', label='Lv Reduced Sf_act')
154plt.plot(model['Cavity']['V'][:, 'cRv']*1e6, model['Cavity']['p'][:, 'cRv']*7.5e-3, 'r--', label='Rv Reduced Sf_act')
155plt.legend()
156
157
158# %% 5. Multipatch and local dynamics
159# Set up a new multipatch model and set an activation delay
160model_multipatch = VanOsta2023()
161
162# The number of patches is specified in the wall. Here, we set 12 Lv patches
163# and 6 Sv patches. Then, we change the dt in these patches.
164model_multipatch['Wall']['n_patch'][2:4] = [12, 6]
165model_multipatch['Patch']['dt'][2:14] = np.linspace(0, 0.01, 12)
166model_multipatch['Patch']['dt'][14:20] = np.linspace(0, 0.01, 6)
167
168# Run beats
169model_multipatch.run(stable=True)
170
171# Plot data
172fig = plt.figure(2)
173
174# In the first subplot, plot all volumes
175ax1 = plt.subplot(2, 2, 1)
176plt.plot(model_multipatch['Solver']['t']*1e3,
177 model_multipatch['Cavity']['V']*1e6,
178 )
179# in the second subplot, plot all pressures
180ax1 = plt.subplot(2, 2, 2)
181plt.plot(model_multipatch['Solver']['t']*1e3,
182 model_multipatch['Cavity']['p']*7.5e-3,
183 )
184# in the third subplot, plot all natural fiber strains.
185ax1 = plt.subplot(2, 2, 3)
186plt.plot(model_multipatch['Solver']['t']*1e3,
187 model_multipatch['Patch']['Ef'][:, 2:20],
188 )
189
190# %% 6. Save and Load
191# Simulations can be saved and loaded using the following code.
192model_reference = VanOsta2023()
193
194# use .npy extension in filename
195model_reference.save('reference.npy')
196
197# ander bestand
198model = VanOsta2023()
199model.load('reference.npy')
200
201# if you want to save and load a structure without writing it to a file, use
202# the follow lines. Note that signals are not filled, you have to run at least
203# 1 beat.
204data = model_reference.model_export()
205model_reference.model_import(data)