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()