PBPK_public/pythonScripts/primerjava modela IVP z modelom Matrix solve.py

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