From 246db855c74ab820645d4ca02b91841bf4977dd6 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 10 Feb 2025 12:44:55 +0000 Subject: [PATCH 01/10] adding pre-commit to toml --- pyproject.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/pyproject.toml b/pyproject.toml index dff9a72..cdccc61 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,6 +26,7 @@ dependencies = [ [project.optional-dependencies] notebooks = ["ipykernel"] +dev = ["pre-commit"] [project.urls] "Homepage" = "https://github.com/alan-turing-institute/ModularCirc" From fa915a10b567454ef14c7808510ea4fd770c5de4 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 10 Feb 2025 12:49:53 +0000 Subject: [PATCH 02/10] adding balck, flake8 and isort packages. creating pre-commit-config.yaml file --- # .pre-commit-config.yaml | 9 +++++++++ pyproject.toml | 2 +- 2 files changed, 10 insertions(+), 1 deletion(-) create mode 100644 # .pre-commit-config.yaml diff --git a/# .pre-commit-config.yaml b/# .pre-commit-config.yaml new file mode 100644 index 0000000..78a66c1 --- /dev/null +++ b/# .pre-commit-config.yaml @@ -0,0 +1,9 @@ +# .pre-commit-config.yaml +repos: +- repo: https://github.com/pre-commit/pre-commit-hooks + rev: v4.0.1 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-yaml + - id: check-json \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index cdccc61..ac1795d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,7 +26,7 @@ dependencies = [ [project.optional-dependencies] notebooks = ["ipykernel"] -dev = ["pre-commit"] +dev = ["pre-commit", "black", "flake8", "isort"] [project.urls] "Homepage" = "https://github.com/alan-turing-institute/ModularCirc" From 8ebf9b459cdc113bf51d9ade1a039966c790f180 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 10 Feb 2025 13:42:23 +0000 Subject: [PATCH 03/10] adding packages to pre-commit file --- # .pre-commit-config.yaml | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/# .pre-commit-config.yaml b/# .pre-commit-config.yaml index 78a66c1..e427938 100644 --- a/# .pre-commit-config.yaml +++ b/# .pre-commit-config.yaml @@ -6,4 +6,19 @@ repos: - id: trailing-whitespace - id: end-of-file-fixer - id: check-yaml - - id: check-json \ No newline at end of file + - id: check-json + +- repo: https://github.com/psf/black + rev: 25.1.0 # Use the latest stable version + hooks: + - id: black + +- repo: https://github.com/PyCQA/flake8 + rev: 7.1.1 # Use the latest stable version + hooks: + - id: flake8 + +- repo: https://github.com/pre-commit/mirrors-isort + rev: 6.0.0 # Use the latest stable version + hooks: + - id: isort \ No newline at end of file From 5bb3a75e5a8bc32feae600c81dde7f18c6a90d31 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 10 Feb 2025 14:05:31 +0000 Subject: [PATCH 04/10] fixing formatting errors by adding an escape character to backslashes --- ModularCirc/Analysis/BaseAnalysis.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/ModularCirc/Analysis/BaseAnalysis.py b/ModularCirc/Analysis/BaseAnalysis.py index 1e7b1b3..99e9903 100644 --- a/ModularCirc/Analysis/BaseAnalysis.py +++ b/ModularCirc/Analysis/BaseAnalysis.py @@ -116,7 +116,7 @@ def plot_t_q_out(self, component:str, ax=None, time_units:str='s', volume_units: ) ax.set_title(component.upper() + ': PV loop') ax.set_xlabel(f'Time (${time_units}$)') - ax.set_ylabel(f'Flux (${volume_units}\cdot {time_units}^' + "{-1}$)") + ax.set_ylabel(f'Flux (${volume_units}\\cdot {time_units}^' + "{-1}$)") return ax def plot_fluxes(self, component:str, ax=None, time_units:str='s', volume_units:str='mL'): @@ -147,7 +147,7 @@ def plot_fluxes(self, component:str, ax=None, time_units:str='s', volume_units:s ) ax.set_title(f"{component.upper()}: Fluxes") ax.set_xlabel(f'Time (${time_units}$)') - ax.set_ylabel(f'Flux (${volume_units}\cdot {time_units}$)') + ax.set_ylabel(f'Flux (${volume_units}\\cdot {time_units}$)') ax.set_xlim(self.tsym[0], self.tsym[-1]) ax.legend() return ax @@ -160,11 +160,11 @@ def plot_valve_opening(self, component:str, ax=None, time_units:str='s'): self.model.commponents[component].PHI.values[self.tind], linestyle='-', linewidth=4, - label=f'{component} $\phi$' + label=f'{component} $\\phi$' ) - ax.set_title(f"{component.upper()}: $\Phi$") + ax.set_title(f"{component.upper()}: $\\Phi$") ax.set_xlabel(f'Time (${time_units}$)') - ax.set_ylabel(f'$\phi$') + ax.set_ylabel(f'$\\phi$') ax.set_xlim(self.tsym[0], self.tsym[-1]) ax.legend() return ax From bce851628c10bda6d2ea91df3fd2adf756ec0ab4 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 10 Feb 2025 14:14:11 +0000 Subject: [PATCH 05/10] updating the pre-commit yaml --- # .pre-commit-config.yaml | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/# .pre-commit-config.yaml b/# .pre-commit-config.yaml index e427938..5d20e14 100644 --- a/# .pre-commit-config.yaml +++ b/# .pre-commit-config.yaml @@ -7,18 +7,20 @@ repos: - id: end-of-file-fixer - id: check-yaml - id: check-json + - id: check-added-large-files + args: ['--maxkb=500'] # set the max file size to 500KB -- repo: https://github.com/psf/black - rev: 25.1.0 # Use the latest stable version - hooks: - - id: black +# - repo: https://github.com/psf/black +# rev: 25.1.0 # Use the latest stable version +# hooks: +# - id: black -- repo: https://github.com/PyCQA/flake8 - rev: 7.1.1 # Use the latest stable version - hooks: - - id: flake8 +# - repo: https://github.com/PyCQA/flake8 +# rev: 7.1.1 # Use the latest stable version +# hooks: +# - id: flake8 -- repo: https://github.com/pre-commit/mirrors-isort - rev: 6.0.0 # Use the latest stable version - hooks: - - id: isort \ No newline at end of file +# - repo: https://github.com/pre-commit/mirrors-isort +# rev: 6.0.0 # Use the latest stable version +# hooks: +# - id: isort \ No newline at end of file From 409de95ac6c62a3d12daca6b7318ffe4388e1e5c Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 10 Feb 2025 14:18:10 +0000 Subject: [PATCH 06/10] updating the precommit file name --- # .pre-commit-config.yaml => .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) rename # .pre-commit-config.yaml => .pre-commit-config.yaml (97%) diff --git a/# .pre-commit-config.yaml b/.pre-commit-config.yaml similarity index 97% rename from # .pre-commit-config.yaml rename to .pre-commit-config.yaml index 5d20e14..8fad403 100644 --- a/# .pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -23,4 +23,4 @@ repos: # - repo: https://github.com/pre-commit/mirrors-isort # rev: 6.0.0 # Use the latest stable version # hooks: -# - id: isort \ No newline at end of file +# - id: isort From b34202151e4321ba841a512fb02266761ec684f9 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 14 Apr 2025 15:55:07 +0100 Subject: [PATCH 07/10] pre-commits unleashed --- .github/workflows/ci.yml | 8 +- .gitignore | 2 +- README.md | 16 +- Tutorials/Tutorial_01/README.md | 4 +- Tutorials/Tutorial_01/circ_utils.py | 22 +-- Tutorials/Tutorial_01/parameters.json | 64 +++--- Tutorials/Tutorial_01/requirements.txt | 2 +- Tutorials/Tutorial_02/README.md | 6 +- Tutorials/Tutorial_03/parameters_01.json | 2 +- Tutorials/Tutorial_03/parameters_02.json | 2 +- src/ModularCirc/Analysis/BaseAnalysis.py | 124 ++++++------ src/ModularCirc/Components/ComponentBase.py | 29 ++- .../Components/HC_constant_elastance.py | 37 ++-- .../Components/HC_mixed_elastance.py | 40 ++-- src/ModularCirc/Components/R_component.py | 19 +- src/ModularCirc/Components/Rc_component.py | 51 +++-- src/ModularCirc/Components/Rlc_component.py | 22 +-- src/ModularCirc/Components/Valve_maynard.py | 26 +-- src/ModularCirc/Components/Valve_non_ideal.py | 17 +- .../Components/Valve_simple_bernoulli.py | 15 +- src/ModularCirc/Components/__init__.py | 2 +- src/ModularCirc/HelperRoutines.py | 72 +++---- src/ModularCirc/Models/KPat5MixedModel.py | 18 +- .../Models/KPat5MixedModel_parameters.py | 32 +-- .../Models/KorakianitisMaynardModel.py | 18 +- .../KorakianitisMaynardModel_parameters.py | 22 +-- .../Models/KorakianitisMixedMaynardModel.py | 18 +- ...orakianitisMixedMaynardModel_parameters.py | 22 +-- .../Models/KorakianitisMixedModel.py | 18 +- .../KorakianitisMixedModel_parameters.py | 24 +-- src/ModularCirc/Models/KorakianitisModel.py | 14 +- .../Models/KorakianitisModel_parameters.py | 26 +-- .../Models/MixedHeartMaynard4eWindkessel.py | 18 +- ...ixedHeartMaynard4eWindkessel_parameters.py | 28 +-- src/ModularCirc/Models/NaghaviModel.py | 82 ++++---- .../Models/NaghaviModelParameters.py | 38 ++-- src/ModularCirc/Models/OdeModel.py | 32 +-- src/ModularCirc/Models/ParametersObject.py | 18 +- src/ModularCirc/Models/README.md | 4 +- src/ModularCirc/Solver.py | 183 ++++++++---------- src/ModularCirc/StateVariable.py | 42 ++-- src/ModularCirc/Time.py | 28 +-- src/ModularCirc/_BatchRunner.py | 66 +++---- src/ModularCirc/__init__.py | 2 +- ...orakianitisMixedModel_expected_output.json | 2 +- .../NaghaviModel_expected_output.json | 2 +- tests/test_KorakianitisMixedModel.py | 44 ++--- tests/test_NaghaviModel.py | 42 ++-- tests/test_Solver.py | 52 ++--- 49 files changed, 726 insertions(+), 751 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index a9e3901..e30a2e2 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -38,15 +38,15 @@ jobs: environment-file: environment.yml python-version: '3.10' auto-activate-base: false - + - name: Install dependencies - run: | + run: | python -m pip install --upgrade pip pip install setuptools --upgrade pip install ./ - name: Run tests - run: | + run: | python -m unittest discover -s tests build-venv: @@ -76,4 +76,4 @@ jobs: - name: Run tests run: | source venv/bin/activate - python -m unittest discover -s tests \ No newline at end of file + python -m unittest discover -s tests diff --git a/.gitignore b/.gitignore index ff2a942..fe1de8e 100644 --- a/.gitignore +++ b/.gitignore @@ -171,4 +171,4 @@ cython_debug/ .pypirc # Ignore vscode settings -.vscode/ \ No newline at end of file +.vscode/ diff --git a/README.md b/README.md index 1882c08..6b6a1f8 100644 --- a/README.md +++ b/README.md @@ -7,7 +7,7 @@ The scope of this package is to provide a framework for building **0D models** and **simulating cardiovascular flow** and **mechanics**. Conceptually, the models can be split into three types of components: 1. **Heart chambers** 2. **Valves** -3. **Vessels** +3. **Vessels** ## Clone the ModularCirc GitHub repo locally @@ -26,7 +26,7 @@ Before installation of the ModularCirc package, please setup a virtual environme Install Conda from https://docs.conda.io/projects/conda/en/stable/user-guide/install/index.html -Run: +Run: ``` conda create --name @@ -37,7 +37,7 @@ Proceed to installing the ModularCirc package. ### Python virtual environment setup -Run `python3 -m venv venv`. This creates a virtual environment called `venv` in your base directory. +Run `python3 -m venv venv`. This creates a virtual environment called `venv` in your base directory. Activate the python environment: `source venv/bin/activate` @@ -45,7 +45,7 @@ Proceed to installing the ModularCirc package. ## Installation -To install the pip package: +To install the pip package: ```bash python -m pip install ModularCirc_LevanBokeria @@ -79,14 +79,14 @@ TEMPLATE_TIME_SETUP_DICT = { 'name' : 'TimeTest', 'ncycles' : 40, 'tcycle' : 1.0, - 'dt' : 0.001, + 'dt' : 0.001, 'export_min' : 1 } ``` Here, `ncycles` indicates the maximum number of heart cycles to run, before the simulation finishes. -If the simulation reaches steady state faster than that, the simulation will end provided the number of cycles is higher than `export_min`. -`tcycle` indicates the duration of the heart beat and `dt` represent the time step size used in the temporal discretization. -These measurements assume that time is measured in **seconds**. +If the simulation reaches steady state faster than that, the simulation will end provided the number of cycles is higher than `export_min`. +`tcycle` indicates the duration of the heart beat and `dt` represent the time step size used in the temporal discretization. +These measurements assume that time is measured in **seconds**. If the units used are different, ensure this is done consistently in line with other parameters. 4. Create an instance of the parameter object and used it to change the default values: diff --git a/Tutorials/Tutorial_01/README.md b/Tutorials/Tutorial_01/README.md index f1afe59..e7da96f 100644 --- a/Tutorials/Tutorial_01/README.md +++ b/Tutorials/Tutorial_01/README.md @@ -6,7 +6,7 @@ Involved in a larger project involving patients with pulmonary arterial hyperten Currently, while the monitors provide a lot of cardiac data, the usefulness of the data that they provide is limited, as it is difficult to interpret and use the information to make approprate changes to patient management. -The project focuses on the development of a digital twin to aid with interpretation of the cardiac data from these monitors. +The project focuses on the development of a digital twin to aid with interpretation of the cardiac data from these monitors. The sensitivity analysis aims to focus in on the parameters from the model that most affect pulmonary arterial pressure and cardiac output. @@ -28,4 +28,4 @@ Using a reduced set of parameters, you can more efficiently fit the model to pat 4) Complete K fold cross validation. 5) Retrain model on all the data with the new reduced number of components. 6) Use the reduced PCA results as output and the original parameter set you created as input for emulation - use [Autoemulate](https://github.com/alan-turing-institute/autoemulate) to find the best emulator for the data, this uses the "step3" notebook. - 7) Conduct a sensitivity analysis using the results from emulation using SAlib. \ No newline at end of file + 7) Conduct a sensitivity analysis using the results from emulation using SAlib. diff --git a/Tutorials/Tutorial_01/circ_utils.py b/Tutorials/Tutorial_01/circ_utils.py index d440782..416a013 100644 --- a/Tutorials/Tutorial_01/circ_utils.py +++ b/Tutorials/Tutorial_01/circ_utils.py @@ -26,8 +26,8 @@ def signal_get_pulse(signal, dt, num=100): Returns: _type_: _description_ """ - - ind = np.argmin(signal) + + ind = np.argmin(signal) ncycle = len(signal) new_signal = np.interp(np.linspace(0, ncycle, num), np.arange(ncycle), np.roll(signal, -ind)) new_dt = ncycle / (num - 1) * dt @@ -80,8 +80,8 @@ def run_case(row, output_path, N_cycles, dt): ) # replace the .. with the correct Class and inputs if applicable solver.setup( - suppress_output=True, - optimize_secondary_sv=False, + suppress_output=True, + optimize_secondary_sv=False, conv_cols=["p_ao"], method='LSODA' ) @@ -138,7 +138,7 @@ def simulation_loader(input_directory): file_series.sort_values('Index', inplace=True) file_series.set_index('Index', inplace=True) - + # Define a dataframe for the values collected from the simulation... signal_df = file_series.apply( lambda row: list(pd.read_csv(os.path.join(input_directory, row['file']), index_col='Index').to_numpy().reshape((-1))), @@ -147,7 +147,7 @@ def simulation_loader(input_directory): signal_df.rename(columns=template_columns, inplace=True) return signal_df - + dict_parameters_condensed_range = dict() dict_parameters_condensed_single = dict() @@ -166,7 +166,7 @@ def condense_dict_parameters(dict_param:dict, prev=""): dict_parameters_condensed_range[new_key] = tuple(np.array(r) * value) else: dict_parameters_condensed_single[new_key] = val[0] - return + return ######## DEFINED A FUNCTION TO PLOT THE VARIANCE @@ -186,7 +186,7 @@ def plot_variance(pca, width=8, dpi=100): cumulative_explained_variance = np.cumsum(explained_variance_ratio) axs[1].semilogy(grid, cumulative_explained_variance, "o-") axs[1].set( - xlabel="Component", title="% Cumulative Variance", + xlabel="Component", title="% Cumulative Variance", ) # Set up figure fig.set(figwidth=8, dpi=100) @@ -207,13 +207,13 @@ def scale_time_parameters_and_asign_to_components(df): # 800 ms in this case df['la.delay'] = df['la.delay'] * df['T'] / 800. - + df['la.t_tr'] = df['la.t_tr'] * df['T'] / 800. df['lv.t_tr'] = df['lv.t_tr'] * df['T'] / 800. - + df['la.tau'] = df['la.tau'] * df['T'] / 800. df['lv.tau'] = df['lv.tau'] * df['T'] / 800. df['la.t_max'] = df['la.t_max'] * df['T'] / 800. df['lv.t_max'] = df['lv.t_max'] * df['T'] / 800. - return \ No newline at end of file + return diff --git a/Tutorials/Tutorial_01/parameters.json b/Tutorials/Tutorial_01/parameters.json index e371504..f622cfb 100644 --- a/Tutorials/Tutorial_01/parameters.json +++ b/Tutorials/Tutorial_01/parameters.json @@ -1,35 +1,35 @@ { "VESSELS" :{ - + "ao" : { - "r" :[240, + "r" :[240, [0.5, 1.5] - ], - "c" :[0.3, + ], + "c" :[0.3, [0.5, 1.5] - ], + ], "l" :[0], "v_ref":[100], "v" :[130] }, "art" : { - "r" :[1125, - [0.5, 1.5] + "r" :[1125, + [0.5, 1.5] + ], + "c" :[3, + [0.5, 1.5] ], - "c" :[3, - [0.5, 1.5] - ], "l" :[0], "v_ref":[900], "v" :[1092] }, "ven" : { - "r" :[9, + "r" :[9, [0.5, 1.5] - ], - "c" :[133.3, + ], + "c" :[133.3, [0.5, 1.5] - ], + ], "l" :[0.0], "v_ref":[2800], "v" :[3780] @@ -37,47 +37,47 @@ }, "VALVES" : { "av": { - "r" :[6, + "r" :[6, [0.5, 1.5] ] }, "mv": { - "r" :[4.1, + "r" :[4.1, [0.5, 1.5] ] } }, "CHAMBERS": { "la" : { - "E_pas":[0.44, + "E_pas":[0.44, [0.5, 1.5] - ], - "E_act":[0.45, + ], + "E_act":[0.45, [0.5, 1.5] - ], - "v_ref":[10, + ], + "v_ref":[10, [0.5, 1.5] ], - "k_pas":[0.05, + "k_pas":[0.05, [0.333, 1.5] ], - "v" :[93], + "v" :[93], "delay":[150], "t_tr" :[225], "tau" :[25], "t_max":[150] }, "lv" : { - "E_pas":[1.0, + "E_pas":[1.0, [0.5, 1.5] - ], - "E_act":[3, + ], + "E_act":[3, [0.5, 1.5] - ], - "v_ref":[10.0, + ], + "v_ref":[10.0, [0.5, 1.5] ], - "k_pas":[0.03, + "k_pas":[0.03, [0.333, 1.5] ], "delay":[0], @@ -86,7 +86,7 @@ "t_max":[280] } }, - "T" :[[800], + "T" :[[800], [0.375, 1.5] - ] -} \ No newline at end of file + ] +} diff --git a/Tutorials/Tutorial_01/requirements.txt b/Tutorials/Tutorial_01/requirements.txt index 3d4fb80..96d61a3 100644 --- a/Tutorials/Tutorial_01/requirements.txt +++ b/Tutorials/Tutorial_01/requirements.txt @@ -4,4 +4,4 @@ tqdm ipykernel jupyter autoemulate -SALib \ No newline at end of file +SALib diff --git a/Tutorials/Tutorial_02/README.md b/Tutorials/Tutorial_02/README.md index 5322243..e45ee58 100644 --- a/Tutorials/Tutorial_02/README.md +++ b/Tutorials/Tutorial_02/README.md @@ -1,6 +1,6 @@ # A short overview of Modularcirc -**ModularCirc** is intended to be a versatile platform of 0D models of the cardiovascular system. +**ModularCirc** is intended to be a versatile platform of 0D models of the cardiovascular system. It is users are envisioned to have different levels of expertise/experience coding. ![Main figure](./Figures/Slides.png) @@ -15,7 +15,5 @@ This will cover: 4) How to generate quick plots of the data/retrieve basic values from the model - ### *Part II*: 1) What is the structure of a model and supporting classes - 2) How to retrieve predefined components + 2) How to retrieve predefined components 3) How to modify existing models with new components? - - \ No newline at end of file diff --git a/Tutorials/Tutorial_03/parameters_01.json b/Tutorials/Tutorial_03/parameters_01.json index 9ca1811..20bcee3 100644 --- a/Tutorials/Tutorial_03/parameters_01.json +++ b/Tutorials/Tutorial_03/parameters_01.json @@ -192,4 +192,4 @@ 1.3 ] ] -} \ No newline at end of file +} diff --git a/Tutorials/Tutorial_03/parameters_02.json b/Tutorials/Tutorial_03/parameters_02.json index e39482b..de22313 100644 --- a/Tutorials/Tutorial_03/parameters_02.json +++ b/Tutorials/Tutorial_03/parameters_02.json @@ -381,4 +381,4 @@ 4 ] ] -} \ No newline at end of file +} diff --git a/src/ModularCirc/Analysis/BaseAnalysis.py b/src/ModularCirc/Analysis/BaseAnalysis.py index ff85823..e7c90d7 100644 --- a/src/ModularCirc/Analysis/BaseAnalysis.py +++ b/src/ModularCirc/Analysis/BaseAnalysis.py @@ -10,63 +10,63 @@ class ValveData(): def __init__(self, name) -> None: - self._name = name - + self._name = name + def set_opening_closing(self, open, closed): self._open = open self._closed = closed - + def __repr__(self) -> str: - return f"Valve {self._name}: \n" + f" - opening ind: {self._open} \n" + f" - closing ind: {self._closed} \n" - + return f"Valve {self._name}: \n" + f" - opening ind: {self._open} \n" + f" - closing ind: {self._closed} \n" + @property def open(self): return self._open - + @property def closed(self): return self._closed - - + + class VentricleData(): def __init__(self, name:str, volume_unit:str='ml') -> None: self._name = name self._vu = volume_unit - + def set_volumes(self, edv, esv): self._edv = edv self._esv = esv - + def __repr__(self) -> str: - return (f"Ventricle {self._name}: \n" + - f" - EDV: {self._edv:.2e} {self._vu}\n" + + return (f"Ventricle {self._name}: \n" + + f" - EDV: {self._edv:.2e} {self._vu}\n" + f" - ESV: {self._esv:.2e} {self._vu}" ) - + @property def edv(self): return self._edv - + @property def esv(self): return self._esv - + class BaseAnalysis(): def __init__(self, model:OdeModel=None) -> None: self.model = model - + self.valves= dict() self.ventricles = dict() - + self.tind = np.arange(start=self.model.time_object.n_t-(self.model.time_object.n_c-1) * self.model.time_object.export_min, stop =self.model.time_object.n_t) self.tsym = self.model.time_object._sym_t.values[self.tind] self.tsym = self.tsym - self.tsym[0] - + def plot_t_v(self, component:str, ax=None, time_units:str='s', volume_units:str='mL', linestyle='-'): if ax is None: _, ax = plt.subplots(figsize=(5,5)) - ax.plot(self.tsym, + ax.plot(self.tsym, self.model.components[component].V.values[self.tind], linewidth=4, label=component, @@ -76,11 +76,11 @@ def plot_t_v(self, component:str, ax=None, time_units:str='s', volume_units:str= ax.set_ylabel(f'Volume (${volume_units}$)') ax.set_xlim(self.tsym[0], self.tsym[-1]) return ax - + def plot_t_p(self, component:str, ax=None, time_units:str='s', pressure_units:str='mmHg', linestyle='-'): if ax is None: _, ax = plt.subplots(figsize=(5,5)) - ax.plot(self.tsym, + ax.plot(self.tsym, self.model.components[component].P.values[self.tind], linewidth=4, label=component, @@ -90,7 +90,7 @@ def plot_t_p(self, component:str, ax=None, time_units:str='s', pressure_units:st ax.set_ylabel(f'Pressure (${pressure_units}$)') ax.set_xlim(self.tsym[0], self.tsym[-1]) return ax - + def plot_p_v_loop(self, component, ax=None, volume_units:str='mL', pressure_units:str="mmHg"): if ax is None: _, ax = plt.subplots(figsize=(5,5)) @@ -103,12 +103,12 @@ def plot_p_v_loop(self, component, ax=None, volume_units:str='mL', pressure_unit ax.set_xlabel(f'Volume (${volume_units}$)') ax.set_ylabel(f'Pressure (${pressure_units}$)') return ax - + def plot_t_q_out(self, component:str, ax=None, time_units:str='s', volume_units:str='mL', linestyle='-'): if ax is None: _, ax = plt.subplots(figsize=(5,5)) ax.plot( - self.tsym, + self.tsym, self.model.components[component].Q_o.values[self.tind], linewidth=4, linestyle=linestyle, @@ -124,34 +124,34 @@ def plot_fluxes(self, component:str, ax=None, time_units:str='s', volume_units:s _, ax = plt.subplots(figsize=(5,5)) ax.plot( self.tsym, - self.model.components[component].Q_i.values[self.tind] - + self.model.components[component].Q_i.values[self.tind] - self.model.components[component].Q_o.values[self.tind], linestyle='-', linewidth=4, alpha=0.6, label=f'{component} $dV/dt$' - ) + ) ax.plot( self.tsym, self.model.components[component].Q_i.values[self.tind], linestyle=':', linewidth=4, label=f'{component} $Q_i$' - ) + ) ax.plot( self.tsym, self.model.components[component].Q_o.values[self.tind], linestyle=':', linewidth=4, label=f'{component} $Q_o$' - ) - ax.set_title(f"{component.upper()}: Fluxes") + ) + ax.set_title(f"{component.upper()}: Fluxes") ax.set_xlabel(f'Time (${time_units}$)') - ax.set_ylabel(f'Flux (${volume_units}\\cdot {time_units}$)') + ax.set_ylabel(f'Flux (${volume_units}\\cdot {time_units}$)') ax.set_xlim(self.tsym[0], self.tsym[-1]) - ax.legend() + ax.legend() return ax - + def plot_valve_opening(self, component:str, ax=None, time_units:str='s'): if ax is None: _, ax = plt.subplots(figsize=(5,5)) @@ -162,21 +162,21 @@ def plot_valve_opening(self, component:str, ax=None, time_units:str='s'): linewidth=4, label=f'{component} $\\phi$' ) - ax.set_title(f"{component.upper()}: $\\Phi$") + ax.set_title(f"{component.upper()}: $\\Phi$") ax.set_xlabel(f'Time (${time_units}$)') - ax.set_ylabel(f'$\\phi$') + ax.set_ylabel(f'$\\phi$') ax.set_xlim(self.tsym[0], self.tsym[-1]) - ax.legend() + ax.legend() return ax - + def compute_opening_closing_valve(self, component:str, shift:float=0.0): """ Method using for computing a valve opening and closing time indexes. Output values are stored in: - `self.valve[component].open` - `self.valve[component].closed` - - ## Inputs + + ## Inputs component : str name of valve that is being measured shift : float @@ -185,7 +185,7 @@ def compute_opening_closing_valve(self, component:str, shift:float=0.0): valve = self.model.components[component] nshift= int(shift/ self.model.time_object.dt) self.valves[component] = ValveData(component) - + if not hasattr(valve, 'PHI'): pi = valve.P_i.values[self.tind] po = valve.P_o.values[self.tind] @@ -195,13 +195,13 @@ def compute_opening_closing_valve(self, component:str, shift:float=0.0): min_phi = np.min(phi) is_open = (phi - min_phi) > 1.0e-2 is_open_shifted = np.roll(is_open, -nshift) - + ind = np.arange(len(is_open))[is_open_shifted] self.valves[component].set_opening_closing(open = ind[0], closed= ind[-1]) return - - + + def compute_ventricle_volume_limits(self, component:str, vic:int, voc:int)->None: """ Method for computing the end diastolic and systolic volumes of ventricles. Results are stored in `self.ventricles[component].edv` and `self.ventricles[component].esv`. @@ -216,41 +216,41 @@ def compute_ventricle_volume_limits(self, component:str, vic:int, voc:int)->None """ ventricle = self.model.components[component] volume = ventricle.V.values[self.tind] - + self.ventricles[component] = VentricleData(component) self.ventricles[component].set_volumes(edv=volume[vic], esv=volume[voc]) return - - + + def compute_cardiac_output(self, component:str)->float: """ Method for estimating the cardiac output. - + ## Inputs component : str name of the component (ideally the aortic valve) for which we are computing the cardiac output - + ## Outputs CO : float - the cardiac ouput value + the cardiac ouput value """ valve = self.model.components[component] dt = self.model.time_object.dt T = self.model.time_object.tcycle / 60.0 - + q = valve.Q_i.values[self.tind] self.CO = q.sum() * dt / T / self.model.time_object.export_min - + return self.CO - + def compute_ejection_fraction(self, component:str)->float: """ Method for estimating the cardiac output. - + ## Inputs component : str name of the component (ideally one of the ventricles) for which we are computing Ejection Fraction - + ## Outputs CO : float the ejection fraction @@ -260,8 +260,8 @@ def compute_ejection_fraction(self, component:str)->float: esv = ventricle.V.values[self.tind].min() self.EF = (edv - esv) / edv return self.EF - - + + def compute_artery_pressure_range(self, component:str)->tuple[float]: """Method for computing the range of pressures in artery @@ -275,27 +275,27 @@ def compute_artery_pressure_range(self, component:str)->tuple[float]: c = self.model.components[component] p = c.P_i.values[self.tind] return (np.min(p), np.max(p)) - + def compute_end_systolic_pressure(self, component:str, upstream:str): """" Method for computing the end systolic pressure in the aorta and pulmonary artery. - + ## Inputs component : str name of the component for which we are computing the end systolic pressure upstream : str name of the upstream component which is used as reference point for when the valve in between components closes - + ## Ouputs - ESP : float + ESP : float vessel end systolic pressure """ c = self.model.components[component] p = c.P_i.values[self.tind] - + uc = self.model.components[upstream] up = uc.P_i.values[self.tind] - - flag = up > p + + flag = up > p ind = np.arange(len(flag))[flag] return p[ind[-1]] diff --git a/src/ModularCirc/Components/ComponentBase.py b/src/ModularCirc/Components/ComponentBase.py index 22ca65a..d154d26 100644 --- a/src/ModularCirc/Components/ComponentBase.py +++ b/src/ModularCirc/Components/ComponentBase.py @@ -10,13 +10,13 @@ def __init__(self, ) -> None: self._name = name self._to = time_object - + self._P_i = StateVariable(name=name+'_P_i', timeobj=time_object) self._Q_i = StateVariable(name=name+'_Q_i', timeobj=time_object) self._P_o = StateVariable(name=name+'_P_o', timeobj=time_object) self._Q_o = StateVariable(name=name+'_Q_o', timeobj=time_object) self._V = StateVariable(name=name+'_V' , timeobj=time_object) - + if v is not None: self.v0 = v self.V.loc[0] = v @@ -24,7 +24,7 @@ def __init__(self, self.v0 = None self.p0 = p return - + def __repr__(self) -> str: var = (f" Component {self._name}" + '\n' + f" - P_i " + str(self._P_i) + '\n' + @@ -33,39 +33,39 @@ def __repr__(self) -> str: f" - Q_o " + str(self._Q_o) + '\n' ) return var - + @property def P_i(self): return self._P_i._u - + @property def P_o(self): return self._P_o._u - + @property def Q_i(self): return self._Q_i._u - + @property def Q_o(self): return self._Q_o._u - + @property def V(self): return self._V._u - + def make_unique_io_state_variable(self, q_flag:bool=False, p_flag:bool=True) -> None: - if q_flag: + if q_flag: self._Q_o = self._Q_i self._Q_o.set_name(self._name + '_Q') - if p_flag: + if p_flag: self._P_o = self._P_i self._P_o.set_name(self._name + '_P') - return - + return + def setup(self) -> None: raise Exception("This is a template class only.") - + def __del__(self): # print('Blah') if hasattr(self, '_P_i'): @@ -78,4 +78,3 @@ def __del__(self): del self._Q_o if hasattr(self, '_V'): del self._V - \ No newline at end of file diff --git a/src/ModularCirc/Components/HC_constant_elastance.py b/src/ModularCirc/Components/HC_constant_elastance.py index 90d7f3e..e6a7bc8 100644 --- a/src/ModularCirc/Components/HC_constant_elastance.py +++ b/src/ModularCirc/Components/HC_constant_elastance.py @@ -8,10 +8,10 @@ import numpy as np class HC_constant_elastance(ComponentBase): - def __init__(self, + def __init__(self, name:str, time_object: TimeClass, - E_pas: float, + E_pas: float, E_act: float, v_ref: float, v : float = None, @@ -24,62 +24,62 @@ def __init__(self, self.E_act = E_act self.v_ref = v_ref self.eps = 1.0e-3 - + def af_parameterised(t): return af(time_shift(t, kwargs['delay'], time_object.tcycle) , **kwargs) self._af = af_parameterised - + self.make_unique_io_state_variable(p_flag=True, q_flag=False) - - @property + + @property def P(self): return self._P_i._u - + def comp_E(self, t:float) -> float: return self._af(t) * self.E_act + (1.0 - self._af(t)) * self.E_pas - + def comp_dEdt(self, t:float) -> float: return (self.comp_E(t + self.eps) - self.comp_E(t - self.eps)) / 2.0 / self.eps - + def comp_p(self, t:float, v:float=None, y:np.ndarray[float]=None) ->float: e = self.comp_E(t) if y is not None: v = y return e * (v - self.v_ref) - + def comp_v(self, t:float=None, p:float=None, y:np.ndarray[float]=None)->float: e = self.comp_E(t) if y is not None: p = y return p / e + self.v_ref - + def comp_dpdt(self, t:float=None, V:float=None, q_i:float=None, q_o:float=None, y:np.ndarray[float]=None) -> float: if y is not None: V, q_i, q_o, = y[:3] dEdt = self.comp_dEdt(t) e = self.comp_E(t) return dEdt * (V - self.v_ref) + e * (q_i - q_o) - + def setup(self) -> None: E_pas = self.E_pas E_act = self.E_act v_ref = self.v_ref eps = self.eps af = self._af - + comp_E = lambda t: af(t) * E_act + (1.0 - af(t)) * E_pas comp_dEdt = lambda t: (comp_E(t + eps) - comp_E(t - eps)) / 2.0 / eps comp_p = lambda t, y: comp_E(t) * (y - v_ref) comp_v = lambda t, y: y / comp_E(t) + v_ref comp_dpdt = lambda t, y: comp_dEdt(t) * (y[0] - v_ref) + comp_E(t) * (y[1] - y[2]) - + self._V.set_dudt_func(chamber_volume_rate_change, function_name='chamber_volume_rate_change') - self._V.set_inputs(pd.Series({'q_in':self._Q_i.name, + self._V.set_inputs(pd.Series({'q_in':self._Q_i.name, 'q_out':self._Q_o.name})) - self._P_i.set_dudt_func(comp_dpdt, function_name='self.comp_dpdt') - self._P_i.set_inputs(pd.Series({'V':self._V.name, - 'q_i':self._Q_i.name, + self._P_i.set_dudt_func(comp_dpdt, function_name='self.comp_dpdt') + self._P_i.set_inputs(pd.Series({'V':self._V.name, + 'q_i':self._Q_i.name, 'q_o':self._Q_o.name})) if self.p0 is None or self.p0 is np.NaN: self._P_i.set_i_func(comp_p, function_name='self.comp_p') @@ -91,4 +91,3 @@ def setup(self) -> None: self._V.set_i_inputs(pd.Series({'p':self._P_i.name})) if (self.v0 is None or self.v0 is np.NaN) and (self.p0 is None or self.p0 is np.NaN): raise Exception("Solver needs at least the initial volume or pressure to be defined!") - \ No newline at end of file diff --git a/src/ModularCirc/Components/HC_mixed_elastance.py b/src/ModularCirc/Components/HC_mixed_elastance.py index ab4a27d..1382706 100644 --- a/src/ModularCirc/Components/HC_mixed_elastance.py +++ b/src/ModularCirc/Components/HC_mixed_elastance.py @@ -46,13 +46,13 @@ def gen_total_dpdt(active_p, passive_p, _af, active_dpdt, passive_dpdt): def total_dpdt(t, y): _af_t = _af(t) _d_af_dt = _af(t, dt=True) - return (_d_af_dt *(active_p(y[0]) - passive_p(y[0])) + - _af_t * active_dpdt( y[1], y[2]) + + return (_d_af_dt *(active_p(y[0]) - passive_p(y[0])) + + _af_t * active_dpdt( y[1], y[2]) + (1. -_af_t) * passive_dpdt(y[0], y[1], y[2])) return total_dpdt def gen_comp_v(E_pas, v_ref, k_pas): - def comp_v(t, y): + def comp_v(t, y): return v_ref + np.log(y[0] / E_pas + 1.0) / k_pas return comp_v @@ -63,16 +63,16 @@ def time_shifter(t): def gen__af(af, time_shifter, kwargs): varnames = [name for name in af.__code__.co_varnames if name != 'coeff' and name != 't'] - kwargs2 = {key: val for key,val in kwargs.items() if key in varnames} + kwargs2 = {key: val for key,val in kwargs.items() if key in varnames} def _af(t, dt=False): return af(time_shifter(t), dt=dt, **kwargs2) return _af class HC_mixed_elastance(ComponentBase): - def __init__(self, + def __init__(self, name:str, time_object: TimeClass, - E_pas: float, + E_pas: float, E_act: float, k_pas: float, v_ref: float, @@ -89,13 +89,13 @@ def __init__(self, self.eps = 1.0e-3 self.kwargs = kwargs self.af = af - + self.make_unique_io_state_variable(p_flag=True, q_flag=False) - - @property + + @property def P(self): return self._P_i._u - + def setup(self) -> None: E_pas = self.E_pas k_pas = self.k_pas @@ -105,27 +105,27 @@ def setup(self) -> None: kwargs= self.kwargs T = self._to.tcycle af = self.af - + time_shifter = gen_time_shifter(delay_=kwargs['delay'], T=T) _af = gen__af(af=af, time_shifter=time_shifter, kwargs=kwargs) - + active_p = gen_active_p(E_act=E_act, v_ref=v_ref) active_dpdt = gen_active_dpdt(E_act=E_act) passive_p = gen_passive_p(E_pas=E_pas, k_pas=k_pas, v_ref=v_ref) passive_dpdt = gen_passive_dpdt(E_pas=E_pas, k_pas=k_pas, v_ref=v_ref) total_p = gen_total_p(_af=_af, active_p=active_p, passive_p=passive_p) - total_dpdt = gen_total_dpdt(active_p=active_p, passive_p=passive_p, + total_dpdt = gen_total_dpdt(active_p=active_p, passive_p=passive_p, _af=_af, active_dpdt=active_dpdt, passive_dpdt=passive_dpdt) comp_v = gen_comp_v(E_pas=E_pas, v_ref=v_ref, k_pas=k_pas) - + self._V.set_dudt_func(chamber_volume_rate_change, function_name='chamber_volume_rate_change') - self._V.set_inputs(pd.Series({'q_in' :self._Q_i.name, + self._V.set_inputs(pd.Series({'q_in' :self._Q_i.name, 'q_out':self._Q_o.name})) - - self._P_i.set_dudt_func(total_dpdt, function_name='total_dpdt') - self._P_i.set_inputs(pd.Series({'v' :self._V.name, - 'q_i':self._Q_i.name, + + self._P_i.set_dudt_func(total_dpdt, function_name='total_dpdt') + self._P_i.set_inputs(pd.Series({'v' :self._V.name, + 'q_i':self._Q_i.name, 'q_o':self._Q_o.name})) if self.p0 is None or self.p0 is np.NaN: self._P_i.set_i_func(total_p, function_name='total_p') @@ -136,4 +136,4 @@ def setup(self) -> None: self._V.set_i_func(comp_v, function_name='comp_v') self._V.set_i_inputs(pd.Series({'p':self._P_i.name})) if (self.v0 is None or self.v0 is np.NaN) and (self.p0 is None or self.p0 is np.NaN): - raise Exception("Solver needs at least the initial volume or pressure to be defined!") \ No newline at end of file + raise Exception("Solver needs at least the initial volume or pressure to be defined!") diff --git a/src/ModularCirc/Components/R_component.py b/src/ModularCirc/Components/R_component.py index 1d6d522..a8fed76 100644 --- a/src/ModularCirc/Components/R_component.py +++ b/src/ModularCirc/Components/R_component.py @@ -1,31 +1,30 @@ from .ComponentBase import ComponentBase from ..Time import TimeClass -from ..HelperRoutines import resistor_upstream_pressure +from ..HelperRoutines import resistor_upstream_pressure import pandas as pd class R_component(ComponentBase): - def __init__(self, + def __init__(self, name:str, time_object: TimeClass, - r:float, + r:float, ) -> None: super().__init__(name=name, time_object=time_object) # allow for pressure gradient but not for flow - self.make_unique_io_state_variable(q_flag=True, p_flag=False) + self.make_unique_io_state_variable(q_flag=True, p_flag=False) # setting the resistance value self.R = r - - @property + + @property def P(self): return self._P_i._u - + def p_i_u_func(self, t, y): return resistor_upstream_pressure(t, y=y, r=self.R) - + def setup(self) -> None: r=self.R self._P_i.set_u_func(lambda t, y: resistor_upstream_pressure(t, y, r=r), function_name='resistor_upstream_pressure') - self._P_i.set_inputs(pd.Series({'q_in' :self._Q_i.name, + self._P_i.set_inputs(pd.Series({'q_in' :self._Q_i.name, 'p_out':self._P_o.name})) - \ No newline at end of file diff --git a/src/ModularCirc/Components/Rc_component.py b/src/ModularCirc/Components/Rc_component.py index f9a79b2..118b111 100644 --- a/src/ModularCirc/Components/Rc_component.py +++ b/src/ModularCirc/Components/Rc_component.py @@ -25,13 +25,13 @@ def q_o_u_func(t, y): return q_o_u_func class Rc_component(ComponentBase): - def __init__(self, - name:str, + def __init__(self, + name:str, time_object:TimeClass, - r:float, - c:float, + r:float, + c:float, v_ref:float, - v:float=None, + v:float=None, p:float=None, ) -> None: # super().__init__(time_object, main_var) @@ -39,61 +39,60 @@ def __init__(self, self.R = r self.C = c self.V_ref = v_ref - + if p is not None: self.p0 = p self.P_i.loc[0] = p else: self.p0 = None - - @property + + @property def P(self): return self._P_i._u - + def q_o_u_func(self, t, y): return resistor_model_flow(t=t, y=y, r=self.R) - + def p_i_dudt_func(self, t, y): return grounded_capacitor_model_dpdt(t, y=y, c=self.C) - + def p_i_i_func(self, t, y): return grounded_capacitor_model_pressure(t, y=y, v_ref=self.V_ref, c=self.C) - + def v_i_func(self, t, y): return grounded_capacitor_model_volume(t, y=y, v_ref=self.V_ref, c=self.C) - + def setup(self) -> None: r = self.R v_ref = self.V_ref c = self.C - # Set the dudt function for the input pressure state variable - self._P_i.set_dudt_func(gen_p_i_dudt_func(C=c), + # Set the dudt function for the input pressure state variable + self._P_i.set_dudt_func(gen_p_i_dudt_func(C=c), function_name='grounded_capacitor_model_dpdt') # Set the mapping betwen the local input names and the global names of the state variables - self._P_i.set_inputs(pd.Series({'q_in' :self._Q_i.name, + self._P_i.set_inputs(pd.Series({'q_in' :self._Q_i.name, 'q_out':self._Q_o.name})) if self.p0 is None or self.p0 is np.NaN: # Set the initialization function for the input pressure state variable - self._P_i.set_i_func(gen_p_i_i_func(v_ref=v_ref, c=c), + self._P_i.set_i_func(gen_p_i_i_func(v_ref=v_ref, c=c), function_name='grounded_capacitor_model_pressure') self._P_i.set_i_inputs(pd.Series({'v':self._V.name})) else: self.P_i.loc[0] = self.p0 # Set the function for computing the flows based on the current pressure values at the nodes of the componet - self._Q_o.set_u_func(gen_q_o_u_func(r=r), + self._Q_o.set_u_func(gen_q_o_u_func(r=r), function_name='resistor_model_flow' ) - self._Q_o.set_inputs(pd.Series({'p_in':self._P_i.name, + self._Q_o.set_inputs(pd.Series({'p_in':self._P_i.name, 'p_out':self._P_o.name})) # Set the dudt function for the compartment volume - self._V.set_dudt_func(chamber_volume_rate_change, + self._V.set_dudt_func(chamber_volume_rate_change, function_name='chamber_volume_rate_change') - self._V.set_inputs(pd.Series({'q_in':self._Q_i.name, + self._V.set_inputs(pd.Series({'q_in':self._Q_i.name, 'q_out':self._Q_o.name})) if self.v0 is None or self.v0 is np.NaN: - # Set the initialization function for the input volume state variable - self._V.set_i_func(self.v_i_func, function='grounded_capacitor_model_volume') - self._V.set_i_inputs(pd.Series({'p':self._P_i.name})) - + # Set the initialization function for the input volume state variable + self._V.set_i_func(self.v_i_func, function='grounded_capacitor_model_volume') + self._V.set_i_inputs(pd.Series({'p':self._P_i.name})) + def __del__(self): super().__del__() - \ No newline at end of file diff --git a/src/ModularCirc/Components/Rlc_component.py b/src/ModularCirc/Components/Rlc_component.py index 6e84469..e980592 100644 --- a/src/ModularCirc/Components/Rlc_component.py +++ b/src/ModularCirc/Components/Rlc_component.py @@ -11,33 +11,33 @@ def q_o_dudt_func(t, y): return q_o_dudt_func class Rlc_component(Rc_component): - def __init__(self, - name:str, + def __init__(self, + name:str, time_object:TimeClass, - r:float, - c:float, + r:float, + c:float, l:float, v_ref:float, v:float=None, - p:float=None, + p:float=None, ) -> None: super().__init__(time_object=time_object, name=name, v=v, p=p, r=r, c=c, v_ref=v_ref) self.L = l - + def q_o_dudt_func(self, t, y): return resistor_impedance_flux_rate(t, y=y, r=self.R, l=self.L) - + def setup(self) -> None: Rc_component.setup(self) if (np.abs(self.L) > 1e-11): L = self.L R = self.R - self._Q_o.set_dudt_func(gen_q_o_dudt_func(r=R, l=L), + self._Q_o.set_dudt_func(gen_q_o_dudt_func(r=R, l=L), function_name='resistor_impedance_flux_rate') - self._Q_o.set_inputs(pd.Series({'p_in':self._P_i.name, + self._Q_o.set_inputs(pd.Series({'p_in':self._P_i.name, 'p_out':self._P_o.name, 'q_out':self._Q_o.name})) self._Q_o.set_u_func(None,None) - + def __del__(self): - super().__del__() \ No newline at end of file + super().__del__() diff --git a/src/ModularCirc/Components/Valve_maynard.py b/src/ModularCirc/Components/Valve_maynard.py index 7463e3e..7926794 100644 --- a/src/ModularCirc/Components/Valve_maynard.py +++ b/src/ModularCirc/Components/Valve_maynard.py @@ -1,6 +1,6 @@ from .ComponentBase import ComponentBase from ..Time import TimeClass -from ..HelperRoutines import maynard_valve_flow, maynard_phi_law, maynard_impedance_dqdt +from ..HelperRoutines import maynard_valve_flow, maynard_phi_law, maynard_impedance_dqdt from ..StateVariable import StateVariable import pandas as pd @@ -21,7 +21,7 @@ def phi_dudt_func(t, y): return phi_dudt_func class Valve_maynard(ComponentBase): - def __init__(self, + def __init__(self, name:str, time_object: TimeClass, Kc:float, @@ -29,12 +29,12 @@ def __init__(self, CQ:float, R :float=0.0, L :float=0.0, - RRA:float=0.0, + RRA:float=0.0, *args, **kwargs ) -> None: super().__init__(name=name, time_object=time_object) # allow for pressure gradient but not for flow - self.make_unique_io_state_variable(q_flag=True, p_flag=False) + self.make_unique_io_state_variable(q_flag=True, p_flag=False) # setting the bernoulli resistance value self.CQ = CQ # setting the resistance value @@ -47,27 +47,27 @@ def __init__(self, self.Kc, self.Ko = Kc, Ko # defining the valve opening factor state variable self._PHI = StateVariable(name=name+'_PHI', timeobj=time_object) - + @property def PHI(self): return self._PHI._u - + def q_i_u_func(self, t, y): return maynard_valve_flow(t, y=y, CQ=self.CQ, RRA=self.RRA) - + def q_i_dudt_func(self, t, y): return maynard_impedance_dqdt(t, y=y, CQ=self.CQ, RRA=self.RRA, L=self.L, R=self.R) - + def phi_dudt_func(self, t, y): return maynard_phi_law(t, y=y, Ko=self.Ko, Kc=self.Kc) - + def setup(self) -> None: CQ = self.CQ RRA = self.RRA if self.L < 1.0e-6: q_i_u_func = gen_q_i_u_func(CQ=CQ, RRA=RRA) self._Q_i.set_u_func(q_i_u_func, function_name='maynard_valve_flow') - self._Q_i.set_inputs(pd.Series({'p_in' : self._P_i.name, + self._Q_i.set_inputs(pd.Series({'p_in' : self._P_i.name, 'p_out': self._P_o.name, 'phi' : self._PHI.name})) else: @@ -75,15 +75,15 @@ def setup(self) -> None: L = self.L q_i_dudt_func = gen_q_i_dudt_func(CQ=CQ, RRA=RRA, L=L, R=R) self._Q_i.set_dudt_func(q_i_dudt_func, function_name='maynard_impedance_dqdt') - self._Q_i.set_inputs(pd.Series({'p_in' : self._P_i.name, + self._Q_i.set_inputs(pd.Series({'p_in' : self._P_i.name, 'p_out': self._P_o.name, 'q_in' : self._Q_i.name, 'phi' : self._PHI.name})) - + Ko = self.Ko Kc = self.Kc phi_dudt_func = gen_phi_dudt_func(Ko=Ko, Kc=Kc) self._PHI.set_dudt_func(phi_dudt_func, function_name='maynard_phi_law') self._PHI.set_inputs(pd.Series({'p_in' : self._P_i.name, 'p_out': self._P_o.name, - 'phi' : self._PHI.name})) \ No newline at end of file + 'phi' : self._PHI.name})) diff --git a/src/ModularCirc/Components/Valve_non_ideal.py b/src/ModularCirc/Components/Valve_non_ideal.py index baae67c..4581c87 100644 --- a/src/ModularCirc/Components/Valve_non_ideal.py +++ b/src/ModularCirc/Components/Valve_non_ideal.py @@ -1,37 +1,36 @@ from .ComponentBase import ComponentBase from ..Time import TimeClass -from ..HelperRoutines import non_ideal_diode_flow +from ..HelperRoutines import non_ideal_diode_flow import pandas as pd def gen_q_i_u_func(r, max_func): - def q_i_u_func(t, y): + def q_i_u_func(t, y): return non_ideal_diode_flow(t, y=y, r=r, max_func=max_func) return q_i_u_func class Valve_non_ideal(ComponentBase): - def __init__(self, + def __init__(self, name:str, time_object: TimeClass, - r:float, + r:float, max_func ) -> None: super().__init__(name=name, time_object=time_object) # allow for pressure gradient but not for flow - self.make_unique_io_state_variable(q_flag=True, p_flag=False) + self.make_unique_io_state_variable(q_flag=True, p_flag=False) # setting the resistance value self.R = r self.max_func = max_func - + def q_i_u_func(self, t, y): return non_ideal_diode_flow(t, y=y, r=self.R, max_func=self.max_func) - + def setup(self) -> None: r = self.R max_func = self.max_func # q_i_u_func = lambda t, y: non_ideal_diode_flow(t, y=y, r=r, max_func=max_func) q_i_u_func = gen_q_i_u_func(r=r, max_func=max_func) self._Q_i.set_u_func(q_i_u_func, function_name='non_ideal_diode_flow + max_func') - self._Q_i.set_inputs(pd.Series({'p_in':self._P_i.name, + self._Q_i.set_inputs(pd.Series({'p_in':self._P_i.name, 'p_out':self._P_o.name})) - \ No newline at end of file diff --git a/src/ModularCirc/Components/Valve_simple_bernoulli.py b/src/ModularCirc/Components/Valve_simple_bernoulli.py index 2b01c97..68fbf9c 100644 --- a/src/ModularCirc/Components/Valve_simple_bernoulli.py +++ b/src/ModularCirc/Components/Valve_simple_bernoulli.py @@ -1,32 +1,31 @@ from .ComponentBase import ComponentBase from ..Time import TimeClass -from ..HelperRoutines import simple_bernoulli_diode_flow +from ..HelperRoutines import simple_bernoulli_diode_flow import pandas as pd class Valve_simple_bernoulli(ComponentBase): - def __init__(self, + def __init__(self, name:str, time_object: TimeClass, CQ:float, - RRA:float=0.0, + RRA:float=0.0, ) -> None: super().__init__(name=name, time_object=time_object) # allow for pressure gradient but not for flow - self.make_unique_io_state_variable(q_flag=True, p_flag=False) + self.make_unique_io_state_variable(q_flag=True, p_flag=False) # setting the resistance value self.CQ = CQ # setting the relative regurgitant area self.RRA = RRA - + def q_i_u_func(self, t, y): return simple_bernoulli_diode_flow(t, y=y, CQ=self.CQ, RRA=self.RRA) - + def setup(self) -> None: CQ = self.CQ RRA = self.RRA q_i_u_func = lambda t, y: simple_bernoulli_diode_flow(t, y=y, CQ=CQ, RRA=RRA) self._Q_i.set_u_func(q_i_u_func, function_name='simple_bernoulli_diode_flow') - self._Q_i.set_inputs(pd.Series({'p_in':self._P_i.name, + self._Q_i.set_inputs(pd.Series({'p_in':self._P_i.name, 'p_out':self._P_o.name})) - \ No newline at end of file diff --git a/src/ModularCirc/Components/__init__.py b/src/ModularCirc/Components/__init__.py index 821dddc..c79b8d9 100644 --- a/src/ModularCirc/Components/__init__.py +++ b/src/ModularCirc/Components/__init__.py @@ -6,4 +6,4 @@ from .Rlc_component import Rlc_component from .Valve_non_ideal import Valve_non_ideal from .Valve_simple_bernoulli import Valve_simple_bernoulli -from .Valve_maynard import Valve_maynard \ No newline at end of file +from .Valve_maynard import Valve_maynard diff --git a/src/ModularCirc/HelperRoutines.py b/src/ModularCirc/HelperRoutines.py index a778f5a..294a482 100644 --- a/src/ModularCirc/HelperRoutines.py +++ b/src/ModularCirc/HelperRoutines.py @@ -1,4 +1,4 @@ -import numpy as np +import numpy as np from .Time import TimeClass import numba as nb @@ -6,7 +6,7 @@ # from numba import jit from collections.abc import Callable -def resistor_model_flow(t:float, +def resistor_model_flow(t:float, p_in:float=None, p_out:float=None, r:float=None, @@ -41,7 +41,7 @@ def resistor_model_dp(q_in:float, r:float) -> float: return q_in * r # @nb.njit(cache=True) -def resistor_impedance_flux_rate(t:float, +def resistor_impedance_flux_rate(t:float, p_in:float=None, p_out:float=None, q_out:float=None, @@ -49,7 +49,7 @@ def resistor_impedance_flux_rate(t:float, l:float=None, y:np.ndarray[float]=None) -> float: """ - Resistor and impedance in series flux rate of change model. + Resistor and impedance in series flux rate of change model. Args: t (float): current time @@ -64,16 +64,16 @@ def resistor_impedance_flux_rate(t:float, """ if y is not None: p_in, p_out, q_out = y[:3] - return (p_in - p_out - q_out * r ) / l + return (p_in - p_out - q_out * r ) / l -def grounded_capacitor_model_pressure(t:float, - v:float=None, +def grounded_capacitor_model_pressure(t:float, + v:float=None, v_ref:float=None, - c:float=None, + c:float=None, y:np.ndarray[float]=None ) -> float: """ - Capacitor model with constant capacitance. + Capacitor model with constant capacitance. Args: ---- @@ -100,10 +100,10 @@ def grounded_capacitor_model_volume(t:float, return v_ref + p * c # @nb.njit(cache=True) -def grounded_capacitor_model_dpdt(t:float, - q_in:float=None, - q_out:float=None, - c:float=None, +def grounded_capacitor_model_dpdt(t:float, + q_in:float=None, + q_out:float=None, + c:float=None, y:np.ndarray[float]=None ) -> float: if y is not None: @@ -111,8 +111,8 @@ def grounded_capacitor_model_dpdt(t:float, return (q_in - q_out) / c # @nb.njit(cache=True) -def chamber_volume_rate_change(t:float, - q_in:float=None, +def chamber_volume_rate_change(t:float, + q_in:float=None, q_out:float=None, y:np.ndarray[float]=None ) -> float: @@ -131,7 +131,7 @@ def chamber_volume_rate_change(t:float, return q_in - q_out # @nb.njit(cache=True) -def relu_max(val:float) -> float: +def relu_max(val:float) -> float: return np.maximum(val, 0.0) def softplus(val:float, alpha:float=0.2) -> float: @@ -141,11 +141,11 @@ def softplus(val:float, alpha:float=0.2) -> float: y = val.copy() y[alpha * y <= 20.0] = 1/ alpha * np.log(1 + np.exp(alpha * y[alpha * y <=20.0])) return y - + def get_softplus_max(alpha:float): """ Method for generating softmax lambda function based on predefined alpha values - + Args: ---- alpha (float): softplus alpha value @@ -156,10 +156,10 @@ def get_softplus_max(alpha:float): """ return lambda val : softplus(val=val, alpha=alpha) -def non_ideal_diode_flow(t:float, - p_in:float=None, - p_out:float=None, - r:float=None, +def non_ideal_diode_flow(t:float, + p_in:float=None, + p_out:float=None, + r:float=None, max_func:Callable[[float],float]=relu_max, y:np.ndarray[float]=None, ) -> float: @@ -181,11 +181,11 @@ def non_ideal_diode_flow(t:float, return (max_func((p_in - p_out)/ r)) # @jit(cache=True, nopython=True) -def simple_bernoulli_diode_flow(t:float, - p_in:float=None, - p_out:float=None, +def simple_bernoulli_diode_flow(t:float, + p_in:float=None, + p_out:float=None, CQ:float=None, - RRA:float=0.0, + RRA:float=0.0, y:np.ndarray[float]=None, ) -> float: """ @@ -204,8 +204,8 @@ def simple_bernoulli_diode_flow(t:float, if y is not None: p_in, p_out = y[:2] dp = p_in - p_out - return np.where(dp >= 0.0, - CQ * np.sqrt(np.abs(dp)), + return np.where(dp >= 0.0, + CQ * np.sqrt(np.abs(dp)), -CQ * RRA *np.sqrt(np.abs(dp))) # @jit(cache=True, nopython=True) @@ -260,7 +260,7 @@ def leaky_diode_flow(p_in:float, p_out:float, r_o:float, r_r:float) -> float: Leaky diode model that outputs the flow rate through a leaky diode Args: - p_in (float): input pressure + p_in (float): input pressure p_out (float): output pressure r_o (float): outflow resistance r_r (float): regurgitant flow resistance @@ -312,7 +312,7 @@ def activation_function_2(t:float, tr:float, td:float, dt: bool=True) -> float: 0.0 ) return result - + def activation_function_3(t:float, tpwb:float, tpww:float, dt: bool=True) -> float: if not dt: result = ( @@ -363,7 +363,7 @@ def chamber_linear_elastic_law(v:float, E:float, v_ref:float, *args, **kwargs) - v (float): volume E (float): Elastance v_ref (float): reference volume - + Returns: float: chamber pressure """ @@ -385,8 +385,8 @@ def chamber_exponential_law(v:float, E:float, k:float, v_ref:float, *args, **kwa return E * np.exp(k * (v - v_ref) - 1) def chamber_pressure_function(t:float, v:float, v_ref:float, E_pas:float, E_act:float, - activation_function = activation_function_1, - active_law = chamber_linear_elastic_law, + activation_function = activation_function_1, + active_law = chamber_linear_elastic_law, passive_law = chamber_linear_elastic_law, *args, **kwargs) ->float: """ @@ -406,9 +406,9 @@ def chamber_pressure_function(t:float, v:float, v_ref:float, E_pas:float, E_act: float: pressure """ a = activation_function(t) - return (a * active_law(v=v, v_ref=v_ref,t=t, E=E_act, **kwargs) + return (a * active_law(v=v, v_ref=v_ref,t=t, E=E_act, **kwargs) + (1 - a) * passive_law(v=v, v_ref=v_ref, t=t, E=E_pas, **kwargs)) - + def time_shift(t:float, shift:float=np.nan, tcycle:float=0.0): if shift is np.nan: return t @@ -423,4 +423,4 @@ def time_shift(t:float, shift:float=np.nan, tcycle:float=0.0): END = '\033[0m' def bold_text(str_:str): - return BOLD + YELLOW + str_ + END \ No newline at end of file + return BOLD + YELLOW + str_ + END diff --git a/src/ModularCirc/Models/KPat5MixedModel.py b/src/ModularCirc/Models/KPat5MixedModel.py index 30aa2e1..feb9644 100644 --- a/src/ModularCirc/Models/KPat5MixedModel.py +++ b/src/ModularCirc/Models/KPat5MixedModel.py @@ -28,12 +28,12 @@ class KPat5MixedModel(OdeModel): def __init__(self, time_setup_dict, parobj:po=k2006, suppress_printing:bool=False) -> None: super().__init__(time_setup_dict) self.name = 'KorakianitisModel' - + if not suppress_printing: print(parobj) - + # The components... for key, name in zip(parobj.components.keys(), FULL_NAMES): - if key in parobj._vessels: + if key in parobj._vessels: class_ = Rlc_component elif key in parobj._valves: class_ = Valve_simple_bernoulli @@ -42,14 +42,14 @@ def __init__(self, time_setup_dict, parobj:po=k2006, suppress_printing:bool=Fals else: raise Exception(f'Component name {key} not in the model list.') self.components[key] = class_(name=name, - time_object=self.time_object, + time_object=self.time_object, **parobj[key].to_dict()) - if key not in parobj._valves: + if key not in parobj._valves: self.set_v_sv(key) # else: # self.set_phi_sv(key) self.components[key].setup() - + self.connect_modules(self.components['lv'], self.components['ao'], plabel='p_lv', @@ -125,12 +125,12 @@ def __init__(self, time_setup_dict, parobj:po=k2006, suppress_printing:bool=Fals self.components['lv'], plabel='p_lv', qlabel='q_mi') - + for component in self.components.values(): component.setup() - + # def set_phi_sv(self, comp_key:str) -> None: # phi_key = 'phi_' + comp_key # self._state_variable_dict[phi_key] = self.components[comp_key]._PHI # self._state_variable_dict[phi_key].set_name(phi_key) - # self.all_sv_data[phi_key] = self.components[comp_key].PHI \ No newline at end of file + # self.all_sv_data[phi_key] = self.components[comp_key].PHI diff --git a/src/ModularCirc/Models/KPat5MixedModel_parameters.py b/src/ModularCirc/Models/KPat5MixedModel_parameters.py index f72ad75..4a47578 100644 --- a/src/ModularCirc/Models/KPat5MixedModel_parameters.py +++ b/src/ModularCirc/Models/KPat5MixedModel_parameters.py @@ -17,12 +17,12 @@ 'po', # pulmonary valve 'pas', # pulmonary artery sinus ############################### - 'pat0', # pulmonary artery - 'pat1', # pulmonary artery - 'pat2', # pulmonary artery - 'pat3', # pulmonary artery + 'pat0', # pulmonary artery + 'pat1', # pulmonary artery + 'pat2', # pulmonary artery + 'pat3', # pulmonary artery 'pat4', # pulmonary artery - ############################### + ############################### 'pvn' # pulmonary vein ] VESSELS = ['sas', 'sat', 'svn', 'pas', 'pat0', 'pat1', 'pat2', 'pat3', 'pat4', 'pvn'] @@ -47,20 +47,20 @@ def __init__(self, name='Korakianitis 2006') -> None: for type_, type_var in [[VESSELS, VESSELS_PAR], [VALVES, VALVES_PAR], [CHAMBERS, CHAMBERS_PAR]]: for key in type_: self[key] = pd.Series(index=type_var, dtype=object) - + self._vessels = VESSELS self._valves = VALVES self._chambers= CHAMBERS - + self.set_chamber_comp('lv', E_pas= 0.1, k_pas=0.01, E_act= 2.5, v_ref=50.0, tr = 0.30, td = 0.450, v=100.) self.set_chamber_comp('la', E_pas= 0.15, k_pas=0.01, E_act= 0.25, v_ref=10.0, tpwb = 0.0, tpww = 0.09, delay=0.08, v=0.0) self.set_chamber_comp('rv', E_pas= 0.1, k_pas=0.01, E_act= 1.15, v_ref=50., tr=0.30, td=0.45, v=100.) self.set_chamber_comp('ra', E_pas= 0.15, k_pas=0.01, E_act= 0.25, v_ref=10., tpwb=0.0, tpww=0.09, delay=0.08, v=0.0) - + self.set_activation_function('lv', af=activation_function_2) self.set_activation_function('rv', af=activation_function_2) - + self.set_activation_function('la', af=activation_function_3) self.set_activation_function('ra', af=activation_function_3) @@ -69,7 +69,7 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_rlc_comp('sas', r=0.003, c=0.08, l=0.000062, v=0.0, v_ref=0.0) self.set_rlc_comp('sat', r=(0.05 + 0.5 + 0.52), c=1.6 , l=0.0017 , v=550.0, v_ref=0.0) self.set_rlc_comp('svn', r=0.075, c=20.5, v=0.0, v_ref=0.0) - + # pulmonary circulation self.set_rlc_comp('pas', r=0.002 , c=0.18, l=0.000052, v=0.0, v_ref=0.0) ############################################################################################# @@ -80,25 +80,25 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_rlc_comp('pat4', r=(0.01/5.+0.05+0.25), c=(3.8/5.) , l=(0.0017/5.), v=0.0, v_ref=0.0) ############################################################################################# self.set_rlc_comp('pvn', r=0.006 , c=20.5 , v=0.0, v_ref=0.0) - + # valves self.set_valve_comp('ao', CQ=350., RRA=0.0) self.set_valve_comp('mi', CQ=400., RRA=0.0) self.set_valve_comp('po', CQ=350., RRA=0.0) self.set_valve_comp('ti', CQ=400., RRA=0.0) - + def set_chamber_comp(self, key, **kwargs): self._set_comp(key=key, set=CHAMBERS, **kwargs) - + def set_activation_function(self, key, af): self._set_comp(key, set=CHAMBERS, af=af) - + def set_rlc_comp(self, key, **kwargs): self._set_comp(key=key, set=VESSELS, **kwargs) - + def set_valve_comp(self, key, **kwargs): self._set_comp(key=key, set=VALVES, **kwargs) - + def set(self, key, **kwargs): if key in CHAMBERS: self._set_comp(key=key, set=CHAMBERS, **kwargs) diff --git a/src/ModularCirc/Models/KorakianitisMaynardModel.py b/src/ModularCirc/Models/KorakianitisMaynardModel.py index 7874e46..b923493 100644 --- a/src/ModularCirc/Models/KorakianitisMaynardModel.py +++ b/src/ModularCirc/Models/KorakianitisMaynardModel.py @@ -24,12 +24,12 @@ class KorakianitisMaynardModel(OdeModel): def __init__(self, time_setup_dict, parobj:po=k2006) -> None: super().__init__(time_setup_dict) self.name = 'KorakianitisModel' - + print(parobj) - + # The components... for key, name in zip(parobj.components.keys(), FULL_NAMES): - if key in parobj._vessels: + if key in parobj._vessels: class_ = Rlc_component elif key in parobj._valves: class_ = Valve_maynard @@ -38,14 +38,14 @@ def __init__(self, time_setup_dict, parobj:po=k2006) -> None: else: raise Exception(f'Component name {key} not in the model list.') self.components[key] = class_(name=name, - time_object=self.time_object, + time_object=self.time_object, **parobj[key].to_dict()) - if key not in parobj._valves: + if key not in parobj._valves: self.set_v_sv(key) else: self.set_phi_sv(key) self.components[key].setup() - + self.connect_modules(self.components['lv'], self.components['ao'], plabel='p_lv', @@ -102,13 +102,13 @@ def __init__(self, time_setup_dict, parobj:po=k2006) -> None: self.components['lv'], plabel='p_lv', qlabel='q_mi') - + for component in self.components.values(): component.setup() - + def set_phi_sv(self, comp_key:str) -> None: phi_key = 'phi_' + comp_key self._state_variable_dict[phi_key] = self.components[comp_key]._PHI self._state_variable_dict[phi_key].set_name(phi_key) self.all_sv_data[phi_key] = self.components[comp_key].PHI - self.components[comp_key]._PHI._u = self.all_sv_data[phi_key] \ No newline at end of file + self.components[comp_key]._PHI._u = self.all_sv_data[phi_key] diff --git a/src/ModularCirc/Models/KorakianitisMaynardModel_parameters.py b/src/ModularCirc/Models/KorakianitisMaynardModel_parameters.py index a55bb8a..6c9a840 100644 --- a/src/ModularCirc/Models/KorakianitisMaynardModel_parameters.py +++ b/src/ModularCirc/Models/KorakianitisMaynardModel_parameters.py @@ -16,7 +16,7 @@ 'rv', # right ventricle 'po', # pulmonary valve 'pas', # pulmonary artery sinus - 'pat', # pulmonary artery + 'pat', # pulmonary artery 'pvn' # pulmonary vein ] VESSELS = ['sas', 'sat', 'svn', 'pas', 'pat', 'pvn'] @@ -41,11 +41,11 @@ def __init__(self, name='Korakianitis 2006') -> None: for type_, type_var in [[VESSELS, VESSELS_PAR], [VALVES, VALVES_PAR], [CHAMBERS, CHAMBERS_PAR]]: for key in type_: self[key] = pd.Series(index=type_var, dtype=object) - + self._vessels = VESSELS self._valves = VALVES self._chambers= CHAMBERS - + self.set_chamber_comp('lv', E_pas= 0.1, E_act= 2.5, v_ref=5.0, k_pas=0.1, tr = 0.30, td = 0.450, v=500.) self.set_chamber_comp('la', E_pas= 0.15, E_act= 0.25, v_ref=4.0, k_pas=0.1, tpwb = 0.0, tpww = 0.09, delay=0.08, v=0.0) self.set_chamber_comp('rv', E_pas= 0.1, E_act= 1.15, v_ref=10., k_pas=0.1, tr=0.30, td=0.45, v=400.) @@ -53,7 +53,7 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_activation_function('lv', af=activation_function_2) self.set_activation_function('rv', af=activation_function_2) - + self.set_activation_function('la', af=activation_function_3) self.set_activation_function('ra', af=activation_function_3) @@ -62,27 +62,27 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_rlc_comp('sas', r=0.003, c=0.08, l=0.000062, v=0.0, v_ref=0.0) self.set_rlc_comp('sat', r=(0.05 + 0.5 + 0.52), c=1.6 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('svn', r=0.075, c=20.5, v=0.0, v_ref=0.0) - + # pulmonary circulation self.set_rlc_comp('pas', r=0.002 , c=0.18, l=0.000052, v=0.0, v_ref=0.0) self.set_rlc_comp('pat', r=(0.01+0.05+0.25), c=3.8 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('pvn', r=0.006 , c=20.5 , v=0.0, v_ref=0.0) - + # valves - dyn = 1333.22 + dyn = 1333.22 self.set_valve_comp('ao', CQ=350., RRA=0.0, Ko = 0.012/dyn, Kc = 0.012/dyn, L=0.0, R=0.0) self.set_valve_comp('mi', CQ=400., RRA=0.0, Ko = 0.03/dyn, Kc = 0.04/dyn, L=0.0, R=0.0 ) self.set_valve_comp('po', CQ=350., RRA=0.0, Ko = 0.02/dyn, Kc = 0.02/dyn, L=0.0, R=0.0 ) self.set_valve_comp('ti', CQ=400., RRA=0.0, Ko = 0.03/dyn, Kc = 0.04/dyn, L=0.0, R=0.0 ) - + def set_chamber_comp(self, key, **kwargs): self._set_comp(key=key, set=CHAMBERS, **kwargs) - + def set_activation_function(self, key, af): self._set_comp(key, set=CHAMBERS, af=af) - + def set_rlc_comp(self, key, **kwargs): self._set_comp(key=key, set=VESSELS, **kwargs) - + def set_valve_comp(self, key, **kwargs): self._set_comp(key=key, set=VALVES, **kwargs) diff --git a/src/ModularCirc/Models/KorakianitisMixedMaynardModel.py b/src/ModularCirc/Models/KorakianitisMixedMaynardModel.py index c8ab554..2e23339 100644 --- a/src/ModularCirc/Models/KorakianitisMixedMaynardModel.py +++ b/src/ModularCirc/Models/KorakianitisMixedMaynardModel.py @@ -24,12 +24,12 @@ class KorakianitisMixedMaynardModel(OdeModel): def __init__(self, time_setup_dict, parobj:po=k2006) -> None: super().__init__(time_setup_dict) self.name = 'KorakianitisModel' - + print(parobj) - + # The components... for key, name in zip(parobj.components.keys(), FULL_NAMES): - if key in parobj._vessels: + if key in parobj._vessels: class_ = Rlc_component elif key in parobj._valves: class_ = Valve_maynard @@ -38,14 +38,14 @@ def __init__(self, time_setup_dict, parobj:po=k2006) -> None: else: raise Exception(f'Component name {key} not in the model list.') self.components[key] = class_(name=name, - time_object=self.time_object, + time_object=self.time_object, **parobj[key].to_dict()) - if key not in parobj._valves: + if key not in parobj._valves: self.set_v_sv(key) else: self.set_phi_sv(key) self.components[key].setup() - + self.connect_modules(self.components['lv'], self.components['ao'], plabel='p_lv', @@ -102,12 +102,12 @@ def __init__(self, time_setup_dict, parobj:po=k2006) -> None: self.components['lv'], plabel='p_lv', qlabel='q_mi') - + for component in self.components.values(): component.setup() - + def set_phi_sv(self, comp_key:str) -> None: phi_key = 'phi_' + comp_key self._state_variable_dict[phi_key] = self.components[comp_key]._PHI self._state_variable_dict[phi_key].set_name(phi_key) - self.all_sv_data[phi_key] = self.components[comp_key].PHI \ No newline at end of file + self.all_sv_data[phi_key] = self.components[comp_key].PHI diff --git a/src/ModularCirc/Models/KorakianitisMixedMaynardModel_parameters.py b/src/ModularCirc/Models/KorakianitisMixedMaynardModel_parameters.py index 6dc1d8f..1dd6299 100644 --- a/src/ModularCirc/Models/KorakianitisMixedMaynardModel_parameters.py +++ b/src/ModularCirc/Models/KorakianitisMixedMaynardModel_parameters.py @@ -16,7 +16,7 @@ 'rv', # right ventricle 'po', # pulmonary valve 'pas', # pulmonary artery sinus - 'pat', # pulmonary artery + 'pat', # pulmonary artery 'pvn' # pulmonary vein ] VESSELS = ['sas', 'sat', 'svn', 'pas', 'pat', 'pvn'] @@ -41,11 +41,11 @@ def __init__(self, name='Korakianitis 2006') -> None: for type_, type_var in [[VESSELS, VESSELS_PAR], [VALVES, VALVES_PAR], [CHAMBERS, CHAMBERS_PAR]]: for key in type_: self[key] = pd.Series(index=type_var, dtype=object) - + self._vessels = VESSELS self._valves = VALVES self._chambers= CHAMBERS - + self.set_chamber_comp('lv', E_pas= 0.1, E_act= 2.5, v_ref=5.0, k_pas=0.1, tr = 0.30, td = 0.450, v=50.) self.set_chamber_comp('la', E_pas= 0.15, E_act= 0.25, v_ref=4.0, k_pas=0.1, tpwb = 0.0, tpww = 0.09, delay=0.08, v=0.0) self.set_chamber_comp('rv', E_pas= 0.1, E_act= 1.15, v_ref=10., k_pas=0.1, tr=0.30, td=0.45, v=100.) @@ -53,7 +53,7 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_activation_function('lv', af=activation_function_2) self.set_activation_function('rv', af=activation_function_2) - + self.set_activation_function('la', af=activation_function_3) self.set_activation_function('ra', af=activation_function_3) @@ -62,27 +62,27 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_rlc_comp('sas', r=0.003, c=0.08, l=0.000062, v=450.0, v_ref=0.0) self.set_rlc_comp('sat', r=(0.05 + 0.5 + 0.52), c=1.6 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('svn', r=0.075, c=20.5, v=0.0, v_ref=0.0) - + # pulmonary circulation self.set_rlc_comp('pas', r=0.002 , c=0.18, l=0.000052, v=200.0, v_ref=0.0) self.set_rlc_comp('pat', r=(0.01+0.05+0.25), c=3.8 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('pvn', r=0.006 , c=20.5 , v=0.0, v_ref=0.0) - + # valves - dyn = 1333.22 + dyn = 1333.22 self.set_valve_comp('ao', CQ=350., RRA=0.0, Ko = 0.012/dyn, Kc = 0.012/dyn, L=0.0, R=0.0) self.set_valve_comp('mi', CQ=400., RRA=0.0, Ko = 0.03/dyn, Kc = 0.04/dyn, L=0.0, R=0.0 ) self.set_valve_comp('po', CQ=350., RRA=0.0, Ko = 0.02/dyn, Kc = 0.02/dyn, L=0.0, R=0.0 ) self.set_valve_comp('ti', CQ=400., RRA=0.0, Ko = 0.03/dyn, Kc = 0.04/dyn, L=0.0, R=0.0 ) - + def set_chamber_comp(self, key, **kwargs): self._set_comp(key=key, set=CHAMBERS, **kwargs) - + def set_activation_function(self, key, af): self._set_comp(key, set=CHAMBERS, af=af) - + def set_rlc_comp(self, key, **kwargs): self._set_comp(key=key, set=VESSELS, **kwargs) - + def set_valve_comp(self, key, **kwargs): self._set_comp(key=key, set=VALVES, **kwargs) diff --git a/src/ModularCirc/Models/KorakianitisMixedModel.py b/src/ModularCirc/Models/KorakianitisMixedModel.py index 525e5d5..cd174f6 100644 --- a/src/ModularCirc/Models/KorakianitisMixedModel.py +++ b/src/ModularCirc/Models/KorakianitisMixedModel.py @@ -25,12 +25,12 @@ class KorakianitisMixedModel(OdeModel): def __init__(self, time_setup_dict, parobj:po=KorakianitisMixedModel_parameters, suppress_printing:bool=False) -> None: super().__init__(time_setup_dict) self.name = 'KorakianitisModel' - + if not suppress_printing: print(parobj) - + # The components... for key, name in zip(parobj.components.keys(), FULL_NAMES): - if key in parobj._vessels: + if key in parobj._vessels: class_ = Rlc_component elif key in parobj._valves: class_ = Valve_simple_bernoulli @@ -39,14 +39,14 @@ def __init__(self, time_setup_dict, parobj:po=KorakianitisMixedModel_parameters, else: raise Exception(f'Component name {key} not in the model list.') self.components[key] = class_(name=name, - time_object=self.time_object, + time_object=self.time_object, **parobj[key].to_dict()) - if key not in parobj._valves: + if key not in parobj._valves: self.set_v_sv(key) # else: # self.set_phi_sv(key) self.components[key].setup() - + self.connect_modules(self.components['lv'], self.components['ao'], plabel='p_lv', @@ -103,12 +103,12 @@ def __init__(self, time_setup_dict, parobj:po=KorakianitisMixedModel_parameters, self.components['lv'], plabel='p_lv', qlabel='q_mi') - + for component in self.components.values(): component.setup() - + # def set_phi_sv(self, comp_key:str) -> None: # phi_key = 'phi_' + comp_key # self._state_variable_dict[phi_key] = self.components[comp_key]._PHI # self._state_variable_dict[phi_key].set_name(phi_key) - # self.all_sv_data[phi_key] = self.components[comp_key].PHI \ No newline at end of file + # self.all_sv_data[phi_key] = self.components[comp_key].PHI diff --git a/src/ModularCirc/Models/KorakianitisMixedModel_parameters.py b/src/ModularCirc/Models/KorakianitisMixedModel_parameters.py index 26e26c9..deffc4a 100644 --- a/src/ModularCirc/Models/KorakianitisMixedModel_parameters.py +++ b/src/ModularCirc/Models/KorakianitisMixedModel_parameters.py @@ -16,7 +16,7 @@ 'rv', # right ventricle 'po', # pulmonary valve 'pas', # pulmonary artery sinus - 'pat', # pulmonary artery + 'pat', # pulmonary artery 'pvn' # pulmonary vein ] VESSELS = ['sas', 'sat', 'svn', 'pas', 'pat', 'pvn'] @@ -49,11 +49,11 @@ def __init__(self, name='Korakianitis 2006') -> None: for type_, type_var in [[VESSELS, VESSELS_PAR], [VALVES, VALVES_PAR], [CHAMBERS, CHAMBERS_PAR]]: for key in type_: self[key] = pd.Series(index=type_var, dtype=object) - + self._vessels = VESSELS self._valves = VALVES self._chambers= CHAMBERS - + self.set_chamber_comp('lv', E_pas= 1.7, E_act= 2.5, v_ref=5.0, k_pas=0.01, tr = 0.30, td = 0.450, v=50.) self.set_chamber_comp('la', E_pas= 0.5, E_act= 0.25, v_ref=4.0, k_pas=0.01, tpwb = 0.0, tpww = 0.09, delay=0.08, v=0.0) self.set_chamber_comp('rv', E_pas= 0.67, E_act= 1.15, v_ref=10., k_pas=0.01, tr=0.30, td=0.45, v=100.) @@ -61,7 +61,7 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_activation_function('lv', af=activation_function_2) self.set_activation_function('rv', af=activation_function_2) - + self.set_activation_function('la', af=activation_function_3) self.set_activation_function('ra', af=activation_function_3) @@ -70,31 +70,31 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_rlc_comp('sas', r=0.003, c=0.08, l=0.000062, v=450.0, v_ref=0.0) self.set_rlc_comp('sat', r=(0.05 + 0.5 + 0.52), c=1.6 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('svn', r=0.075, c=20.5, v=0.0, v_ref=0.0) - + # pulmonary circulation self.set_rlc_comp('pas', r=0.002 , c=0.18, l=0.000052, v=200.0, v_ref=0.0) self.set_rlc_comp('pat', r=(0.01+0.05+0.25), c=3.8 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('pvn', r=0.006 , c=20.5 , v=0.0, v_ref=0.0) - + # valves - dyn = 1333.22 + dyn = 1333.22 self.set_valve_comp('ao', CQ=350., RRA=0.0) self.set_valve_comp('mi', CQ=400., RRA=0.0) self.set_valve_comp('po', CQ=350., RRA=0.0) self.set_valve_comp('ti', CQ=400., RRA=0.0) - + def set_chamber_comp(self, key, **kwargs): self._set_comp(key=key, set=CHAMBERS, **kwargs) - + def set_activation_function(self, key, af): self._set_comp(key, set=CHAMBERS, af=af) - + def set_rlc_comp(self, key, **kwargs): self._set_comp(key=key, set=VESSELS, **kwargs) - + def set_valve_comp(self, key, **kwargs): self._set_comp(key=key, set=VALVES, **kwargs) - + def set(self, key, **kwargs): if key in CHAMBERS: self._set_comp(key=key, set=CHAMBERS, **kwargs) diff --git a/src/ModularCirc/Models/KorakianitisModel.py b/src/ModularCirc/Models/KorakianitisModel.py index 1592487..91db341 100644 --- a/src/ModularCirc/Models/KorakianitisModel.py +++ b/src/ModularCirc/Models/KorakianitisModel.py @@ -24,12 +24,12 @@ class KorakianitisModel(OdeModel): def __init__(self, time_setup_dict, parobj:po=k2006, suppress_printing:bool=False) -> None: super().__init__(time_setup_dict) self.name = 'KorakianitisModel' - + if not suppress_printing: print(parobj) - + # The components... for key, name in zip(parobj.components.keys(), FULL_NAMES): - if key in parobj._vessels: + if key in parobj._vessels: class_ = Rlc_component elif key in parobj._valves: class_ = Valve_simple_bernoulli @@ -38,11 +38,11 @@ def __init__(self, time_setup_dict, parobj:po=k2006, suppress_printing:bool=Fals else: raise Exception(f'Component name {key} not in the model list.') self.components[key] = class_(name=name, - time_object=self.time_object, + time_object=self.time_object, **parobj[key].to_dict()) if key not in parobj._valves: self.set_v_sv(key) self.components[key].setup() - + self.connect_modules(self.components['lv'], self.components['ao'], plabel='p_lv', @@ -99,6 +99,6 @@ def __init__(self, time_setup_dict, parobj:po=k2006, suppress_printing:bool=Fals self.components['lv'], plabel='p_lv', qlabel='q_mi') - + for component in self.components.values(): - component.setup() \ No newline at end of file + component.setup() diff --git a/src/ModularCirc/Models/KorakianitisModel_parameters.py b/src/ModularCirc/Models/KorakianitisModel_parameters.py index e4414f5..d40013a 100644 --- a/src/ModularCirc/Models/KorakianitisModel_parameters.py +++ b/src/ModularCirc/Models/KorakianitisModel_parameters.py @@ -16,7 +16,7 @@ 'rv', # right ventricle 'po', # pulmonary valve 'pas', # pulmonary artery sinus - 'pat', # pulmonary artery + 'pat', # pulmonary artery 'pvn' # pulmonary vein ] VESSELS = ['sas', 'sat', 'svn', 'pas', 'pat', 'pvn'] @@ -42,13 +42,13 @@ def __init__(self, name='Korakianitis 2006') -> None: for type_, type_var in [[VESSELS, VESSELS_PAR], [VALVES, VALVES_PAR], [CHAMBERS, CHAMBERS_PAR]]: for key in type_: self[key] = pd.Series(index=type_var, dtype=object) - - # self.timings = {key : pd.Series(index=TIMINGS) for key in CHAMBERS} - + + # self.timings = {key : pd.Series(index=TIMINGS) for key in CHAMBERS} + self._vessels = VESSELS self._valves = VALVES self._chambers= CHAMBERS - + self.set_chamber_comp('lv', E_pas= 0.1, E_act= 2.5, v_ref=5.0, tr = 0.30, td = 0.450, v=500.) self.set_chamber_comp('la', E_pas= 0.15, E_act= 0.25, v_ref=4.0, tpwb = 0.0, tpww = 0.09, delay=0.08, v=0.0) self.set_chamber_comp('rv', E_pas= 0.1, E_act= 1.15, v_ref=10., tr=0.30, td=0.45, v=400.) @@ -56,7 +56,7 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_activation_function('lv', af=activation_function_2) self.set_activation_function('rv', af=activation_function_2) - + self.set_activation_function('la', af=activation_function_3) self.set_activation_function('ra', af=activation_function_3) @@ -65,30 +65,30 @@ def __init__(self, name='Korakianitis 2006') -> None: self.set_rlc_comp('sas', r=0.003, c=0.08, l=0.000062, v=0.0, v_ref=0.0) self.set_rlc_comp('sat', r=(0.05 + 0.5 + 0.52), c=1.6 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('svn', r=0.075, c=20.5, v=0.0, v_ref=0.0) - + # pulmonary circulation self.set_rlc_comp('pas', r=0.002 , c=0.18, l=0.000052, v=0.0, v_ref=0.0) self.set_rlc_comp('pat', r=(0.01+0.05+0.25), c=3.8 , l=0.0017 , v=0.0, v_ref=0.0) self.set_rlc_comp('pvn', r=0.006 , c=20.5 , v=0.0, v_ref=0.0) - + # valves self.set_valve_comp('ao', CQ=350., RRA=0.0) self.set_valve_comp('mi', CQ=400., RRA=0.0) self.set_valve_comp('po', CQ=350., RRA=0.0) self.set_valve_comp('ti', CQ=400., RRA=0.0) - + def set_chamber_comp(self, key, **kwargs): self._set_comp(key=key, set=CHAMBERS, **kwargs) - + def set_activation_function(self, key, af): self._set_comp(key, set=CHAMBERS, af=af) - + def set_rlc_comp(self, key, **kwargs): self._set_comp(key=key, set=VESSELS, **kwargs) - + def set_valve_comp(self, key, **kwargs): self._set_comp(key=key, set=VALVES, **kwargs) - + def set(self, key, **kwargs): if key in CHAMBERS: self._set_comp(key=key, set=CHAMBERS, **kwargs) diff --git a/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel.py b/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel.py index 49a1203..d1fdbf3 100644 --- a/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel.py +++ b/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel.py @@ -28,12 +28,12 @@ class MixedHeartMaynard4eWindkessel(OdeModel): def __init__(self, time_setup_dict, parobj:po=MHM4W_parobj) -> None: super().__init__(time_setup_dict) self.name = 'MixedHeartMaynard4eWindkessel' - + print(parobj) - + # The components... for key, name in zip(parobj.components.keys(), FULL_NAMES): - if key in parobj._vessels: + if key in parobj._vessels: class_ = Rlc_component elif key in parobj._imp or key in parobj._cap: class_ = R_component @@ -44,15 +44,15 @@ def __init__(self, time_setup_dict, parobj:po=MHM4W_parobj) -> None: else: raise Exception(f'Component name {key} not in the model list.') self.components[key] = class_(name=name, - time_object=self.time_object, + time_object=self.time_object, **parobj[key].to_dict()) - - if key not in parobj._valves + parobj._cap + parobj._imp: + + if key not in parobj._valves + parobj._cap + parobj._imp: self.set_v_sv(key) # else: # self.set_phi_sv(key) self.components[key].setup() - + self.connect_modules(self.components['lv'], self.components['ao'], plabel='p_lv', @@ -124,9 +124,9 @@ def __init__(self, time_setup_dict, parobj:po=MHM4W_parobj) -> None: qlabel='q_mi') for component in self.components.values(): component.setup() - + def set_phi_sv(self, comp_key:str) -> None: phi_key = 'phi_' + comp_key self._state_variable_dict[phi_key] = self.components[comp_key]._PHI self._state_variable_dict[phi_key].set_name(phi_key) - self.all_sv_data[phi_key] = self.components[comp_key].PHI \ No newline at end of file + self.all_sv_data[phi_key] = self.components[comp_key].PHI diff --git a/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel_parameters.py b/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel_parameters.py index 388363a..6a5c0e5 100644 --- a/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel_parameters.py +++ b/src/ModularCirc/Models/MixedHeartMaynard4eWindkessel_parameters.py @@ -17,7 +17,7 @@ 'rv', # right ventricle 'po', # pulmonary valve 'pai', # pulmonary artery impedance - 'pa', # pulmonary artery + 'pa', # pulmonary artery 'pc', # pulmonary capilary 'pv', # pulmonary vein ] @@ -57,13 +57,13 @@ def __init__(self, name='MixedHeartMaynard4eWindkessel_parameters') -> None: for type_, type_var in OBJ_PAR_PAIRS: for key in type_: self[key] = pd.Series(index=type_var, dtype=object) - + self._vessels = VESSELS self._valves = VALVES self._chambers= CHAMBERS self._imp = IMPEDANCES self._cap = CAPILARIES - + self.set_chamber_comp('lv', E_pas= 0.1, E_act= 2.5, k_pas=0.03, v_ref=5.0, tr = 0.30, td = 0.450, v=50.) self.set_chamber_comp('la', E_pas= 0.15, E_act= 0.25, k_pas=0.03, v_ref=4.0, tpwb = 0.0, tpww = 0.09, delay=0.08, v=0.0) self.set_chamber_comp('rv', E_pas= 0.1, E_act= 1.15, k_pas=0.03, v_ref=10., tr =0.30, td=0.45, v=50.) @@ -71,7 +71,7 @@ def __init__(self, name='MixedHeartMaynard4eWindkessel_parameters') -> None: self.set_activation_function('lv', af=activation_function_2) self.set_activation_function('rv', af=activation_function_2) - + self.set_activation_function('la', af=activation_function_3) self.set_activation_function('ra', af=activation_function_3) @@ -79,19 +79,19 @@ def __init__(self, name='MixedHeartMaynard4eWindkessel_parameters') -> None: # systemic circulation self.set_rlc_comp('sa', r=0.05, c=1.6 , l=0.0017 , v=450.0, v_ref=0.0) self.set_rlc_comp('sv', r=0.075, c=20.5, v=0.0, v_ref=0.0) - + # set impedances self.set_resistance('sai', r = 0.003) self.set_resistance('pai', r = 0.002) - + # pulmonary circulation self.set_rlc_comp('pa', r=0.01, c=3.8 , l=0.0017 , v=250.0, v_ref=0.0) self.set_rlc_comp('pv', r=0.006, c=20.5 , v=0.0, v_ref=0.0) - + # set capilary resistances self.set_resistance('sc', r = 0.5 + 0.52) self.set_resistance('pc', r = 0.05 + 0.25) - + # valves ##################################################### # self.set_valve_comp('ao', r=0.01, max_func=relu_max) @@ -114,19 +114,19 @@ def __init__(self, name='MixedHeartMaynard4eWindkessel_parameters') -> None: self.set_valve_comp('po', CQ=350., RRA=0.0) self.set_valve_comp('ti', CQ=400., RRA=0.0) ##################################################### - - + + def set_chamber_comp(self, key, **kwargs): self._set_comp(key=key, set=CHAMBERS, **kwargs) - + def set_activation_function(self, key, af): self._set_comp(key, set=CHAMBERS, af=af) - + def set_rlc_comp(self, key, **kwargs): self._set_comp(key=key, set=VESSELS, **kwargs) - + def set_valve_comp(self, key, **kwargs): self._set_comp(key=key, set=VALVES, **kwargs) - + def set_resistance(self, key, **kwargs): self._set_comp(key=key, set=IMPEDANCES + CAPILARIES, **kwargs) diff --git a/src/ModularCirc/Models/NaghaviModel.py b/src/ModularCirc/Models/NaghaviModel.py index 9e88a1a..2f46d71 100644 --- a/src/ModularCirc/Models/NaghaviModel.py +++ b/src/ModularCirc/Models/NaghaviModel.py @@ -7,56 +7,56 @@ class NaghaviModel(OdeModel): def __init__(self, time_setup_dict, parobj:NaghaviModelParameters=NaghaviModelParameters(), suppress_printing:bool=False) -> None: super().__init__(time_setup_dict) self.name = 'NaghaviModel' - + if not suppress_printing: print(parobj) - + # Defining the aorta object self.components['ao'] = Rlc_component( name='Aorta', - time_object=self.time_object, + time_object=self.time_object, r = parobj['ao']['r'], - c = parobj['ao']['c'], + c = parobj['ao']['c'], l = parobj['ao']['l'], v_ref = parobj['ao']['v_ref'], v = parobj['ao']['v'], ) self.set_v_sv('ao') - + # Defining the arterial system object self.components['art'] = Rlc_component(name='Arteries', time_object=self.time_object, r = parobj['art']['r'], - c = parobj['art']['c'], + c = parobj['art']['c'], l = parobj['art']['l'], v_ref = parobj['art']['v_ref'], v = parobj['art']['v'], ) self.set_v_sv('art') - + # Defining the venous system object self.components['ven'] = Rlc_component(name='VenaCava', time_object=self.time_object, r = parobj['ven']['r'], - c = parobj['ven']['c'], + c = parobj['ven']['c'], l = parobj['ven']['l'], v_ref = parobj['ven']['v_ref'], v = parobj['ven']['v'], ) self.set_v_sv('ven') - + # Defining the aortic valve object self.components['av'] = Valve_non_ideal(name='AorticValve', time_object=self.time_object, r=parobj['av']['r'], max_func=parobj['av']['max_func'] ) - + # Defining the mitral valve object self.components['mv'] = Valve_non_ideal(name='MitralValve', time_object=self.time_object, r=parobj['mv']['r'], max_func=parobj['mv']['max_func'] ) - + # Defining the left atrium activation function # def la_af(t, t_max=parobj['la']['t_max'], t_tr=parobj['la']['t_tr'], tau=parobj['la']['tau'], af=parobj['la']['activation_function']): # return af(time_shift(t, 100., self.time_object), t_max=t_max, t_tr=t_tr, tau=tau) @@ -69,18 +69,18 @@ def __init__(self, time_setup_dict, parobj:NaghaviModelParameters=NaghaviModelPa k_pas=parobj['la']['k_pas'], v =parobj['la']['v'], af =parobj['la']['activation_function'], - t_tr =parobj['la']['t_tr'], + t_tr =parobj['la']['t_tr'], t_max=parobj['la']['t_max'], - tau =parobj['la']['tau'], + tau =parobj['la']['tau'], delay=parobj['la']['delay'] ) self._state_variable_dict['v_la'] = self.components['la']._V self._state_variable_dict['v_la'].set_name('v_la') self.all_sv_data['v_la'] = self.components['la'].V self.components['la']._V._u = self.all_sv_data['v_la'] - + # Defining the left ventricle activation function - self.components['lv'] = HC_mixed_elastance(name='LeftVentricle', + self.components['lv'] = HC_mixed_elastance(name='LeftVentricle', time_object=self.time_object, E_pas=parobj['lv']['E_pas'], E_act=parobj['lv']['E_act'], @@ -88,7 +88,7 @@ def __init__(self, time_setup_dict, parobj:NaghaviModelParameters=NaghaviModelPa v_ref=parobj['lv']['v_ref'], v =parobj['lv']['v'], af =parobj['lv']['activation_function'], - t_tr =parobj['la']['t_tr'], + t_tr =parobj['la']['t_tr'], t_max=parobj['la']['t_max'], tau =parobj['la']['tau'], delay=parobj['lv']['delay'] @@ -97,48 +97,48 @@ def __init__(self, time_setup_dict, parobj:NaghaviModelParameters=NaghaviModelPa self._state_variable_dict['v_lv'].set_name('v_lv') self.all_sv_data['v_lv'] = self.components['lv'].V self.components['lv']._V._u = self.all_sv_data['v_lv'] - + for component in self.components.values(): component.setup() - + # connect the left ventricle class to the aortic valve - self.connect_modules(self.components['lv'], - self.components['av'], - plabel='p_lv', + self.connect_modules(self.components['lv'], + self.components['av'], + plabel='p_lv', qlabel='q_av', ) # connect the aortic valve to the aorta - self.connect_modules(self.components['av'], - self.components['ao'], - plabel='p_ao', + self.connect_modules(self.components['av'], + self.components['ao'], + plabel='p_ao', qlabel='q_av') # connect the aorta to the arteries - self.connect_modules(self.components['ao'], - self.components['art'], - plabel= 'p_art', + self.connect_modules(self.components['ao'], + self.components['art'], + plabel= 'p_art', qlabel= 'q_ao') # connect the arteries to the veins - self.connect_modules(self.components['art'], - self.components['ven'], - plabel= 'p_ven', + self.connect_modules(self.components['art'], + self.components['ven'], + plabel= 'p_ven', qlabel= 'q_art') # connect the veins to the left atrium - self.connect_modules(self.components['ven'], - self.components['la'], - plabel= 'p_la', + self.connect_modules(self.components['ven'], + self.components['la'], + plabel= 'p_la', qlabel='q_ven', ) # connect the left atrium to the mitral valve - self.connect_modules(self.components['la'], - self.components['mv'], - plabel= 'p_la', + self.connect_modules(self.components['la'], + self.components['mv'], + plabel= 'p_la', qlabel='q_mv') # connect the mitral valve to the left ventricle - self.connect_modules(self.components['mv'], - self.components['lv'], - plabel='p_lv', + self.connect_modules(self.components['mv'], + self.components['lv'], + plabel='p_lv', qlabel='q_mv', ) - + for component in self.components.values(): - component.setup() \ No newline at end of file + component.setup() diff --git a/src/ModularCirc/Models/NaghaviModelParameters.py b/src/ModularCirc/Models/NaghaviModelParameters.py index 8df07bb..48e0175 100644 --- a/src/ModularCirc/Models/NaghaviModelParameters.py +++ b/src/ModularCirc/Models/NaghaviModelParameters.py @@ -4,7 +4,7 @@ TEMPLATE_TIME_SETUP_DICT = { "name": "TimeTest", - "ncycles": 30, + "ncycles": 30, "tcycle": 1000., "dt": 1., "export_min": 1, @@ -20,27 +20,27 @@ def __init__(self, name='Naghavi Model') -> None: for key in ['av', 'mv']: self.components[key] = pd.Series(index=['r', 'max_func'], dtype=object) for key in ['la', 'lv']: - self.components[key] = pd.Series(index=['E_pas', - 'E_act', - 'v_ref', + self.components[key] = pd.Series(index=['E_pas', + 'E_act', + 'v_ref', 'k_pas', - 'activation_function', - 't_tr', + 'activation_function', + 't_tr', 't_max', - 'tau', + 'tau', 'delay', 'v', 'p', 'activation_func'], dtype=object) - + self.set_rlc_comp(key='ao', r=240., c=0.3, l=0.0, v_ref=100., v=0.025*5200.0, p=None) self.set_rlc_comp(key='art', r=1125., c=3.0, l=0.0, v_ref=900., v=0.21 *5200.0, p=None) self.set_rlc_comp(key='ven', r=9.0 , c=133.3,l=0.0, v_ref=2800., v=0.727*5200.0, p=None) - + self.set_valve_comp(key='av', r=6. , max_func=relu_max) - self.set_valve_comp(key='mv', r=4.1, max_func=relu_max) - - + self.set_valve_comp(key='mv', r=4.1, max_func=relu_max) + + # original self.set_chamber_comp('la', E_pas=0.44, E_act=0.45, v_ref=10., k_pas=0.05, # 0.05 activation_function=activation_function_1, @@ -49,19 +49,19 @@ def __init__(self, name='Naghavi Model') -> None: self.set_chamber_comp('lv', E_pas=1.0, E_act=3.0, v_ref=10., k_pas=0.027, # 0.027 activation_function=activation_function_1, t_tr=420., t_max=280., tau=25., delay=None, v=0.02*5200.0, p=None) - + def set_rc_comp(self, key:str, **kwargs): self._set_comp(key=key, set=['ao','art', 'ven'], **kwargs) - + def set_rlc_comp(self, key:str, **kwargs): self._set_comp(key=key, set=['ao','art', 'ven'], **kwargs) - + def set_valve_comp(self, key:str, **kwargs): self._set_comp(key=key, set=['av', 'mv'], **kwargs) - - + + def set_chamber_comp(self, key:str, **kwargs): self._set_comp(key=key, set=['lv', 'la'], **kwargs) - + def set_activation_function(self, key:str, activation_func=activation_function_2): - self._set_comp(key=key, set=['lv', 'la'], activation_func=activation_func) \ No newline at end of file + self._set_comp(key=key, set=['lv', 'la'], activation_func=activation_func) diff --git a/src/ModularCirc/Models/OdeModel.py b/src/ModularCirc/Models/OdeModel.py index 6925e0e..d273f63 100644 --- a/src/ModularCirc/Models/OdeModel.py +++ b/src/ModularCirc/Models/OdeModel.py @@ -11,19 +11,19 @@ def __init__(self, time_setup_dict) -> None: self.all_sv_data = pd.DataFrame(index=self.time_object.time.index, dtype='float64') self.components = dict() self.name = 'Template' - - - def connect_modules(self, - module1:ComponentBase, + + + def connect_modules(self, + module1:ComponentBase, module2:ComponentBase, pvariable:StateVariable = None, qvariable:StateVariable = None, - plabel:str=None, + plabel:str=None, qlabel:str=None - ) ->None: + ) ->None: """ Method for connecting Component modules. - + Inputs ------ module1 (Component) : upstream Component @@ -41,7 +41,7 @@ def connect_modules(self, else: module2._Q_i = qvariable module1._Q_o = qvariable - + if pvariable is None: if module1._P_o._ode_sys_mapping['u_func'] is not None or module1._P_o._ode_sys_mapping['dudt_func'] is not None: module2._P_i = module1._P_o @@ -52,7 +52,7 @@ def connect_modules(self, else: module2._P_i = pvariable module1._P_o = pvariable - + if plabel is not None and plabel not in self._state_variable_dict.keys(): module1._P_o.set_name(plabel) self._state_variable_dict[plabel] = module1._P_o @@ -62,12 +62,12 @@ def connect_modules(self, self._state_variable_dict[qlabel] = module1._Q_o self.all_sv_data[qlabel] = module1.Q_o return - - + + @property def state_variable_dict(self): return self._state_variable_dict - + def __str__(self) -> str: out = f'Model {self.name} \n\n' for component in self.components.values(): @@ -77,15 +77,15 @@ def __str__(self) -> str: for key in self._state_variable_dict.keys(): out += f" - {key} \n" return out - + def set_v_sv(self, comp_key:str): v_key = 'v_' + comp_key self._state_variable_dict[v_key] = self.components[comp_key]._V self._state_variable_dict[v_key].set_name(v_key) self.all_sv_data[v_key] = self.components[comp_key].V ##### test # self.components[comp_key]._V._u = self.all_sv_data[v_key] - - + + def __del__(self): if hasattr(self, 'name'): del self.name @@ -98,4 +98,4 @@ def __del__(self): del sv del self._state_variable_dict del self.time_object - return \ No newline at end of file + return diff --git a/src/ModularCirc/Models/ParametersObject.py b/src/ModularCirc/Models/ParametersObject.py index 2616cd7..2d709dc 100644 --- a/src/ModularCirc/Models/ParametersObject.py +++ b/src/ModularCirc/Models/ParametersObject.py @@ -11,11 +11,11 @@ def __init__(self, name='ParametersObject') -> None: self._cap = [] self._imp = [] self.timings = {} - - + + def __repr__(self) -> str: out = self._name + ' parameters set: \n' - for comp in self.components.keys(): + for comp in self.components.keys(): out += f" * Component - {bold_text(str(comp))}" + '\n' for key, item in self.components[comp].to_dict().items(): if isinstance(item, float): @@ -24,13 +24,13 @@ def __repr__(self) -> str: out += (f" - {bold_text(str(key)):<20} : {item}") + '\n' out += '\n' return out - + def __getitem__(self, key): return self.components[key] - + def __setitem__(self, key, val): self.components[key] = val - + def _set_comp(self, key:str, set=None, **kwargs): if key not in set and set is not None: raise Exception(f'Wrong key! {key} not in {set}.') @@ -38,7 +38,7 @@ def _set_comp(self, key:str, set=None, **kwargs): if val is None: continue assert k in self[key].index.values, f" {key}: {k} not in {self[key].index.values}" self[key].loc[k] = val - + def set_timings(self, key:str, **kwargs): if key not in self._chambers: raise Exception(f'Wrong key! {key} not in {self._chambers}.') @@ -46,6 +46,6 @@ def set_timings(self, key:str, **kwargs): if val is None: continue assert k in self.timings[key].index.values, f" {key}: {k} not in {self[key].index.values}" self[key].loc[k] = val - + def add_parameter_to_component(self, key:str, par:str, val=0.0): - self.components[key][par] = val \ No newline at end of file + self.components[key][par] = val diff --git a/src/ModularCirc/Models/README.md b/src/ModularCirc/Models/README.md index 9c1b30d..2a8fac6 100644 --- a/src/ModularCirc/Models/README.md +++ b/src/ModularCirc/Models/README.md @@ -53,7 +53,7 @@ This model is comprised of the following components: - tricuspid valve: simple Bernoulli model - **2 parameters** - right ventricle: linear time-varying elastance model - - **5 parameters** + - **5 parameters** - pulmonary valve: simple Bernoulli model - **2 parameters** - pulmonary artery sinus (RLC 3 component windkessel) @@ -109,4 +109,4 @@ Here, the following components were used: - valves: **Maynard** valves (4 parameters) - arteries: RLC with sinuses represented as a separate R component and capilaries similarly by R. This equivalent to the systemic and pulmonary arterial system being simulated as 2 4-component Windkessels in series with one another. (7/8 parameters systemic + 7/8 parameters pulmonary) -**Total parameter count: 58 or 56 (if you choose to add the capilary resistance to the resistance of the previous component)** \ No newline at end of file +**Total parameter count: 58 or 56 (if you choose to add the capilary resistance to the resistance of the previous component)** diff --git a/src/ModularCirc/Solver.py b/src/ModularCirc/Solver.py index d75936f..4a518ca 100644 --- a/src/ModularCirc/Solver.py +++ b/src/ModularCirc/Solver.py @@ -21,10 +21,10 @@ import warnings class Solver(): - def __init__(self, + def __init__(self, model:OdeModel=None, ) -> None: - + self.model = model @@ -87,12 +87,12 @@ def __init__(self, # Number of sub-iterations for the solver. self._n_sub_iter = 1 - + # flag for checking if the model is converged or not... self.converged = False - - - def setup(self, + + + def setup(self, optimize_secondary_sv:bool=False, suppress_output:bool=False, step_tol:float=1e-2, @@ -104,7 +104,7 @@ def setup(self, )->None: """ Method for detecting which are the principal variables and which are the secondary ones. - + ## Inputs optimize_secondary_sv : boolean flag used to switch on the optimization for secondary variable computations, this flag needs to be @@ -120,14 +120,14 @@ def setup(self, # Loop over the state variables and check if they have an update function, - # This code ensures that each state variable's update function is correctly assigned and indexed, + # This code ensures that each state variable's update function is correctly assigned and indexed, # allowing the solver to update the state variables during the simulation for key, component in self._vd.items(): - + # Get the index of the state variable. mkey = self._global_sv_id[key] # _global_sv_id - maps the state variable names to their indexes. - - # initialization function for a state variable. This function is used to initialize the state + + # initialization function for a state variable. This function is used to initialize the state # variable at the beginning of the simulation. if component.i_func is not None: if not suppress_output: print(f" -- Variable {bold_text(key)} added to the init list.") @@ -136,8 +136,8 @@ def setup(self, self._initialize_by_function[key] = component self._global_sv_init_fun[mkey] = component.i_func self._global_sv_init_ind[mkey] = [self._global_sv_id[key2] for key2 in component.i_inputs.to_list()] - - # derivative function for a state variable. This function is used to update the state variable + + # derivative function for a state variable. This function is used to update the state variable # during the numerical integration process. if component.dudt_func is not None: if not suppress_output: print(f" -- Variable {bold_text(key)} added to the principal variable key list.") @@ -148,14 +148,14 @@ def setup(self, self._global_psv_update_ind[mkey] = [self._global_sv_id[key2] for key2 in component.inputs.to_list()] # Pad the index array to the length of the state variable array. - self._global_psv_update_ind[mkey] = np.pad(self._global_psv_update_ind[mkey], - (0, self._N_sv-len(self._global_psv_update_ind[mkey])), + self._global_psv_update_ind[mkey] = np.pad(self._global_psv_update_ind[mkey], + (0, self._N_sv-len(self._global_psv_update_ind[mkey])), mode='constant', constant_values=-1) - + # Add the state variable name to the global primary state variable list. self._global_psv_names.append(key) - - # updated function for the secondary state variable. This function is used to update the state variable + + # updated function for the secondary state variable. This function is used to update the state variable # based on the current values of the primary state variables using algebraic relationships. elif component.u_func is not None: if not suppress_output: print(f" -- Variable {bold_text(key)} added to the secondary variable key list.") @@ -164,8 +164,8 @@ def setup(self, self._global_ssv_update_fun[mkey] = component.u_func self._global_ssv_update_fun_n[mkey] = component.u_name self._global_ssv_update_ind[mkey] = [self._global_sv_id[key2] for key2 in component.inputs.to_list()] - self._global_ssv_update_ind[mkey] = np.pad(self._global_ssv_update_ind[mkey], - (0, self._N_sv-len(self._global_ssv_update_ind[mkey])), + self._global_ssv_update_ind[mkey] = np.pad(self._global_ssv_update_ind[mkey], + (0, self._N_sv-len(self._global_ssv_update_ind[mkey])), mode='constant', constant_values=-1) else: continue @@ -173,8 +173,8 @@ def setup(self, if not suppress_output: print(' ') self.generate_dfdt_functions() - - + + if self._conv_cols is None: # If no specific columns for convergence (_conv_cols) are provided, # automatically select columns from the DataFrame (_asd) whose names @@ -185,8 +185,8 @@ def setup(self, self._cols = self._conv_cols # End the method without returning any specific value. - return None - + return None + def generate_dfdt_functions(self): @@ -203,39 +203,39 @@ def initialize_by_function(y:np.ndarray[float]) -> np.ndarray[float]: """ Initialize the state variables using a set of initialization functions. - This function applies a list of initialization functions (`funcs1`) to - specific subsets of the input array `y`, as determined by the indices - in `ids1`. Each function is called with `t=0.0` and the corresponding + This function applies a list of initialization functions (`funcs1`) to + specific subsets of the input array `y`, as determined by the indices + in `ids1`. Each function is called with `t=0.0` and the corresponding subset of `y`, and the results are combined into a single NumPy array. - The input array `y` is usually the initial state variable array, so + The input array `y` is usually the initial state variable array, so the 0th row of the self._asd DataFrame. Args: - y (np.ndarray[float]): A 1D NumPy array representing the state - variables to be initialized. Each subset of `y` is passed to + y (np.ndarray[float]): A 1D NumPy array representing the state + variables to be initialized. Each subset of `y` is passed to the corresponding initialization function. Returns: - np.ndarray[float]: A 1D NumPy array containing the initialized + np.ndarray[float]: A 1D NumPy array containing the initialized state variables, with the same length as the input array `y`. Note: - - Each function in `funcs1` is expected to accept two arguments: - `t` (a float, representing time) and `y` (a NumPy array, + - Each function in `funcs1` is expected to accept two arguments: + `t` (a float, representing time) and `y` (a NumPy array, representing the subset of state variables). - + Example use: >>> initialize_by_function(y=self._asd.iloc[0].to_numpy()) """ - return np.fromiter([fun(t=0.0, y=y[inds]) for fun, inds in zip(funcs1, ids1)], + return np.fromiter([fun(t=0.0, y=y[inds]) for fun, inds in zip(funcs1, ids1)], dtype=np.float64) # Function to update the secondary state variables based on the primary state variables. funcs2 = np.array(list(self._global_ssv_update_fun.values())) ids2 = np.stack(list(self._global_ssv_update_ind.values())) - - # @nb.njit(cache=True) + + # @nb.njit(cache=True) def s_u_update(t, y:np.ndarray[float]) -> np.ndarray[float]: """ Updates the secondary state variables based on the current values of the primary state variables. @@ -250,21 +250,21 @@ def s_u_update(t, y:np.ndarray[float]) -> np.ndarray[float]: Example use: >>> s_u_update(t=0.0, y=self._asd.iloc[0].to_numpy()) """ - return np.fromiter([fi(t=t, y=yi) for fi, yi in zip(funcs2, y[ids2])], + return np.fromiter([fi(t=t, y=yi) for fi, yi in zip(funcs2, y[ids2])], dtype=np.float64) - + def s_u_residual(y, yall, keys): """ Function to compute the residual of the secondary state variables.""" yall[keys] = y return (y - s_u_update(0.0, yall)) - + def optimize(y:np.ndarray, keys): """ Function to optimize the secondary state variables.""" yk = y[keys] sol = least_squares( # root - s_u_residual, - yk, - args=(y, keys), + s_u_residual, + yk, + args=(y, keys), ftol=1.0e-5, xtol=1.0e-15, loss='linear', @@ -285,7 +285,7 @@ def optimize(y:np.ndarray, keys): # indexes of the primary state variables dependencies. ids3 = np.stack(list(self._global_psv_update_ind.values())) - + T = self._to.tcycle N_zeros_0 = len(self._global_sv_id) _n_sub_iter = self._n_sub_iter @@ -316,9 +316,9 @@ def optimize(y:np.ndarray, keys): keys3_dict2[key].update({val,}) else: keys3_dict2[key].update(keys4_dict[val]) - + sparsity_map = dict() - + for i, key in enumerate(keys3): sparsity_map[i] = set() for val in keys3_dict2[key]: @@ -335,7 +335,7 @@ def optimize(y:np.ndarray, keys): perm_mat = np.zeros((len(perm), len(perm))) for i,j in enumerate(perm): perm_mat[i,j] = 1 - + self.perm_mat = perm_mat # reorders the sparse matrix to reduce the bandwidth @@ -349,7 +349,7 @@ def optimize(y:np.ndarray, keys): self.lband = lband self.uband = uband - + def pv_dfdt_update(t, y:np.ndarray[float]) -> np.ndarray[float]: """ Function to compute the derivatives of the primary state variables over time.""" @@ -372,20 +372,20 @@ def pv_dfdt_update(t, y:np.ndarray[float]) -> np.ndarray[float]: # updates the secondary state variables, and optimises them if necessary for _ in range(_n_sub_iter): - y_temp[keys4] = s_u_update(t, y_temp) + y_temp[keys4] = s_u_update(t, y_temp) if _optimize_secondary_sv: y_temp[keys4] = optimize(y_temp, keys4) # returns the derivatives of the primary state variables, reordered back to the original order return perm_mat @ np.fromiter([fi(t=ht, y=yi) for fi, yi in zip(funcs3, y_temp[ids3])], dtype=np.float64) - - - self.initialize_by_function = initialize_by_function + + + self.initialize_by_function = initialize_by_function self.pv_dfdt_global = pv_dfdt_update self.s_u_update = s_u_update - + self.optimize = optimize self.s_u_residual = s_u_residual - + def advance_cycle(self, y0, cycleID, step = 1): @@ -399,9 +399,9 @@ def advance_cycle(self, y0, cycleID, step = 1): with warnings.catch_warnings(): warnings.simplefilter("ignore") if self._method != 'LSODA': - res = solve_ivp(fun=self.pv_dfdt_global, - t_span=(t[0], t[-1]), - y0=self.perm_mat @ y0, + res = solve_ivp(fun=self.pv_dfdt_global, + t_span=(t[0], t[-1]), + y0=self.perm_mat @ y0, t_eval=t, max_step=self.dt, method=self._method, @@ -409,9 +409,9 @@ def advance_cycle(self, y0, cycleID, step = 1): rtol=self._rtol, ) else: - res = solve_ivp(fun=self.pv_dfdt_global, - t_span=(t[0], t[-1]), - y0=self.perm_mat @ y0, + res = solve_ivp(fun=self.pv_dfdt_global, + t_span=(t[0], t[-1]), + y0=self.perm_mat @ y0, t_eval=t, method=self._method, atol=self._atol, @@ -431,29 +431,29 @@ def advance_cycle(self, y0, cycleID, step = 1): ids = list(self._global_psv_update_fun.keys()) inds= list(range(len(ids))) self._asd.iloc[cycleID*n_t:(end_cycle)*n_t+1, ids] = y[inds, 0:n_t*step+1].T - + if cycleID == 0: return False cycleP = end_cycle - 1 cs = self._asd[self._cols].iloc[cycleP*n_t:end_cycle*n_t, :].values cp = self._asd[self._cols].iloc[(cycleP-1) *n_t:(cycleP)*n_t, :].values - + cp_ptp = np.max(np.abs(cp), axis=0) cp_r = np.max(np.abs(cs - cp), axis=0) - + test = cp_r / cp_ptp test[cp_ptp <= 1e-10] = cp_r[cp_ptp <= 1e-10] if np.max(test) > self._step_tol : return False return True - - + + def solve(self): - + # initialize the solution fields self._asd.loc[0, self._initialize_by_function.index] = \ self.initialize_by_function(y=self._asd.loc[0].to_numpy()).T - + # Solve the main system of ODEs.. for i in range(0, self._to.ncycles, self.step): # step is a pulse, we might wabnt to do it in all pulses @@ -474,66 +474,55 @@ def solve(self): if i + self.step - 1 == self._to.ncycles - 1: self._Nconv = i + self.step - 1 self.converged = False - + self._to.n_t = (self.Nconv+1)*(self._to.n_c-1) + 1 - + self._asd = self._asd.iloc[:self._to.n_t] self._to._sym_t = self._to._sym_t.head(self._to.n_t) self._to._cycle_t = self._to._cycle_t.head(self._to.n_t) - - - keys4 = np.array(list(self._global_ssv_update_fun.keys())) - temp = np.zeros(self._asd.iloc[:,keys4].shape) + + + keys4 = np.array(list(self._global_ssv_update_fun.keys())) + temp = np.zeros(self._asd.iloc[:,keys4].shape) for i, line in enumerate(self._asd.values) : line[keys4] = self.s_u_update(0.0, line) if self._optimize_secondary_sv: temp[i,:] = self.optimize(line, keys4) else: - temp[i,:] = line[keys4] - self._asd.iloc[:,keys4] = temp - + temp[i,:] = line[keys4] + self._asd.iloc[:,keys4] = temp + for key in self._vd.keys(): - self._vd[key]._u = self._asd[key] + self._vd[key]._u = self._asd[key] + - @property def vd(self) -> Series[StateVariable]: return self._vd - + @property def dt(self) -> float: return self._to.dt - + @property def Nconv(self) -> float: return self._Nconv - + @property def optimize_secondary_sv(self)->bool: return self._optimize_secondary_sv - - + + @property def n_sub_iter(self)->int: - return self._n_sub_iter - - + return self._n_sub_iter + + @n_sub_iter.setter def n_sub_iter(self, value): assert isinstance(value, int) assert value > 0 self._n_sub_iter = value - - - - - - - - - - - \ No newline at end of file diff --git a/src/ModularCirc/StateVariable.py b/src/ModularCirc/StateVariable.py index d93a827..586e9ae 100644 --- a/src/ModularCirc/StateVariable.py +++ b/src/ModularCirc/StateVariable.py @@ -9,7 +9,7 @@ def __init__(self, name:str, timeobj:TimeClass) -> None: self._to = timeobj self._u = pd.Series(np.zeros((timeobj.n_t,)), name=name, dtype='float64') self._cv = 0.0 - + self._ode_sys_mapping = pd.Series({ 'dudt_func' : None, 'u_func' : None, @@ -19,75 +19,75 @@ def __init__(self, name:str, timeobj:TimeClass) -> None: 'i_func' : None, # initialization function 'i_inputs' : None, # initialization function inputs }) - + def __repr__(self) -> str: return f" > variable name: {self._name}" - + def set_dudt_func(self, function, function_name:str)->None: self._ode_sys_mapping['dudt_func'] = function self._ode_sys_mapping['dudt_name'] = function_name return - + def set_u_func(self, function, function_name:str)->None: self._ode_sys_mapping['u_func'] = function self._ode_sys_mapping['u_name'] = function_name return - + def set_inputs(self, inputs:Series[str]): self._ode_sys_mapping['inputs'] = inputs return - + def set_name(self, name)->None: self._name = name return - + def set_i_func(self, function, function_name:str)->None: self._ode_sys_mapping['i_func'] = function self._ode_sys_mapping['i_name'] = function_name - + def set_i_inputs(self, inputs:Series[str])->None: self._ode_sys_mapping['i_inputs'] = inputs - + @property def name(self): return self._name - + @property def dudt_func(self): return self._ode_sys_mapping['dudt_func'] - + @property def dudt_name(self) -> str: return self._ode_sys_mapping['dudt_name'] - + @property def inputs(self) -> Series[str]: return self._ode_sys_mapping['inputs'] - + @property def u_func(self): return self._ode_sys_mapping['u_func'] - + @property def u_name(self) -> str: return self._ode_sys_mapping['u_name'] - + @property def u(self) -> list[float]: return self._u - + @property def i_func(self): return self._ode_sys_mapping['i_func'] - + @property def i_name(self): return self._ode_sys_mapping['i_name'] - + @property def i_inputs(self): return self._ode_sys_mapping['i_inputs'] - + def __del__(self): if hasattr(self, '_name'): del self._name @@ -105,6 +105,4 @@ def __del__(self): del subitem del item # print(self._ode_sys_mapping) - del self._ode_sys_mapping - - \ No newline at end of file + del self._ode_sys_mapping diff --git a/src/ModularCirc/Time.py b/src/ModularCirc/Time.py index 72e07ae..252f439 100644 --- a/src/ModularCirc/Time.py +++ b/src/ModularCirc/Time.py @@ -14,35 +14,35 @@ def __init__(self, time_setup_dict) -> None: self._time_setup_dict = time_setup_dict self._initialize_time_array() self.cti = 0 # current time step index - + @property def ncycles(self): if 'ncycles' in self._time_setup_dict.keys(): return self._time_setup_dict['ncycles'] - else: + else: return None - - @property + + @property def tcycle(self): if 'tcycle' in self._time_setup_dict.keys(): return self._time_setup_dict['tcycle'] else: return None - + @property def dt(self): if 'dt' in self._time_setup_dict.keys(): return self._time_setup_dict['dt'] else: return None - + @property def export_min(self): if 'export_min' in self._time_setup_dict.keys(): return self._time_setup_dict['export_min'] else: return None - + def _initialize_time_array(self): # discretization of on heart beat, used as template self._one_cycle_t = pd.Series(np.linspace( @@ -51,25 +51,25 @@ def _initialize_time_array(self): num = int(self.tcycle / self.dt)+1, dtype= np.float64 )) - + # discretization of the entire simulation duration self._sym_t = pd.Series( [t+cycle*self.tcycle for cycle in range(self.ncycles) for t in self._one_cycle_t[:-1]] + [self._one_cycle_t.values[-1]+(self.ncycles-1)*self.tcycle,] ) - + # array of the current time within the heart cycle self._cycle_t = pd.Series( [t for _ in range(self.ncycles) for t in self._one_cycle_t[:-1]] + [self._one_cycle_t.values[-1],] ) - + self.time = pd.DataFrame({'cycle_t' : self._cycle_t, 'sym_t' : self._sym_t}) - + # the total number of time steps including initial time step self.n_t = len(self._sym_t) - + # the number of time steps in a cycle self.n_c = len(self._one_cycle_t) return - + def new_time_step(self): - self.cti += 1 \ No newline at end of file + self.cti += 1 diff --git a/src/ModularCirc/_BatchRunner.py b/src/ModularCirc/_BatchRunner.py index 28fd924..03935f3 100644 --- a/src/ModularCirc/_BatchRunner.py +++ b/src/ModularCirc/_BatchRunner.py @@ -24,21 +24,21 @@ def __init__(self, sampler:str='LHS', seed=DEFAULT_RANDOM_SEED) -> None: self._sample_generator = sampler_dictionary[sampler] self._seed = seed return - + def setup_sampler(self, template_json): with open(template_json) as file: self._template : dict = json.load(file) self._parameters_2_sample = dict() self._parameters_constant = dict() - + self.condense_dict_parameters(self._template) self._d = len(self._parameters_2_sample) self._l_bounds = [bounds[0] for bounds in self._parameters_2_sample.values()] self._u_bounds = [bounds[1] for bounds in self._parameters_2_sample.values()] - - self._sampler : QMCEngine = self._sample_generator(d=self._d, seed=self._seed) + + self._sampler : QMCEngine = self._sample_generator(d=self._d, seed=self._seed) return - + def condense_dict_parameters(self, dict_param, prev=''): for key, val in dict_param.items(): if len(prev) > 0: @@ -53,8 +53,8 @@ def condense_dict_parameters(self, dict_param, prev=''): self._parameters_2_sample[new_key] = tuple(np.array(r) * value) else: self._parameters_constant[new_key] = val[0] - return - + return + def sample(self, nsamples:int): samples = self._sampler.random(nsamples) samples_scaled = scale(sample=samples, l_bounds=self._l_bounds, u_bounds=self._u_bounds) @@ -62,11 +62,11 @@ def sample(self, nsamples:int): for key, val in self._parameters_constant.items(): self._samples[key] = val return - + @property def samples(self): return self._samples.copy() - + def map_sample_timings(self, map:dict = dict(), ref_time=1.0): for key, mappings in map.items(): for key2 in mappings: @@ -75,27 +75,27 @@ def map_sample_timings(self, map:dict = dict(), ref_time=1.0): if key not in mappings: self._samples.drop(key, inplace=True, axis=1) self._ref_time = ref_time return - + def map_vessel_volume(self): vessels = list(self._template['VESSELS'].keys()) vessels_c_names =[vessel + '.c' for vessel in vessels] - + samples_vessels_c = self._samples[vessels_c_names] tot_vessels_c = samples_vessels_c.sum(axis=1) - + for vessel in vessels: self._samples[vessel + '.v'] = self._samples['v_tot'] * self._samples[vessel + '.c'] / tot_vessels_c self._samples.drop('v_tot', axis=1, inplace=True) return - - + + def setup_model(self, model:OdeModel, po:ParametersObject, time_setup:dict): self._model_generator = model self._po_generator = po self._tst = time_setup return - - + + def run_batch(self, n_jobs=1, **kwargs): if n_jobs == 1: success = self._samples.apply(lambda row : self._run_case(row, **kwargs), axis=0) @@ -105,12 +105,12 @@ def run_batch(self, n_jobs=1, **kwargs): for _, row in tqdm(self._samples.iterrows(), total=len(self._samples)) ) return success - + def _run_case(self, row, **kwargs): time_setup = self._tst.copy() time_setup['tcycle'] = row['T'] time_setup['dt']*= row['T'] / self._ref_time - + po : ParametersObject = self._po_generator() for key, val in row.items(): if key == "T": @@ -119,33 +119,33 @@ def _run_case(self, row, **kwargs): obj, param = key.split(".") except: raise Exception(key) - + po._set_comp(obj,[obj,], **{param: val}) model : OdeModel = self._model_generator(time_setup_dict = time_setup, parobj=po, suppress_printing=True) - + solver = Solver(model=model) - + if 'optimize_secondary_sv' in kwargs: optimize_secondary_sv = kwargs['optimize_secondary_sv'] else: optimize_secondary_sv = False - + if 'conv_cols' in kwargs: conv_cols = kwargs['conv_cols'] else: conv_cols = None - + if 'method' in kwargs: method = kwargs['method'] else: method = 'LSODA' - + if 'out_cols' in kwargs: out_cols = kwargs['out_cols'] else: out_cols = None - + if 'output_path' in kwargs: output_path = kwargs['output_path'] else: @@ -156,25 +156,21 @@ def _run_case(self, row, **kwargs): optimize_secondary_sv=optimize_secondary_sv, conv_cols=conv_cols, method=method - ) + ) solver.solve() - + if not solver.converged: return False - + if out_cols is None: raw_signal = solver._asd.copy() else: raw_signal = solver._asd[out_cols].copy() - + raw_signal_short = raw_signal.tail(model.time_object.n_c).copy() raw_signal_short.index = pd.MultiIndex.from_tuples([(row.name, i) for i in range(len(raw_signal_short))], names= ['realization', 'time_ind']) raw_signal_short.loc[:,'T'] = model.time_object._sym_t.values[-model.time_object.n_c:] - + if output_path is not None: raw_signal_short.loc[row.name].to_csv(os.path.join(output_path, f'all_outputs_{row.name}.csv')) - + return raw_signal_short - - - - diff --git a/src/ModularCirc/__init__.py b/src/ModularCirc/__init__.py index 510915f..55936d3 100644 --- a/src/ModularCirc/__init__.py +++ b/src/ModularCirc/__init__.py @@ -1 +1 @@ -from ._BatchRunner import _BatchRunner as BatchRunner \ No newline at end of file +from ._BatchRunner import _BatchRunner as BatchRunner diff --git a/tests/expected_outputs/KorakianitisMixedModel_expected_output.json b/tests/expected_outputs/KorakianitisMixedModel_expected_output.json index 7328d40..e89cd67 100644 --- a/tests/expected_outputs/KorakianitisMixedModel_expected_output.json +++ b/tests/expected_outputs/KorakianitisMixedModel_expected_output.json @@ -298,4 +298,4 @@ } } } -} \ No newline at end of file +} diff --git a/tests/expected_outputs/NaghaviModel_expected_output.json b/tests/expected_outputs/NaghaviModel_expected_output.json index 2a3a28f..617506d 100644 --- a/tests/expected_outputs/NaghaviModel_expected_output.json +++ b/tests/expected_outputs/NaghaviModel_expected_output.json @@ -158,4 +158,4 @@ } } } -} \ No newline at end of file +} diff --git a/tests/test_KorakianitisMixedModel.py b/tests/test_KorakianitisMixedModel.py index 0d21846..b7b776b 100644 --- a/tests/test_KorakianitisMixedModel.py +++ b/tests/test_KorakianitisMixedModel.py @@ -19,20 +19,20 @@ def setUp(self): """ Set up the test environment for the KorakianitisMixedModel. - This method initializes the necessary components and configurations + This method initializes the necessary components and configurations required for testing the KorakianitisMixedModel. It performs the following: - + - Sets a random seed for reproducibility. - Defines the base directory for file paths. - - Configures the simulation time setup parameters, including the number + - Configures the simulation time setup parameters, including the number of cycles, cycle duration, time step size, and minimum export cycles. - Initializes the parameter object for the KorakianitisMixedModel. - - Creates an instance of the KorakianitisMixedModel with the specified + - Creates an instance of the KorakianitisMixedModel with the specified time setup and parameter object, suppressing console output. - - Initializes the solver for the model and configures it with the + - Initializes the solver for the model and configures it with the "LSODA" method, suppressing console output. - Loads expected output values from a JSON file for validation during tests. - - Verifies the existence of the expected output file and raises an + - Verifies the existence of the expected output file and raises an assertion error if the file is not found. Raises: @@ -57,14 +57,14 @@ def setUp(self): # Initializing the parameter object self.parobj = KorakianitisMixedModel_parameters() - # Initializing the model - self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, - parobj=self.parobj, + # Initializing the model + self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, suppress_printing=True) - + # Initializing the solver self.solver = Solver(model=self.model) - + # Solver is being setup: switching off console printing and setting the solver method to "LSODA" self.solver.setup(suppress_output=True, method='LSODA',step=1) @@ -85,7 +85,7 @@ def test_model_initialization(self): self.assertIsInstance(self.model, KorakianitisMixedModel) # Verify model has attribute self.assertTrue(hasattr(self.solver.model, 'components')) - # Verify is a component + # Verify is a component self.assertIn('lv', self.solver.model.components) # Verify correct assignment of parameters from parobj to model self.assertEqual(self.solver.model.components['lv'].E_pas, self.parobj.components['lv']['E_pas']) @@ -105,9 +105,9 @@ def test_solver_run(self): 2. The solver converges successfully for each step size. 3. The state variables of the model components are updated correctly after solving. 4. The computed results match the expected values within a specified tolerance. - + """ - + cycle_step_sizes = [1, 3, 5, 7] # Define the step sizes to test for i_cycle_step_size in cycle_step_sizes: @@ -120,11 +120,11 @@ def test_solver_run(self): # Initializing the parameter object self.parobj = KorakianitisMixedModel_parameters() - # Initializing the model - self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, - parobj=self.parobj, + # Initializing the model + self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, suppress_printing=True) - + # Initializing the solver self.solver = Solver(model=self.model) @@ -132,10 +132,10 @@ def test_solver_run(self): self.solver.setup(suppress_output=True, method='LSODA', step=i_cycle_step_size) # Running the model - self.solver.solve() + self.solver.solve() # Verify the solver converged - self.assertTrue(self.solver.converged or self.solver._Nconv is not None) + self.assertTrue(self.solver.converged or self.solver._Nconv is not None) # Verifying the model changed the state variables stored within components. self.assertTrue(len(self.solver.model.components['lv'].V.values) > 0) @@ -161,8 +161,8 @@ def test_solver_run(self): [self.expected_values["results"][str(i_cycle_step_size)][key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()] ) new_ndarray = np.array([new_dict[key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()]) - test_ndarray = np.where(np.abs(expected_ndarray) > 1e-6, - np.abs((expected_ndarray - new_ndarray) / expected_ndarray), + test_ndarray = np.where(np.abs(expected_ndarray) > 1e-6, + np.abs((expected_ndarray - new_ndarray) / expected_ndarray), np.abs((expected_ndarray - new_ndarray))) self.assertTrue((test_ndarray < RELATIVE_TOLERANCE).all()) diff --git a/tests/test_NaghaviModel.py b/tests/test_NaghaviModel.py index c5ffbaf..ccbaa62 100644 --- a/tests/test_NaghaviModel.py +++ b/tests/test_NaghaviModel.py @@ -17,8 +17,8 @@ class TestNaghaviModel(unittest.TestCase): """ Unit tests for the NaghaviModel and its associated Solver. - This test suite verifies the correct initialization, configuration, and functionality of the - NaghaviModel and Solver classes. It ensures that the model behaves as expected under various + This test suite verifies the correct initialization, configuration, and functionality of the + NaghaviModel and Solver classes. It ensures that the model behaves as expected under various conditions and that the computed results match the expected outputs. Classes: TestNaghaviModel: A unittest.TestCase subclass containing tests for the NaghaviModel and Solver. @@ -47,7 +47,7 @@ def setUp(self): Raises: AssertionError: If the expected output file does not exist at the specified path. """ - + # Set a random seed for reproducibility np.random.seed(42) @@ -64,14 +64,14 @@ def setUp(self): self.parobj = NaghaviModelParameters() - self.model = NaghaviModel(time_setup_dict=self.time_setup_dict, - parobj=self.parobj, + self.model = NaghaviModel(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, suppress_printing=True) - + # Initializing the solver self.solver = Solver(model=self.model) - - # Solver is being setup: switching off console printing and setting the solver method to "LSODA" + + # Solver is being setup: switching off console printing and setting the solver method to "LSODA" self.solver.setup( optimize_secondary_sv=False, suppress_output=True, @@ -79,7 +79,7 @@ def setUp(self): conv_cols=['p_lv', 'v_lv'], method='LSODA', step=1 - ) + ) # Load expected values from a JSON file output_file_path = os.path.join(self.base_dir, 'expected_outputs', 'NaghaviModel_expected_output.json') @@ -98,7 +98,7 @@ def test_model_initialization(self): self.assertIsInstance(self.model, NaghaviModel) # Verify model has attribute self.assertTrue(hasattr(self.solver.model, 'components')) - # Verify is a component + # Verify is a component self.assertIn('lv', self.solver.model.components) # Verify correct assignment of parameters from parobj to model self.assertEqual(self.solver.model.components['lv'].E_pas, self.parobj.components['lv']['E_pas']) @@ -120,9 +120,9 @@ def test_solver_run(self): 2. The solver converges successfully for each step size. 3. The state variables of the model components are updated after solving. 4. The computed results match the expected values within a specified tolerance. - + """ - + cycle_step_sizes = [1, 3, 5, 7] # Define the step sizes to test for i_cycle_step_size in cycle_step_sizes: @@ -135,11 +135,11 @@ def test_solver_run(self): # Initializing the parameter object self.parobj = NaghaviModelParameters() - # Initializing the model - self.model = NaghaviModel(time_setup_dict=self.time_setup_dict, - parobj=self.parobj, + # Initializing the model + self.model = NaghaviModel(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, suppress_printing=True) - + # Initializing the solver self.solver = Solver(model=self.model) @@ -151,13 +151,13 @@ def test_solver_run(self): conv_cols=['p_lv', 'v_lv'], method='LSODA', step=i_cycle_step_size - ) + ) # Running the model - self.solver.solve() + self.solver.solve() # Verify the solver converged - self.assertTrue(self.solver.converged or self.solver._Nconv is not None) + self.assertTrue(self.solver.converged or self.solver._Nconv is not None) # Verifying the model changed the state variables stored within components. self.assertTrue(len(self.solver.model.components['lv'].V.values) > 0) @@ -183,8 +183,8 @@ def test_solver_run(self): [self.expected_values["results"][str(i_cycle_step_size)][key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()] ) new_ndarray = np.array([new_dict[key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()]) - test_ndarray = np.where(np.abs(expected_ndarray) > 1e-6, - np.abs((expected_ndarray - new_ndarray) / expected_ndarray), + test_ndarray = np.where(np.abs(expected_ndarray) > 1e-6, + np.abs((expected_ndarray - new_ndarray) / expected_ndarray), np.abs((expected_ndarray - new_ndarray))) self.assertTrue((test_ndarray < RELATIVE_TOLERANCE).all()) diff --git a/tests/test_Solver.py b/tests/test_Solver.py index 9dcad00..1979e63 100644 --- a/tests/test_Solver.py +++ b/tests/test_Solver.py @@ -16,31 +16,31 @@ class TestSolver(unittest.TestCase): """ - TestSolver is a unittest.TestCase class designed to test the functionality of the Solver class and its + TestSolver is a unittest.TestCase class designed to test the functionality of the Solver class and its integration with the KorakianitisMixedModel. Methods ------- setUp(): - Initializes the test environment, including setting a random seed, defining the time setup dictionary, + Initializes the test environment, including setting a random seed, defining the time setup dictionary, initializing the parameter object, model, and solver, and setting up the solver. test_solver_initialization(): Verifies that the solver is correctly initialized and that the model is correctly assigned to the solver. test_solver_setup(): Verifies the setup attributes of the solver and ensures that the generated functions are not None. test_initialize_by_function(): - Tests the initialize_by_function() method of the solver, ensuring it accepts the expected input and + Tests the initialize_by_function() method of the solver, ensuring it accepts the expected input and returns the expected output. test_optimize(): Tests the optimize() method of the solver, ensuring it accepts the expected input and can run with it. test_pv_dfdt_update(): - Tests the pv_dfdt_update() method of the solver, ensuring it accepts the expected input and the output + Tests the pv_dfdt_update() method of the solver, ensuring it accepts the expected input and the output matches the expected output. test_s_u_update(): - Tests the s_u_update() method of the solver, ensuring it accepts the expected input and the output + Tests the s_u_update() method of the solver, ensuring it accepts the expected input and the output matches the expected output. test_solver_solve(): - Tests the solve() method of the solver, ensuring the solver converges and the output matches the + Tests the solve() method of the solver, ensuring the solver converges and the output matches the expected values. """ @@ -80,22 +80,22 @@ def setUp(self): self.parobj = KorakianitisMixedModel_parameters() # Initialize the model - self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, - parobj=self.parobj, + self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, suppress_printing=True) - + # Initialize the solver self.solver = Solver(model=self.model) # Setup the solver - self.solver.setup(suppress_output=True, method='LSODA', step=1) + self.solver.setup(suppress_output=True, method='LSODA', step=1) def test_solver_initialization(self): """ Test the initialization of the solver. - This test verifies that the solver is correctly initialized as an instance of the + This test verifies that the solver is correctly initialized as an instance of the Solver class and that the model is correctly assigned to the solver. Assertions: @@ -129,8 +129,8 @@ def test_solver_setup(self): # Verify the generated functions self.assertIsNotNone(self.solver.pv_dfdt_global) self.assertIsNotNone(self.solver.s_u_update) - self.assertIsNotNone(self.solver.optimize) - self.assertIsNotNone(self.solver.initialize_by_function) + self.assertIsNotNone(self.solver.optimize) + self.assertIsNotNone(self.solver.initialize_by_function) def test_initialize_by_function(self): @@ -184,10 +184,10 @@ def test_optimize(self): 2. Verify that the optimize() method can run with the expected input. Expected input: - - y_temp: Temporary y values for optimization. Called y_temp to reflect the + - y_temp: Temporary y values for optimization. Called y_temp to reflect the name of the variable in the optimize() method. - keys4: Keys required for optimization. Called keys4 to reflect the name of the - variable in the optimize() method. + variable in the optimize() method. Raises: AssertionError: If the optimize() method does not accept the expected input @@ -244,7 +244,7 @@ def test_pv_dfdt_update(self): y0 = np.load(input_file_path) # Verify the function can run with the expected input - pv_dfdt_result = self.solver.pv_dfdt_global(t=0, y=y0) + pv_dfdt_result = self.solver.pv_dfdt_global(t=0, y=y0) # Verify the output matches the expected output @@ -313,11 +313,11 @@ def test_solver_solve(self): 1. Configures the solver with the given step size. 2. Solves the system using the solver. 3. Verifies that the solver has converged. - 7. Compares the new solution dictionary values with the expected values and asserts that they are + 7. Compares the new solution dictionary values with the expected values and asserts that they are within an acceptable tolerance. Raises: - AssertionError: If the solver did not converge or if the computed values do not match the expected + AssertionError: If the solver did not converge or if the computed values do not match the expected values within the tolerance. """ @@ -343,10 +343,10 @@ def test_solver_solve(self): self.parobj = KorakianitisMixedModel_parameters() # Initialize the model - self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, - parobj=self.parobj, + self.model = KorakianitisMixedModel(time_setup_dict=self.time_setup_dict, + parobj=self.parobj, suppress_printing=True) - + # Initialize the solver self.solver = Solver(model=self.model) @@ -354,10 +354,10 @@ def test_solver_solve(self): self.solver.setup(suppress_output=True, method='LSODA', step=i_cycle_step_size) # Running the model - self.solver.solve() + self.solver.solve() # Verify the solver converged - self.assertTrue(self.solver.converged or self.solver._Nconv is not None) + self.assertTrue(self.solver.converged or self.solver._Nconv is not None) # Redefine tind based on how many heart cycle have actually been necessary to reach steady state tind_fin = np.arange(start=self.model.time_object.n_t-self.model.time_object.n_c, @@ -372,7 +372,7 @@ def test_solver_solve(self): 'P_i': value.P_i.values[tind_fin].mean(), 'Q_i': value.Q_i.values[tind_fin].mean() } - + # Check that the values are the same as the expected values expected_ndarray = np.array( [self.expected_values["results"][str(i_cycle_step_size)][key1][key2] for key1 in new_dict.keys() for key2 in new_dict[key1].keys()] @@ -383,9 +383,9 @@ def test_solver_solve(self): np.abs((expected_ndarray - new_ndarray) / expected_ndarray), np.abs((expected_ndarray - new_ndarray)) ) - self.assertTrue((test_ndarray < RELATIVE_TOLERANCE).all(), + self.assertTrue((test_ndarray < RELATIVE_TOLERANCE).all(), f"Test failed for step size {i_cycle_step_size}: {test_ndarray}") if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main() From 0dc2a1f69cd8d3ca909cdd4d88b795f747408431 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 14 Apr 2025 16:17:12 +0100 Subject: [PATCH 08/10] fixed pre-commit yaml --- .pre-commit-config.yaml | 31 +++++++++++++++++++++++++++++-- 1 file changed, 29 insertions(+), 2 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 8fad403..a22db30 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,7 +1,10 @@ -# .pre-commit-config.yaml +ci: + autoupdate_commit_msg: "chore: update pre-commit hooks" + autofix_commit_msg: "style: pre-commit fixes" + repos: - repo: https://github.com/pre-commit/pre-commit-hooks - rev: v4.0.1 + rev: v5.0.0 hooks: - id: trailing-whitespace - id: end-of-file-fixer @@ -9,6 +12,11 @@ repos: - id: check-json - id: check-added-large-files args: ['--maxkb=500'] # set the max file size to 500KB + - id: check-case-conflict + - id: check-merge-conflict + - id: check-symlinks + # - id: debug-statements + - id: mixed-line-ending # - repo: https://github.com/psf/black # rev: 25.1.0 # Use the latest stable version @@ -24,3 +32,22 @@ repos: # rev: 6.0.0 # Use the latest stable version # hooks: # - id: isort + +# - repo: https://github.com/astral-sh/ruff-pre-commit +# rev: "v0.11.5" +# hooks: +# # first, lint + autofix +# - id: ruff +# types_or: [python, pyi, jupyter] +# args: ["--fix", "--show-fixes"] +# # then, format +# - id: ruff-format + +# - repo: https://github.com/pre-commit/mirrors-mypy +# rev: "v1.15.0" +# hooks: +# - id: mypy +# files: src +# args: [] +# additional_dependencies: +# - pytest From d1b59cc3938614a00479491f60ce47c9de501d7b Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 14 Apr 2025 16:37:44 +0100 Subject: [PATCH 09/10] addint CONTRIBUTING.md --- CONTRIBUTING.md | 47 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 47 insertions(+) create mode 100644 CONTRIBUTING.md diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 0000000..ceb4a73 --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,47 @@ +See the [Scientific Python Developer Guide][spc-dev-intro] for a detailed +description of best practices for developing scientific packages. + +[spc-dev-intro]: https://learn.scientific-python.org/development/ + +# Setting up a development environment manually + +You can set up a development environment by running: + +```zsh +python3 -m venv venv # create a virtualenv called venv +source ./venv/bin/activate # now `python` points to the virtualenv python +pip install -v -e ".[dev]" # -v for verbose, -e for editable, [dev] for dev dependencies +``` + +# Post setup + +You should prepare pre-commit, which will help you by checking that commits pass +required checks: + +```bash +pip install pre-commit # or brew install pre-commit on macOS +pre-commit install # this will install a pre-commit hook into the git repo +``` + +You can also/alternatively run `pre-commit run` (changes only) or +`pre-commit run --all-files` to check even without installing the hook. + +# Testing + +This repo uses `unittest` for testing. You can run locally the tests by running the following command: + +```bash + python -m unittest discover -s tests +``` +there is also a autamtated test pipeline that runs the tests on every push to the repository (see [here](.github/workflows/ci.yml)). + + +# Pre-commit + +This project uses pre-commit for all style checking. Install pre-commit and run: + +```bash +pre-commit run -a +``` + +to check all files. From c70965c2d98f4c9acd60f8e00dd4e5a43de7cc88 Mon Sep 17 00:00:00 2001 From: Levan Bokeria Date: Mon, 14 Apr 2025 16:48:34 +0100 Subject: [PATCH 10/10] explaining the pre-commit hooks --- .pre-commit-config.yaml | 21 ++++++++++++--------- 1 file changed, 12 insertions(+), 9 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index a22db30..7e521e3 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -6,17 +6,20 @@ repos: - repo: https://github.com/pre-commit/pre-commit-hooks rev: v5.0.0 hooks: - - id: trailing-whitespace - - id: end-of-file-fixer - - id: check-yaml - - id: check-json - - id: check-added-large-files + - id: trailing-whitespace # remove trailing whitespace + - id: end-of-file-fixer # ensure files end with a newline + - id: check-yaml # check YAML files for syntax errors + - id: check-json # check JSON files for syntax errors + - id: check-added-large-files # check for large files args: ['--maxkb=500'] # set the max file size to 500KB - - id: check-case-conflict - - id: check-merge-conflict - - id: check-symlinks + - id: check-case-conflict # check for case conflicts in filenames. + - id: check-merge-conflict # This hook checks for merge conflict markers in files. + # It ensures that there are no unresolved merge conflicts in the codebase. + - id: check-symlinks # check for broken symlinks # - id: debug-statements - - id: mixed-line-ending + - id: mixed-line-ending # check for mixed line endings, meaning that + # a file contains both CRLF and LF line endings. This can cause issues + # when working with files across different operating systems. # - repo: https://github.com/psf/black # rev: 25.1.0 # Use the latest stable version