Upload files to 'pythonScripts'

pythonScripts/primerjava modela IVP z modelom Matrix solve.py

@@ -0,0 +1,220 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import os
+import cModel
+import time
+import pdfkit
+import tempfile
+import runSolver
+import io
+import importlib
+# Primerjava obeh modelov in razlika obeh modelov
+# Če je venousInput":{"function":"venousInput"}, potem modela IVP pa Matrix_solve ne bosta enaka, ker je fora v A = P**T D P in to če je D= D(t) ne bo delovalo
+#run solver                          
+fh=os.path.expanduser('~')                      #Nastavitev poti in konfiguracijskih datotek
+
+importlib.reload(cModel)
+importlib.reload(runSolver)
+
+def getModel(solution, modelName):
+    """
+    Load and parse the model from the given solution dictionary using the modelName.
+    """
+    Q = solution[modelName]  # Retrieve the solution dictionary for the specified model
+    model = cModel.model()  # Create an instance of the cModel class
+    setupFile = Q['setup']  # Get the setup file path from the solution
+    modelFile = Q['model']  # Get the model file path from the solution
+    parameterFile = Q['parameters']  # Get the parameter file path from the solution
+    print('modelFile: {} {}'.format(modelFile, os.path.isfile(modelFile)))  # Print the model file path and its existence status
+    print('parameterFile: {} {}'.format(parameterFile, os.path.isfile(parameterFile)))  # Print the parameter file path and its existence status
+    model.parse(modelFile, parameterFile)  # Parse the model file with the given parameters
+    return model  # Return the loaded model
+
+def mergeSolutions(seq):
+    """
+    Merge multiple solution dictionaries into a single dictionary.
+    """
+    out = {}
+    for v in ['t', 'sol', 'se', 'qt', 'sOut']:  # Merge arrays from each solution
+        out[v] = np.concatenate([x[v] for x in seq])
+    for v in ['lut', 'lutSE', 'setup', 'model', 'parameters', 'qt', 'sOut']:  # Use data from the last solution for these keys
+        out[v] = seq[-1][v]
+    return out  # Return the merged dictionary
+
+def plot_results(solution, ivp_model_name, matrix_model_name, epsilon=1e-10):
+    """
+    Plot results for the specified IVP and matrix models, including their difference.
+    """
+    ivp_model = getModel(solution, ivp_model_name)
+    matrix_model = getModel(solution, matrix_model_name)
+    
+    ivp_setup = solution[ivp_model_name]['setup']
+    matrix_setup = solution[matrix_model_name]['setup']
+    
+    ivp_tscale = runSolver.getScale(ivp_setup)
+    matrix_tscale = runSolver.getScale(matrix_setup)
+    
+    print(f'IVP tscale={ivp_tscale}')
+    print(f'Matrix tscale={matrix_tscale}')
+    
+    name = ['venous', 'adipose', 'brain', 'heart', 'kidney', 'liver', 'lung', 'muscle', 'skin', 'stomach', 'splanchnic', 'excrement']
+    tmax = solution[ivp_model_name]['t'][-1]
+    
+    compartment_max = [-1] * len(name)  # Renamed variable to avoid conflict with built-in max function
+
+    models = {
+        ivp_model_name: {'color': 'yellow', 'shadeColor': 'gold'},
+        matrix_model_name: {'color': 'blue', 'shadeColor': 'lightgreen', 'linestyle': '--'},  
+        'Matrix-IVP': {'color': 'red', 'shadeColor': 'pink'}
+    }
+
+    fig, axs = plt.subplots(4, 3, figsize=(15, 20))
+    
+    for i in range(len(name)):
+        row = i // 3
+        col = i % 3
+        ax = axs[row, col]
+        
+        for model_name in [ivp_model_name, matrix_model_name]:
+            model = solution[model_name]
+            v = name[i]
+            
+            try:
+                j = model['lut'][v]
+            except KeyError:
+                try:
+                    v1 = alias[v]
+                    j = model['lut'][v1]
+                except KeyError:
+                    print(f'No data for {v}/{v1}')
+                    continue
+
+            fy = model['sol'][:, j]
+            fe = model['se'][:, j]
+            t = model['t']
+
+            ax.plot(
+                t / (ivp_tscale if model_name == ivp_model_name else matrix_tscale),
+                fy,
+                color=models[model_name]['color'],
+                linestyle=models[model_name].get('linestyle', '-'),  # Dodajmo slog črte za matriko
+                label=f'{model_name} Mean'
+            )
+            ax.fill_between(
+                t / (ivp_tscale if model_name == ivp_model_name else matrix_tscale),
+                fy - fe, fy + fe,
+                color=models[model_name]['shadeColor'],
+                alpha=0.1
+            )
+            ax.plot(
+                t / (ivp_tscale if model_name == ivp_model_name else matrix_tscale),
+                fy - fe,
+                color=models[model_name]['shadeColor'],
+                linewidth=1,
+                alpha=0.2
+            )
+            ax.plot(
+                t / (ivp_tscale if model_name == ivp_model_name else matrix_tscale),
+                fy + fe,
+                color=models[model_name]['shadeColor'],
+                linewidth=1,
+                alpha=0.2
+            )
+
+        # Calculate and plot the difference of matrix and IVP
+        try:
+            ivp_fy = solution[ivp_model_name]['sol'][:, j]
+            matrix_fy = solution[matrix_model_name]['sol'][:, j]
+            difference = matrix_fy - ivp_fy
+            difference_fe = np.sqrt(solution[matrix_model_name]['se'][:, j]**2 + solution[ivp_model_name]['se'][:, j]**2)
+            ax.plot(
+                t / ivp_tscale,
+                difference,
+                color=models['Matrix-IVP']['color'],
+                label='Matrix - IVP Difference'
+            )
+            ax.fill_between(
+                t / ivp_tscale,
+                difference - difference_fe,
+                difference + difference_fe,
+                color=models['Matrix-IVP']['shadeColor'],
+                alpha=0.1
+            )
+            ax.plot(
+                t / ivp_tscale,
+                difference - difference_fe,
+                color=models['Matrix-IVP']['shadeColor'],
+                linewidth=1,
+                alpha=0.2
+            )
+            ax.plot(
+                t / ivp_tscale,
+                difference + difference_fe,
+                color=models['Matrix-IVP']['shadeColor'],
+                linewidth=1,
+                alpha=0.2
+            )
+        except KeyError:
+            print(f'No difference data for {name[i]}')
+
+        if compartment_max[i] > 0:
+            ax.set_ylim([0, compartment_max[i]])
+        ax.set_xlim([0, 1.1 * tmax / ivp_tscale])
+        ax.set_xlabel(ivp_setup['tUnit'])
+        ax.set_title(name[i])
+        ax.legend()
+
+    output_file = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'results', 'combined_plot_matrix_ivp.png')
+    plt.savefig(output_file)
+    plt.show()
+
+
+
+def main():
+    """
+    Main function to load data, analyze model, and plot results.
+    """
+    fh = os.path.expanduser('~')
+
+    # Define job configurations for loading
+    job_configs = {
+        'cDiazepam_IVP': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam_IVP')
+        },
+        'cDiazepam1_IVP': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam1_IVP')
+        },
+        'cDiazepamF_IVP': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamF_IVP')
+        },
+        'cDiazepamB_IVP': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamB_IVP')
+        },
+        'cDiazepam_Matrix': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam_Matrix')
+        },
+        'cDiazepam1_Matrix': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam1_Matrix')
+        },
+        'cDiazepamF_Matrix': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamF_Matrix')
+        },
+        'cDiazepamB_Matrix': {
+            'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamB_Matrix')
+        }
+    }
+
+    # Load and merge solutions
+    solution = {}
+    for modelName, job in job_configs.items():
+        seq = [runSolver.loadSolutionFromDir(job['jobDir'], True)]
+        solution[modelName] = mergeSolutions(seq)
+
+    # Analyze and plot results
+    for ivp_model_name in [name for name in job_configs if 'IVP' in name]:
+        matrix_model_name = ivp_model_name.replace('IVP', 'Matrix')
+        if matrix_model_name in job_configs:
+            plot_results(solution,ivp_model_name, matrix_model_name)
+
+if __name__ == "__main__":
+    main()