risanje grafov_IVP.py 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. import os
  4. import cModel
  5. import time
  6. import pdfkit
  7. import tempfile
  8. import runSolver
  9. import io
  10. import importlib
  11. #Nariše graf iz IVP modelov
  12. importlib.reload(cModel)
  13. importlib.reload(runSolver)
  14. def getModel(solution, modelName):
  15. """
  16. Load and parse the model from the given solution dictionary using the modelName.
  17. """
  18. Q = solution[modelName] # Retrieve the solution dictionary for the specified model
  19. model = cModel.model() # Create an instance of the cModel class
  20. setupFile = Q['setup'] # Get the setup file path from the solution
  21. modelFile = Q['model'] # Get the model file path from the solution
  22. parameterFile = Q['parameters'] # Get the parameter file path from the solution
  23. print('modelFile: {} {}'.format(modelFile, os.path.isfile(modelFile))) # Print the model file path and its existence status
  24. print('parameterFile: {} {}'.format(parameterFile, os.path.isfile(parameterFile))) # Print the parameter file path and its existence status
  25. model.parse(modelFile, parameterFile) # Parse the model file with the given parameters
  26. return model # Return the loaded model
  27. def mergeSolutions(seq):
  28. """
  29. Merge multiple solution dictionaries into a single dictionary.
  30. """
  31. out = {}
  32. for v in ['t', 'sol', 'se', 'qt', 'sOut']: # Merge arrays from each solution
  33. out[v] = np.concatenate([x[v] for x in seq])
  34. for v in ['lut', 'lutSE', 'setup', 'model', 'parameters', 'qt', 'sOut']: # Use data from the last solution for these keys
  35. out[v] = seq[-1][v]
  36. return out # Return the merged dictionary
  37. # Define paths
  38. fh = os.path.expanduser('~')
  39. # Define job configurations for loading
  40. job_configs = {
  41. 'cDiazepam_IVP': {
  42. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam_IVP')
  43. },
  44. 'cDiazepam1_IVP': {
  45. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam1_IVP')
  46. },
  47. 'cDiazepamF_IVP': {
  48. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamF_IVP')
  49. },
  50. 'cDiazepamB_IVP': {
  51. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamB_IVP')
  52. }
  53. }
  54. # Load and merge solutions
  55. solution = {}
  56. for modelName, job in job_configs.items():
  57. seq = [runSolver.loadSolutionFromDir(job['jobDir'], True)]
  58. solution[modelName] = mergeSolutions(seq)
  59. def analyze_model(model, setup):
  60. """
  61. Perform matrix operations and print results.
  62. """
  63. tscale = runSolver.getScale(setup) # Get the time scale from the setup
  64. print(f'tscale={tscale}') # Print the time scale
  65. model.inspect() # Inspect the model
  66. print("***********done************")
  67. print(model.M(1).shape) # Print the shape of the matrix from the model
  68. print(model.m) # Print the model parameters
  69. nt = setup['nt'] # Number of time points
  70. qtmax = 1 # Maximum time (minute)
  71. qt = np.linspace(0, qtmax, nt) # Generate time points
  72. par = 'venousInput' # Parameter to plot
  73. # Plot parameters
  74. try:
  75. hw = model.get(par) # Get the parameter from the model
  76. print(hw) # Print the parameter information
  77. ht = [10 * hw['value'](x) for x in qt] # Calculate parameter values over time
  78. plt.plot(qt / tscale, ht) # Plot the parameter values
  79. except (KeyError, TypeError): # Handle errors if the parameter is not found
  80. print(f'Troubles getting {par}')
  81. pass
  82. # Measure performance of matrix operations
  83. start_time = time.time() # Record start time
  84. for i in range(100000): # Perform matrix operations 100,000 times
  85. model.M(1e7) # Call the matrix operation
  86. end_time = time.time() # Record end time
  87. print('Time: {:.3f} s'.format(end_time - start_time)) # Print elapsed time
  88. fM = model.M(0) # Get the matrix from the model
  89. print('Rank: {}/{}'.format(np.linalg.matrix_rank(fM), fM.shape)) # Print the rank of the matrix
  90. np.set_printoptions(suppress=True, precision=2, linewidth=150) # Set NumPy print options
  91. print(f'{fM}') # Print the matrix
  92. v, Q = np.linalg.eig(fM) # Perform eigenvalue decomposition
  93. np.set_printoptions(suppress=False) # Reset NumPy print options
  94. print(Q[:, 2:4]) # Print selected eigenvectors
  95. Q1 = np.linalg.inv(Q) # Calculate the inverse of the eigenvector matrix
  96. D = np.diag(v) # Create a diagonal matrix of eigenvalues
  97. t = np.linspace(0, 100, 101) # Generate time points for plotting
  98. fy = [model.u(x)[14] for x in t] # Calculate model output over time
  99. print('{} {}'.format(len(fy), len(t))) # Print lengths of output and time points
  100. plt.figure() # Create a new figure
  101. plt.plot(fy) # Plot the model output
  102. fu = model.u(0) # Get the model output at time 0
  103. print(Q1 @ fu) # Print the result of the matrix multiplication
  104. def plot_results(solution, modelName):
  105. """
  106. Plot results for the specified model from the solution.
  107. """
  108. model = getModel(solution, modelName) # Load the model
  109. setup = solution[modelName]['setup'] # Get the setup from the solution
  110. tscale = runSolver.getScale(setup) # Get the time scale
  111. print(f'tscale={tscale}') # Print the time scale
  112. model.inspect() # Inspect the model
  113. name = ['venous', 'adipose', 'brain', 'heart', 'kidney', 'liver', 'lung', 'muscle', 'skin', 'stomach', 'splanchnic', 'excrement'] # Compartment names
  114. tmax = solution[modelName]['t'][-1] # Get the maximum time from the solution
  115. # Define colors and shading for models
  116. models = {
  117. 'cDiazepam_IVP': {'color': 'orange', 'shadeColor': 'red'},
  118. 'cDiazepamF_IVP': {'color': 'blue', 'shadeColor': 'green'},
  119. 'cDiazepam1_IVP': {'color': 'black', 'shadeColor': 'yellow'},
  120. 'cDiazepamB_IVP': {'color': 'purple', 'shadeColor': 'gold'}
  121. }
  122. fig, axs = plt.subplots(4, 3, figsize=(15, 20)) # Create a grid of subplots
  123. for i in range(len(name)): # Loop over each compartment
  124. row = i // 3 # Determine row index
  125. col = i % 3 # Determine column index
  126. ax = axs[row, col] # Get the subplot
  127. for m in models: # Loop over each model
  128. fM = models[m] # Get model configuration
  129. Q = solution[m] # Get solution for the model
  130. v = name[i] # Get the compartment name
  131. try:
  132. j = Q['lut'][v] # Get the index for the compartment
  133. except KeyError:
  134. try:
  135. v1 = alias[v] # Handle alias if needed
  136. j = Q['lut'][v1] # Get the index for the alias
  137. except KeyError:
  138. print(f'No data for {v}/{v1}') # Print error if compartment not found
  139. continue
  140. fy = Q['sol'][:, j] # Get the solution for the compartment
  141. fe = Q['se'][:, j] # Get the standard error for the compartment
  142. t = Q['t'] # Get the time points
  143. # Plot the mean values and confidence intervals
  144. ax.plot(t / tscale, fy, color=fM['color'], label=f'{m} Mean')
  145. ax.fill_between(t / tscale, fy - fe, fy + fe, color=fM['shadeColor'], alpha=0.1) # Shaded area for confidence intervals
  146. ax.plot(t / tscale, fy - fe, color=fM['shadeColor'], linewidth=1, alpha=0.2) # Lower bound
  147. ax.plot(t / tscale, fy + fe, color=fM['shadeColor'], linewidth=1, alpha=0.2) # Upper bound
  148. # Find the maximum y-value including confidence intervals
  149. max_y_value = max(np.max(fy + fe), np.max(fy - fe))
  150. ax.set_ylim([0, min(max_y_value, 10000)]) # Set y-axis limit to max value or 5500, whichever is smaller
  151. ax.set_xlim([0, 1.1 * tmax / tscale]) # Set x-axis limits
  152. ax.set_xlabel(setup['tUnit']) # Set x-axis label
  153. ax.set_title(name[i]) # Set plot title
  154. ax.legend() # Add legend
  155. output_file = os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'results', 'plot.png') # Define output file path
  156. plt.savefig(output_file) # Save the plot to file
  157. plt.show() # Display the plot
  158. def main():
  159. """
  160. Main function to load data, analyze model, and plot results.
  161. """
  162. fh = os.path.expanduser('~') # Define file path root
  163. # Define job configurations for loading
  164. job_configs = {
  165. 'cDiazepam_IVP': {
  166. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam_IVP')
  167. },
  168. 'cDiazepam1_IVP': {
  169. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepam1_IVP')
  170. },
  171. 'cDiazepamF_IVP': {
  172. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamF_IVP')
  173. },
  174. 'cDiazepamB_IVP': {
  175. 'jobDir': os.path.join(fh, 'Documents', 'Sola', 'IJS', 'PBPK_public', 'cDiazepamB_IVP')
  176. }
  177. }
  178. # Load and merge solutions
  179. solution = {}
  180. for modelName, job in job_configs.items():
  181. seq = [runSolver.loadSolutionFromDir(job['jobDir'], True)] # Load solution from directory
  182. solution[modelName] = mergeSolutions(seq) # Merge solutions
  183. # Capture print output
  184. output_stream = io.StringIO()
  185. analyze_model(getModel(solution, 'cDiazepamF_IVP'), solution['cDiazepamF_IVP']['setup']) # Analyze the model
  186. analysis_results = analyze_model(getModel(solution, 'cDiazepamF_IVP'), solution['cDiazepamF_IVP']['setup'])
  187. plot_results(solution, 'cDiazepamF_IVP') # Plot the results
  188. if __name__ == "__main__":
  189. main() # Execute the main function